package hep.aida.ref.pdf;

/**
 *
 * @author The FreeHEP team @ SLAC
 *
 */
public abstract class FunctionDerivative {
    
    /**
     * Computes the derivative of a function with respect to one argument.
     * Adapted from Numerical Recipes in FORTRAN dfridr, pp. 182-183.
     * @param f a JAIDA IFunction
     * @param x a double[] with the point where the derivative is evaluated
     * @param ind the index of the argument of f we differentiate wrto
     * @param h a typical step size for x[ind]
     * @return the value of df/dx[ind] at x[]
     */
    public static double derivative(Function f, Variable x, double h) {
        final double BIG  = Double.MAX_VALUE;
        final double CON  = 1.4;
        final double CON2 = CON*CON;
        final double SAFE = 2.;
        
        final int NTAB = 10;
        double[][] a = new double[NTAB][NTAB];
        
        double res = Double.NaN;
        double err = BIG;
                
        a[0][0] = evaluateMatrixElement(f,x,h);

        for (int i=1; i= SAFE*err) {
                return res;
            }
        }
        return res;
    }
    
    private static double evaluateMatrixElement(Function f, Variable x, double increment) {
        double xVal = x.value();
        
        double result = 0;
        
        double xPlus = xVal + increment;
        double xMinus = xVal - increment;
        
        x.setValue(xPlus);
        result += f.value();
        
        x.setValue(xMinus);
        result -= f.value();

        result /= 2*increment;
        
        // Set the dependents back to its original value
        x.setValue(xVal);
        return result;
    }
    
    
}