using System;

namespace Science.Mathematics.VectorCalculus
{
    public class Differentiation
    {
        private Function.ToLastType<double,double> del;
        private double at1, result;
        public Function.ToLastType<double, double> Function
        {
            set { del = value; }
        }
        public double At
        {
            set { at1 = value; }
        }
        public double Result
        {
            get { return result; }
        }		
       
        public Differentiation(Function.ToLastType<double, double> f, double at)
		{
			del = f;
            at1 = at;
		}
	
		public void Compute()
		{
            RiddersMethod df = new RiddersMethod();
            result = df.Find(del, at1, 0.001);       	
		}       
    }

    class RiddersMethod
    {
        private double err;
        public double Error
        {
            get { return err; }
        }

        public RiddersMethod()
        {
        }
        
        private double CON = 1.4;
        private double CON2 = 1.4 * 1.4;
        private double BIG = 1.0E+30;
        private int NTAB = 10;
        private double SAFE = 2.0;

        public double Find(Function.ToLastType<double, double> func, double x, double h)
        {
            int i, j;
            double errt, fac, hh, ans = 0.0;
            double[,] a = new double[NTAB, NTAB];

            hh = h;
            a[0, 0] = (func(x + hh) - func(x - hh)) / (2.0 * hh);
            err = BIG;
            for (i = 1; i < NTAB; i++)
            {
                hh /= CON;
                a[0, i] = (func(x + hh) - func(x - hh)) / (2.0 * hh);
                fac = CON2;
                for (j = 1; j <= i; j++)
                {
                    a[j, i] = (a[j - 1, i] * fac - a[j - 1, i - 1]) / (fac - 1.0);
                    fac = CON2 * fac;
                    errt = Math.Max(Math.Abs(a[j, i] - a[j - 1, i]), Math.Abs(a[j, i] - a[j - 1, i - 1]));
                    if (errt <= err)
                    {
                        err = errt;
                        ans = a[j, i];
                    }
                }
                if (Math.Abs(a[i, i] - a[i - 1, i - 1]) >= SAFE * err)
                {
                    return ans;
                }
            }
            return ans;
        }  
    }
}