using System;

namespace Science.Mathematics.VectorCalculus
{
    public class IntegrationMidpoint
	{
        private Function.ToLastType<double,double> del;
        private double from1, to1, result;
        public IntegrationMidpoint()
        {           
        }
        public IntegrationMidpoint(Function.ToLastType<double, double> f)
		{
			del = f;			
		}
        public IntegrationMidpoint(Function.ToLastType<double, double> f, double from, double to)
        {
            del = f;
            from1 = from;
            to1 = to;
        }
		public void Compute()
		{
            Midpnt obj = new Midpnt(del, from1, to1);
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            result = obj.next();
        }

        public Function.ToLastType<double, double> Function
        {
            set{ del = value; }
        }
	    public double From
		{
			set{from1 = value;}
		}
		public double To
		{
			set{to1 = value;}
		}
		public double Result
		{
			get{return result;}
		}     
	}

    class Midpnt : Quadrature
    {
        protected double a, b, s = 0.0;
        protected Function.ToLastType<double, double> funk;
        public virtual double func(double x)
        {
            return funk(x);
        }
        public Midpnt(Function.ToLastType<double, double> funcc, double aa, double bb)
        {
            funk = funcc;
            a = aa;
            b = bb;
            n = 0;
        }
        public override double next()
        {
            int it, j;
            double x, tnm, sum, del, ddel;
            n++;
            if (n == 1)
            {
                return (s = (b - a) * func(0.5 * (a + b)));
            }
            else
            {
                for (it = 1, j = 1; j < n - 1; j++) it *= 3;
                tnm = it;
                del = (b - a) / (3.0 * tnm);
                ddel = del + del;
                x = a + 0.5 * del;
                sum = 0.0;
                for (j = 0; j < it; j++)
                {
                    sum += func(x);
                    x += ddel;
                    sum += func(x);
                    x += del;
                }
                s = (s + (b - a) * sum / tnm) / 3.0;
                return s;
            }
        }
    }
}