using System;

namespace Science.Mathematics.VectorCalculus
{
    public class Surface
    {
        private Function.ToLastType<double[],double[]> g;
        public Surface()
        {
        }
        public Surface(Function.ToLastType<double[], double[]> parametrization)
        {
            g = parametrization;
        }
        public Function.ToLastType<double[], double[]> Parametrization
        {
            get { return g; }
            set { g = value; }
        }

        private double start1, end1;
        private double start2, end2;
        public double Parameter1StartValue
        {
            get { return start1; }
            set { start1 = value; }
        }
        public double Parameter1EndValue
        {
            get { return end1; }
            set { end1 = value; }
        }
        public double Parameter2StartValue
        {
            get { return start2; }
            set { start2 = value; }
        }
        public double Parameter2EndValue
        {
            get { return end2; }
            set { end2 = value; }
        }

        public double Area
        {
            get{return FindArea();}
        }

        private double FindArea()
        {
            double[] f = {start1, start2};
            double[] t = {end1, end2};
            Function.ToLastType<double[], double> norm
                = new Function.ToLastType<double[], double>(Func);
            IntegrationMultiD obj = new IntegrationMultiD(norm, f, t);
            obj.Compute();            
            return obj.BestEstimation;
        }
        private double Func(double[] u)
        {
            TangentVectorOnSurface tv 
                = new TangentVectorOnSurface(this,u);
            tv.Compute();
            Vector n = tv.Result[0] % tv.Result[1];
            return n.Norm;
        }

        public double CurvatureGauss(double[] at)
        {
            FindCurvature(at);
            return (l * n - m * m) / W;
        }

        public double CurvatureMean(double[] at)
        {
            FindCurvature(at);
            return (G*l + E*n - 2.0*F* m) / 2.0/ W;
        }

        private double W, E, F, G, l, m, n;
        private int i = 0;
        private void FindCurvature(double[] u)
        {
            TangentVectorOnSurface tv
                = new TangentVectorOnSurface(this, u);
            tv.Compute();
            E = tv.Result[0] * tv.Result[0];
            F = tv.Result[0] * tv.Result[1];
            G = tv.Result[1] * tv.Result[1];
            Vector v = tv.Result[0] % tv.Result[1];
            W = v * v;
            UnitVector N = new UnitVector(v);

            Function.ToLastType<double[], double> one
                = new Function.ToLastType<double[], double>(gg);
            int dim = g(u).Length;
            double sum00 = 0.0, sum01 = 0.0, sum11 = 0.0;            
            for (int j = 0; j < dim; j++)
            {
                i = j;
                IteratedPartialDerivative obj = new IteratedPartialDerivative(one, u, 2);
                obj.WithRespectToFirst = 0;
                obj.WithRespectToSecond = 0;
                obj.Compute();
                sum00 += N[j] * obj.Result;
                obj.WithRespectToFirst = 0;
                obj.WithRespectToSecond = 1;
                obj.Compute();
                sum01 += N[j] * obj.Result;
                obj.WithRespectToFirst = 1;
                obj.WithRespectToSecond = 1;
                obj.Compute();
                sum11 += N[j] * obj.Result;
            }
            l = sum00;
            m = sum01;
            n = sum11;
        }
        private double gg(double[] u)
        {
            return g(u)[i];
        }
    }
}