using System;

namespace Science.Physics.GeneralPhysics
{
	/// <summary>
	/// VectorCalculus contains many static methods:
	/// Gradient, Divergence, Curl, LineIntegral, SurfaceIntegral, VolumeIntegral.
    /// We consider only three-dimensional vectors.
	/// </summary>
	public class VectorCalculus
	{
		public VectorCalculus()
		{
		}
        private static Vector.Field f;
        private static Vector.FunctionOfTime ft;
        private static Vector.FunctionOfPosition fp;
        private static Scalar.Field g;
        private static Scalar.FunctionOfTime gt;
        private static Scalar.FunctionOfPosition gp;

        private static Line.Parameterization ptp;
        private static Surface.Parameterization p2tp;
        private static Volume.Parameterization p3tp;

        private static Position p = new Position();
        private static Time tm = new Time();

        /// <summary>
        /// Divergence for a vector defined by a function of position.
        /// </summary> 
        public static Scalar Divergence(Vector v, Position x)
        {
            fp = v.VectorFunctionOfPosition;
            double[] pp = { x.X, x.Y, x.Z };

            Scalar res = new Scalar();
            double sum = 0.0;

            sum += Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpx), pp, 0);
            sum += Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpy), pp, 1);
            sum += Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpz), pp, 2);

            res.Magnitude = sum;
            return res;
        }
        private static double fpx(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return fp(p).X;
        }
        private static double fpy(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return fp(p).Y;
        }
        private static double fpz(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return fp(p).Z;
        }

        /// <summary>
        /// Divergence for a vector defined by a field.
        /// </summary> 
		public static Scalar Divergence(Vector v, Position x, Time t)
		{
		    f = v.VectorField;
            tm.s = t.s;
            double[] pp = { x.X, x.Y, x.Z };

            Scalar res = new Scalar();
            double sum = 0.0;

            sum += Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fx), pp, 0);
            sum += Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fy), pp, 1);
            sum += Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fz), pp, 2);

            res.Magnitude = sum;
            return res;
		}
        private static double fx(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return f(p, tm).X;
        }
        private static double fy(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return f(p, tm).Y;
        }
        private static double fz(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return f(p, tm).Z;
        }

        /// <summary>
        /// Curl for a vector defined by a function of position.
        /// </summary>      
        public static Vector Curl(Vector v, Position x)
        {
            fp = v.VectorFunctionOfPosition;
            double[] pp = { x.X, x.Y, x.Z };

            double[,] res = new Double[3, 3];

            res[1, 0] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpx), pp, 1);
            res[2, 0] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpx), pp, 2);

            res[0, 1] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpy), pp, 0);
            res[2, 1] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpy), pp, 2);

            res[0, 2] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpz), pp, 0);
            res[1, 2] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fpz), pp, 1);
            
            Vector w = new Vector();
            w.X = res[1, 2] - res[2, 1];
            w.Y = res[2, 0] - res[0, 2];
            w.Z = res[0, 1] - res[1, 0];
            return w;
        }

        /// <summary>
        /// Curl for a vector defined by a field.
        /// </summary> 
		public static Vector Curl(Vector v, Position x, Time t)
		{
            f = v.VectorField;
            tm.s = t.s;
            double[] pp = { x.X, x.Y, x.Z };

			double[,] res = new Double[3,3];

            res[1, 0] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fx), pp, 1);
            res[2, 0] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fx), pp, 2);

            res[0, 1] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fy), pp, 0);
            res[2, 1] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fy), pp, 2);

            res[0, 2] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fz), pp, 0);
            res[1, 2] = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(fz), pp, 1);

			Vector w = new Vector();
            w.X = res[1,2]-res[2,1];
            w.Y = res[2,0]-res[0,2];
            w.Z = res[0,1]-res[1,0];
			return w;
		}

       
        /// <summary>
        /// Gradient for a scalar defined by a function of position.
        /// </summary> 
        public static Vector Gradient(Scalar s, Position x)
        {
            gp = s.ScalarFunctionOfPosition;

            Vector res = new Vector();
            double[] pp = { x.X, x.Y, x.Z };
            
            res.X = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(gps), pp, 0);
            res.Y = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(gps), pp, 1);
            res.Z = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(gps), pp, 2);

            return res;
        }
        private static double gps(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return gp(p).Magnitude;
        }

        /// <summary>
        /// Gradient for a scalar defined by a field.
        /// </summary>
        public static Vector Gradient(Scalar s, Position x, Time t)
        {
            g = s.ScalarField;
            tm.s = t.s;

            Vector res = new Vector();
            double[] pp = { x.X, x.Y, x.Z };

            res.X = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(gs), pp, 0);
            res.Y = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(gs), pp, 1);
            res.Z = Calculus.PartialDifferentiation(new Calculus.Function<double[], double>(gs), pp, 2);

            return res;
        }
        private static double gs(double[] x)
        {
            p.X = x[0];
            p.Y = x[1];
            p.Z = x[2];
            return g(p, tm).Magnitude;
        }


        /// <summary>
        /// Time Partial Derivative for a vector defined by a field.
        /// </summary>
        public static Vector TimePartialDerivative(Vector v, Position x, Time t)
        {
            f = v.VectorField;
            p.X = x.X;
            p.Y = x.Y;
            p.Z = x.Z;

            Vector res = new Vector();

            res.X = Calculus.Differentiation(new Calculus.Function<double, double>(fx1), t.s);
            res.Y = Calculus.Differentiation(new Calculus.Function<double, double>(fy1), t.s);
            res.Z = Calculus.Differentiation(new Calculus.Function<double, double>(fz1), t.s);

            return res;
        }
        private static double fx1(double t)
        {
            tm.s = t;
            return f(p, tm).X;
        }
        private static double fy1(double t)
        {
            tm.s = t;
            return f(p, tm).Y;
        }
        private static double fz1(double t)
        {
            tm.s = t;
            return f(p, tm).Z;
        }


        /// <summary>
        /// Time Derivative for a vector defined by a function of time.
        /// </summary>
		public static Vector TimeDerivative(Vector v, Time t)
		{
            ft = v.VectorFunctionOfTime;
      
            Vector res = new Vector();

            res.X = Calculus.Differentiation(new Calculus.Function<double, double>(ftx), t.s);
            res.Y = Calculus.Differentiation(new Calculus.Function<double, double>(fty), t.s);
            res.Z = Calculus.Differentiation(new Calculus.Function<double, double>(ftz), t.s);

            return res;
		}
        private static double ftx(double t)
        {
            tm.s = t;
            return ft(tm).X;
        }
        private static double fty(double t)
        {
            tm.s = t;
            return ft(tm).Y;
        }
        private static double ftz(double t)
        {
            tm.s = t;
            return ft(tm).Z;
        }


        /// <summary>
        /// Time Partial Derivative for a scalar defined by a field.
        /// </summary>
        public static Scalar TimePartialDerivative(Scalar s, Position x, Time t)
        {
            g = s.ScalarField;
            p.X = x.X;
            p.Y = x.Y;
            p.Z = x.Z;

            Scalar res = new Scalar();
            res.Magnitude = Calculus.Differentiation(new Calculus.Function<double, double>(gs1), t.s);

            return res;
        }
        private static double gs1(double t)
        {
            tm.s = t;
            return g(p, tm).Magnitude;
        }   
        

        /// <summary>
        /// Time Derivative for a scalar defined by a function of time.
        /// </summary>
        public static Scalar TimeDerivative(Scalar s, Time t)
        {
            gt = s.ScalarFunctionOfTime;

            Scalar res = new Scalar();
            res.Magnitude = Calculus.Differentiation(new Calculus.Function<double, double>(gts), t.s);

            return res;
        }
        private static double gts(double t)
        {
            tm.s = t;
            return gt(tm).Magnitude;
        }



        /// <summary>
        /// Time Integral for a vector defined by a function of time.
        /// </summary>
        public static Vector TimeIntegral(Vector v, Time ti, Time tf)
        {
            ft = v.VectorFunctionOfTime;

            Vector res = new Vector();

            res.X = Calculus.Integration1D(new Calculus.Function<double, double>(ftx), ti.s, tf.s);
            res.Y = Calculus.Integration1D(new Calculus.Function<double, double>(fty), ti.s, tf.s);
            res.Z = Calculus.Integration1D(new Calculus.Function<double, double>(ftz), ti.s, tf.s);

            return res;
        }

        /// <summary>
        /// Time Integral for a vector defined by a field.
        /// </summary>
        public static Vector TimeIntegral(Vector v, Position x, Time ti, Time tf)
        {
            f = v.VectorField;
            p.X = x.X;
            p.Y = x.Y;
            p.Z = x.Z;

            Vector res = new Vector();
        
            res.X = Calculus.Integration1D(new Calculus.Function<double, double>(fx1), ti.s, tf.s);
            res.Y = Calculus.Integration1D(new Calculus.Function<double, double>(fy1), ti.s, tf.s);
            res.Z = Calculus.Integration1D(new Calculus.Function<double, double>(fz1), ti.s, tf.s);

            return res;
        }

        /// <summary>
        /// Time Integral for a scalar defined by a function of time.
        /// </summary>
        public static Scalar TimeIntegral(Scalar s, Time ti, Time tf)
        {
            gt = s.ScalarFunctionOfTime;

            Scalar res = new Scalar();
            res.Magnitude = Calculus.Integration1D(new Calculus.Function<double, double>(gts), ti.s, tf.s);

            return res;
        }

        /// <summary>
        /// Time Integral for a scalar defined by a function of time.
        /// </summary>
        public static Scalar TimeIntegral(Scalar s, Position x, Time ti, Time tf)
        {
            g = s.ScalarField;
            p.X = x.X;
            p.Y = x.Y;
            p.Z = x.Z;

            Scalar res = new Scalar();
            res.Magnitude = Calculus.Integration1D(new Calculus.Function<double, double>(gs1), ti.s, tf.s);

            return res;
        }



        /// <summary>
        /// Line Integral for a scalar defined by a function of position.
        /// </summary>
        public static Vector LineIntegral(Scalar s, Line l)
        {
            gp = s.ScalarFunctionOfPosition;
            ptp = l.ParameterToPosition;
     
            Vector res = new Vector();

            res.X = Calculus.Integration1D(new Calculus.Function<double, double>(gpvx), l.ParameterStartValue, l.ParameterEndValue);
            res.Y = Calculus.Integration1D(new Calculus.Function<double, double>(gpvy), l.ParameterStartValue, l.ParameterEndValue);
            res.Z = Calculus.Integration1D(new Calculus.Function<double, double>(gpvz), l.ParameterStartValue, l.ParameterEndValue);
            
            return res;
        }
        private static double gpvx(double para)
        {
            Scalar ff = gp(ptp(para));            
            return ff.Magnitude * Calculus.Differentiation(new Calculus.Function<double, double>(givex), para);
        }
        private static double gpvy(double para)
        {
            Scalar ff = gp(ptp(para));           
            return ff.Magnitude * Calculus.Differentiation(new Calculus.Function<double, double>(givey), para);
        }
        private static double gpvz(double para)
        {
            Scalar ff = gp(ptp(para));          
            return ff.Magnitude * Calculus.Differentiation(new Calculus.Function<double, double>(givez), para);       
        } 
        private static double givex(double para)
        {
            Position x = ptp(para);
            return x.X;
        }
        private static double givey(double para)
        {
            Position x = ptp(para);
            return x.Y;
        }
        private static double givez(double para)
        {
            Position x = ptp(para);
            return x.Z;
        }


        /// <summary>
        /// Line Integral for a scalar defined by a field.
        /// </summary>
        public static Vector LineIntegral(Scalar s, Line l, Time t)
        {
            g = s.ScalarField;
            ptp = l.ParameterToPosition;
            tm.s = t.s;

            Vector res = new Vector();

            Calculus.Function<double, double> gvx1 = new Calculus.Function<double, double>(gvx);
            res.X = Calculus.Integration1D(gvx1, l.ParameterStartValue, l.ParameterEndValue);

            Calculus.Function<double, double> gvy1 = new Calculus.Function<double, double>(gvy);
            res.Y = Calculus.Integration1D(gvy1, l.ParameterStartValue, l.ParameterEndValue);

            Calculus.Function<double, double> gvz1 = new Calculus.Function<double, double>(gvz);
            res.Z = Calculus.Integration1D(gvz1, l.ParameterStartValue, l.ParameterEndValue);

            return res;
        }
        private static double gvx(double para)
        {
            Scalar ff = g(ptp(para), tm);
            Calculus.Function<double, double> givex1 = new Calculus.Function<double, double>(givex);
            return ff.Magnitude * Calculus.Differentiation(givex1, para);
        }
        private static double gvy(double para)
        {
            Scalar ff = g(ptp(para), tm);
            Calculus.Function<double, double> givey1 = new Calculus.Function<double, double>(givey);
            return ff.Magnitude * Calculus.Differentiation(givey1, para);
        }
        private static double gvz(double para)
        {
            Scalar ff = g(ptp(para), tm);
            Calculus.Function<double, double> givez1 = new Calculus.Function<double, double>(givez);
            return ff.Magnitude * Calculus.Differentiation(givez1, para);
        }


        /// <summary>
        /// Dot Line Integral for a vector defined by a function of position.
        /// </summary>
        public static Scalar LineIntegralDot(Vector v, Line l)
        {
            fp = v.VectorFunctionOfPosition;
            ptp = l.ParameterToPosition;
  
            Scalar res = new Scalar();
            Calculus.Function<double, double> fpdotv1 = new Calculus.Function<double, double>(fpdotv);
            res.Magnitude = Calculus.IntegrationMidpoint(fpdotv1,l.ParameterStartValue,l.ParameterEndValue);
            return res;
        }
        private static double fpdotv(double para)
        {
            Scalar ff = givev(para) * fp(ptp(para));
            return ff.Magnitude;
        }
        private static Vector givev(double para)
        {
            Vector v = new Vector();

            Calculus.Function<double, double> givex1 = new Calculus.Function<double, double>(givex);
            v.X = Calculus.Differentiation(givex1, para);

            Calculus.Function<double, double> givey1 = new Calculus.Function<double, double>(givey);
            v.Y = Calculus.Differentiation(givey1, para);

            Calculus.Function<double, double> givez1 = new Calculus.Function<double, double>(givez);
            v.Z = Calculus.Differentiation(givez1, para);

            return v;
        }


        /// <summary>
        /// Dot Line Integral for a vector defined by a field.
        /// </summary>
        public static Scalar LineIntegralDot(Vector v, Line l, Time t)
        {
            f = v.VectorField;
            ptp = l.ParameterToPosition;
            tm.s = t.s;

            Scalar res = new Scalar();
            Calculus.Function<double, double> fdotv1 = new Calculus.Function<double, double>(fdotv);
            res.Magnitude = Calculus.IntegrationMidpoint(fdotv1, l.ParameterStartValue, l.ParameterEndValue);
            return res;
        }
        private static double fdotv(double para)
        {
            Scalar ff = givev(para) * f(ptp(para),tm);
            return ff.Magnitude;
        }


        /// <summary>
        /// Cross Line Integral for a vector defined by a function of position. ( \int \vec a \times  d \vec l )
        /// </summary>       
        public static Vector LineIntegralCross(Vector v, Line l)
        {
            fp = v.VectorFunctionOfPosition;
            ptp = l.ParameterToPosition;
   
            Vector res = new Vector();
            res.X = Calculus.IntegrationMidpoint(new Calculus.Function<double, double>(vcrossfpx), 
                l.ParameterStartValue, l.ParameterEndValue);
            res.Y = Calculus.IntegrationMidpoint(new Calculus.Function<double, double>(vcrossfpy),
                l.ParameterStartValue, l.ParameterEndValue);
            res.Z = Calculus.IntegrationMidpoint(new Calculus.Function<double, double>(vcrossfpz),
                l.ParameterStartValue, l.ParameterEndValue);
            return res;
        }
        private static double vcrossfpx(double para)
        {
            Vector ff = fp(ptp(para)) % givev(para);
            return ff.X;
        } 
        private static double vcrossfpy(double para)
        {
            Vector ff = fp(ptp(para)) % givev(para);
            return ff.Y;
        }
        private static double vcrossfpz(double para)
        {
            Vector ff = fp(ptp(para)) % givev(para);
            return ff.Z;
        }


        /// <summary>
        /// Cross Line Integral for a vector defined by a field. 
        /// </summary>
        public static Vector LineIntegralCross(Vector v, Line l, Time t)
        {
            f = v.VectorField;
            ptp = l.ParameterToPosition;
            tm = t;

            Vector res = new Vector();

            res.X = Calculus.IntegrationMidpoint(new Calculus.Function<double, double>(vcrossfx),
                l.ParameterStartValue, l.ParameterEndValue);
            res.Y = Calculus.IntegrationMidpoint(new Calculus.Function<double, double>(vcrossfy),
                l.ParameterStartValue, l.ParameterEndValue);
            res.Z = Calculus.IntegrationMidpoint(new Calculus.Function<double, double>(vcrossfz),
                l.ParameterStartValue, l.ParameterEndValue);
            return res;
        }
        private static double vcrossfx(double para)
        {
            Vector ff = f(ptp(para), tm) % givev(para);
            return ff.X;
        }
        private static double vcrossfy(double para)
        {
            Vector ff = f(ptp(para), tm) % givev(para);
            return ff.Y;
        }
        private static double vcrossfz(double para)
        {
            Vector ff = f(ptp(para), tm) % givev(para);
            return ff.Z;
        }


        /// <summary>
        /// Surface Integral for a scalar defined by a function of position.
        /// </summary>
        public static Vector SurfaceIntegral(Scalar s, Surface sf)
        {
            gp = s.ScalarFunctionOfPosition;
            p2tp = sf.ParameterToPosition;
    
            Vector res = new Vector();
            
            double[] from = new double[2] { sf.Parameter1StartValue, sf.Parameter2StartValue };
            double[] to = new double[2] { sf.Parameter1EndValue, sf.Parameter2EndValue };
            
            double[] inx = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gpjx), from, to);
            res.X = inx[0];
            res.StandardDeviationOfX = inx[1];

            double[] iny = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gpjy), from, to);
            res.Y = iny[0];
            res.StandardDeviationOfY = iny[1];

            double[] inz = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gpjz), from, to);
            res.Z = inz[0];
            res.StandardDeviationOfZ = inz[1];

            return res;
        }
        private static double gx(double[] para)
        {
            Position x = p2tp(para[0], para[1]);
            return x.X;
        }
        private static double gy(double[] para)
        {
            Position x = p2tp(para[0], para[1]);
            return x.Y;
        }
        private static double gz(double[] para)
        {
            Position x = p2tp(para[0], para[1]);
            return x.Z;
        }
        private static double rd(double[] para, int xyz, int uv)   // xyz=0,1,2    uv = 0,1
        {
            Calculus.Function<double[], double> xyzofp2 = new Calculus.Function<double[], double>(gx);
            if (xyz == 1) xyzofp2 = new Calculus.Function<double[], double>(gy);
            if (xyz == 2) xyzofp2 = new Calculus.Function<double[], double>(gz);
            return Calculus.PartialDifferentiation(xyzofp2, para, uv);
        }
        private static double gpjx(double[] para)
        {
            return gp(p2tp(para[0], para[1])).Magnitude
                * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0));

        }
        private static double gpjy(double[] para)
        {
            return gp(p2tp(para[0], para[1])).Magnitude
                * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0));
        }
        private static double gpjz(double[] para)
        {
            return gp(p2tp(para[0], para[1])).Magnitude
                * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0));
        }


        /// <summary>
        /// Surface Integral for a scalar defined by a field.
        /// </summary>
        public static Vector SurfaceIntegral(Scalar s, Surface sf, Time t)
        {
            g = s.ScalarField;
            p2tp = sf.ParameterToPosition;
            tm.s = t.s;

            Vector res = new Vector();
            
            double[] from = new double[2]{sf.Parameter1StartValue, sf.Parameter2StartValue};
            double[] to = new double[2]{sf.Parameter1EndValue, sf.Parameter2EndValue};
            
            double[] inx = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gjx), from, to);
            res.X = inx[0];
            res.StandardDeviationOfX = inx[1];

            double[] iny = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gjy), from, to);
            res.Y = iny[0];
            res.StandardDeviationOfY = iny[1];

            double[] inz = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gjz), from, to);
            res.Z = inz[0];
            res.StandardDeviationOfZ = inz[1];

            return res;
        }
        private static double gjx(double[] para)
        {
            return g(p2tp(para[0], para[1]), tm).Magnitude 
                   * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0));
        }
        private static double gjy(double[] para)
        {
            return g(p2tp(para[0], para[1]), tm).Magnitude
                   * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0));
        }
        private static double gjz(double[] para)
        {
            return g(p2tp(para[0], para[1]), tm).Magnitude
                   * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0));
        }



        /// <summary>
        /// Dot Surface Integral for a vector defined by a function of position.
        /// </summary>        
        public static Scalar SurfaceIntegralDot(Vector v, Surface sf)
        {
            fp = v.VectorFunctionOfPosition;
            p2tp = sf.ParameterToPosition;
    
            Scalar res = new Scalar();
            
            double[] from = new double[2] { sf.Parameter1StartValue, sf.Parameter2StartValue };
            double[] to = new double[2] { sf.Parameter1EndValue, sf.Parameter2EndValue };

            double[] ins = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpdotsf), from, to);
            res.Magnitude = ins[0];
            res.StandardDeviation = ins[1];

            return res;
        }
        private static double fpdotsf(double[] para)
        {
            return fp(p2tp(para[0], para[1])).X
                   * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0)) +
                   fp(p2tp(para[0], para[1])).Y
                   * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0)) +
                   fp(p2tp(para[0], para[1])).Z
                   * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0));
        }


        /// <summary>
        /// Dot Surface Integral for a vector defined by a field.
        /// </summary> 
        public static Scalar SurfaceIntegralDot(Vector v, Surface sf, Time t)
        {
            f = v.VectorField;
            p2tp = sf.ParameterToPosition;
            tm.s = t.s;

            Scalar res = new Scalar();
            
            double[] from = new double[2] { sf.Parameter1StartValue, sf.Parameter2StartValue };
            double[] to = new double[2] { sf.Parameter1EndValue, sf.Parameter2EndValue };
       
            double[] ins = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fdotsf), from, to);
            res.Magnitude = ins[0];
            res.StandardDeviation = ins[1];
   
            return res;
        }
        private static double fdotsf(double[] para)
        {
            return f(p2tp(para[0], para[1]),tm).X
                   * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0)) +
                   f(p2tp(para[0], para[1]),tm).Y
                   * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0)) +
                   f(p2tp(para[0], para[1]),tm).Z
                   * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0));
        }


        /// <summary>
        /// Cross Surface Integral for a vector defined by a function of position.
        /// </summary>         
        public static Vector SurfaceIntegralCross(Vector v, Surface sf)
        {
            fp = v.VectorFunctionOfPosition;
            p2tp = sf.ParameterToPosition;
  
            Vector res = new Vector();
            
            double[] from = new double[2] { sf.Parameter1StartValue, sf.Parameter2StartValue };
            double[] to = new double[2] { sf.Parameter1EndValue, sf.Parameter2EndValue };

            double[] inx = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpcsx), from, to);
            res.X = inx[0];
            res.StandardDeviationOfX = inx[1];

            double[] iny = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpcsy), from, to);
            res.Y = iny[0];
            res.StandardDeviationOfY = iny[1];

            double[] inz = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpcsz), from, to);
            res.Z = inz[0];
            res.StandardDeviationOfZ = inz[1];

            return res;
        }
        private static double fpcsx(double[] para)
        {
            return fp(p2tp(para[0], para[1])).Y
                * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0))
                - fp(p2tp(para[0], para[1])).Z
                * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0));
        }
        private static double fpcsy(double[] para)
        {
            return fp(p2tp(para[0], para[1])).Z
                * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0))
                - fp(p2tp(para[0], para[1])).X
                * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0));
        }
        private static double fpcsz(double[] para)
        {
            return fp(p2tp(para[0], para[1])).X
                * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0))
                - fp(p2tp(para[0], para[1])).Y
                * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0));
        }



        /// <summary>
        /// Cross Surface Integral for a vector defined by a field.
        /// </summary>
        public static Vector SurfaceIntegralCross(Vector v, Surface sf, Time t)
        {
            f = v.VectorField;
            p2tp = sf.ParameterToPosition;
            tm.s = t.s;

            Vector res = new Vector();
           
            double[] from = new double[2] { sf.Parameter1StartValue, sf.Parameter2StartValue };
            double[] to = new double[2] { sf.Parameter1EndValue, sf.Parameter2EndValue };
         
            double[] inx = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fcsx), from, to);
            res.X = inx[0];
            res.StandardDeviationOfX = inx[1];

            double[] iny = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fcsy), from, to);
            res.Y = iny[0];
            res.StandardDeviationOfY = iny[1];

            double[] inz = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fcsz), from, to);
            res.Z = inz[0];
            res.StandardDeviationOfZ = inz[1];

            return res;
        }
        private static double fcsx(double[] para)
        {
            return f(p2tp(para[0], para[1]), tm).Y
                * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0))
                - f(p2tp(para[0], para[1]), tm).Z
                * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0));
        }
        private static double fcsy(double[] para)
        {
            return f(p2tp(para[0], para[1]), tm).Z
                * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0))
                - f(p2tp(para[0], para[1]), tm).X
                * (rd(para, 0, 0) * rd(para, 1, 1) - rd(para, 0, 1) * rd(para, 1, 0));
        }
        private static double fcsz(double[] para)
        {
            return f(p2tp(para[0], para[1]), tm).X
                * (rd(para, 2, 0) * rd(para, 0, 1) - rd(para, 2, 1) * rd(para, 0, 0))
                - f(p2tp(para[0], para[1]), tm).Y
                * (rd(para, 1, 0) * rd(para, 2, 1) - rd(para, 1, 1) * rd(para, 2, 0));
        }




        /// <summary>
        /// Volume Integral for a scalar defined by a function of position.
        /// </summary>
        public static Scalar VolumeIntegral(Scalar s, Volume vol)
        {
            gp = s.ScalarFunctionOfPosition;
            p3tp = vol.ParameterToPosition;
    
            Scalar res = new Scalar();
          
            double[] from = new double[3] { vol.Parameter1StartValue, vol.Parameter2StartValue, vol.Parameter3StartValue };
            double[] to = new double[3] { vol.Parameter1EndValue, vol.Parameter2EndValue, vol.Parameter3EndValue };
        
            double[] ins = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gpvol), from, to);
            res.Magnitude = ins[0];
            res.StandardDeviation = ins[1];

            return res;
        }
        private static double hx(double[] para)
        {
            Position x = p3tp(para[0], para[1], para[2]);
            return x.X;
        }
        private static double hy(double[] para)
        {
            Position x = p3tp(para[0], para[1], para[2]);
            return x.Y;
        }
        private static double hz(double[] para)
        {
            Position x = p3tp(para[0], para[1], para[2]);
            return x.Z;
        }
        private static double rd3(double[] para, int xyz, int uvw)   // xyz=0,1,2    uvw = 0,1,2
        {
            Calculus.Function<double[], double> xyzofp3 = new Calculus.Function<double[], double>(hx);
            if (xyz == 1) xyzofp3 = new Calculus.Function<double[], double>(hy);
            if (xyz == 2) xyzofp3 = new Calculus.Function<double[], double>(hz);
            return Calculus.PartialDifferentiation(xyzofp3, para, uvw);
        }
        private static double ja3(double[] para)
        {
            return (rd3(para, 0, 0) * rd3(para, 1, 1) * rd3(para, 2, 2)
                  + rd3(para, 0, 1) * rd3(para, 1, 2) * rd3(para, 2, 0)
                  + rd3(para, 1, 0) * rd3(para, 2, 1) * rd3(para, 0, 2)
                  - rd3(para, 0, 2) * rd3(para, 1, 1) * rd3(para, 2, 0)
                  - rd3(para, 1, 2) * rd3(para, 2, 1) * rd3(para, 0, 0)
                  - rd3(para, 0, 1) * rd3(para, 1, 0) * rd3(para, 2, 2));
        }
        private static double gpvol(double[] para)
        {
            return gp(p3tp(para[0], para[1], para[2])).Magnitude * Math.Abs(ja3(para));
        }

        /// <summary>
        /// Volume Integral for a scalar defined by a field.
        /// </summary>
        public static Scalar VolumeIntegral(Scalar s, Volume vol, Time t)
        {
            g = s.ScalarField;
            p3tp = vol.ParameterToPosition;
            tm.s = t.s;

            Scalar res = new Scalar();
            
            double[] from = new double[3] { vol.Parameter1StartValue, vol.Parameter2StartValue, vol.Parameter3StartValue};
            double[] to = new double[3] { vol.Parameter1EndValue, vol.Parameter2EndValue, vol.Parameter3EndValue};
        
            double[] ins = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(gvol), from, to);
            res.Magnitude = ins[0];
            res.StandardDeviation = ins[1];

            return res;
        }
        private static double gvol(double[] para)
        {
            return g(p3tp(para[0], para[1], para[2]),tm).Magnitude * Math.Abs(ja3(para));
        }



        /// <summary>
        /// Volume Integral for a vector defined by a function of position.
        /// </summary>
        public static Vector VolumeIntegral(Vector v, Volume vol)
        {
            fp = v.VectorFunctionOfPosition;
            p3tp = vol.ParameterToPosition;
     
            Vector res = new Vector();
          
            double[] from = new double[3] { vol.Parameter1StartValue, vol.Parameter2StartValue, vol.Parameter3StartValue };
            double[] to = new double[3] { vol.Parameter1EndValue, vol.Parameter2EndValue, vol.Parameter3EndValue };

            double[] inx = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpxvol), from, to);
            res.X = inx[0];
            res.StandardDeviationOfX = inx[1];

            double[] iny = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpyvol), from, to);
            res.Y = iny[0];
            res.StandardDeviationOfY = iny[1];

            double[] inz = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fpzvol), from, to);
            res.Z = inz[0];
            res.StandardDeviationOfZ = inz[1];

            return res;
        }
        private static double fpxvol(double[] para)
        {
            return fp(p3tp(para[0], para[1], para[2])).X * Math.Abs(ja3(para));
        }
        private static double fpyvol(double[] para)
        {
            return fp(p3tp(para[0], para[1], para[2])).Y * Math.Abs(ja3(para));
        }
        private static double fpzvol(double[] para)
        {
            return fp(p3tp(para[0], para[1], para[2])).Z * Math.Abs(ja3(para));
        }




        /// <summary>
        /// Volume Integral for a vector defined by a field.
        /// </summary>
        public static Vector VolumeIntegral(Vector v, Volume vol, Time t)
        {
            f = v.VectorField;
            p3tp = vol.ParameterToPosition;
            tm.s = t.s;

            Vector res = new Vector();
       
            double[] from = new double[3] { vol.Parameter1StartValue, vol.Parameter2StartValue, vol.Parameter3StartValue };
            double[] to = new double[3] { vol.Parameter1EndValue, vol.Parameter2EndValue, vol.Parameter3EndValue };
      
            double[] inx = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fxvol), from, to);
            res.X = inx[0];
            res.StandardDeviationOfX = inx[1];

            double[] iny = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fyvol), from, to);
            res.Y = iny[0];
            res.StandardDeviationOfY = iny[1];

            double[] inz = Calculus.IntegrationMultiD(new Calculus.Function<double[], double>(fzvol), from, to);
            res.Z = inz[0];
            res.StandardDeviationOfZ = inz[1];

            return res;
        }
        private static double fxvol(double[] para)
        {
            return f(p3tp(para[0], para[1], para[2]), tm).X * Math.Abs(ja3(para));
        }
        private static double fyvol(double[] para)
        {
            return f(p3tp(para[0], para[1], para[2]), tm).Y * Math.Abs(ja3(para));
        }
        private static double fzvol(double[] para)
        {
            return f(p3tp(para[0], para[1], para[2]), tm).Z * Math.Abs(ja3(para));
        } 
	}
}