using System;

namespace Science.Physics.GeneralPhysics
{
	/// <summary>
	/// Tensor
	/// </summary>
	public class Tensor : Dimension
	{
        public delegate Tensor Field(Position x, Time t);
        public delegate Tensor FunctionOfPosition(Position x);
        public delegate Tensor FunctionOfTime(Time t);

		private double xx=0.0,xy=0.0,xz=0.0;
		private double yx=0.0,yy=0.0,yz=0.0;
		private double zx=0.0,zy=0.0,zz=0.0;
        private double dxx = 0.0, dxy = 0.0, dxz = 0.0;
        private double dyx = 0.0, dyy = 0.0, dyz = 0.0;
        private double dzx = 0.0, dzy = 0.0, dzz = 0.0;
		private bool xxv=false, xyv=false, xzv=false;
		private bool yxv=false, yyv=false, yzv=false;
		private bool zxv=false, zyv=false, zzv=false;
		
		public Tensor()
		{	
		}	

		public Tensor(Vector a, Vector b)  // OuterProduct
		{
			xx = a.X * b.X;
			xy = a.X * b.Y;
			xz = a.X * b.Z;
			yx = a.Y * b.X;
			yy = a.Y * b.Y;
			yz = a.Y * b.Z;
			zx = a.Z * b.X;
			zy = a.Z * b.Y;
			zz = a.Z * b.Z;
		}

        private Tensor.Field g;
        public Tensor(Tensor.Field f)
		{
			g = f;
		}
        public Tensor.Field TensorField
        {
            get { return g; }
            set { g = value; }
        }

        private Tensor.FunctionOfPosition gp;
        public Tensor(Tensor.FunctionOfPosition fp)
		{
			gp = fp;
		}
        public Tensor.FunctionOfPosition TensorFunctionOfPosition
        {
            get { return gp; }
            set { gp = value; }
        }

        private Tensor.FunctionOfTime gt;
        public Tensor(Tensor.FunctionOfTime ft)
        {
            gt = ft;
        }
        public Tensor.FunctionOfTime TensorFunctionOfTime
        {
            get { return gt; }
            set { gt = value; }
        }

		public void Set(Position p, Time tt)
		{
            Tensor t = g(p, tt);
            xx = t.XX;
            xy = t.XY;
            xz = t.XZ;
            yx = t.YX;
            yy = t.YY;
            yz = t.YZ;
            zx = t.ZX;
            zy = t.ZY;
            zz = t.ZZ;
		}

		public void Set(Position p)
        {
            Tensor t = gp(p);
            xx = t.XX;
            xy = t.XY;
            xz = t.XZ;
            yx = t.YX;
            yy = t.YY;
            yz = t.YZ;
            zx = t.ZX;
            zy = t.ZY;
            zz = t.ZZ;
		}

        public void Set(Time t)
        {
            Tensor tt = gt(t);
            xx = tt.XX;
            xy = tt.XY;
            xz = tt.XZ;
            yx = tt.YX;
            yy = tt.YY;
            yz = tt.YZ;
            zx = tt.ZX;
            zy = tt.ZY;
            zz = tt.ZZ;
        }


		public double XX
		{
			get{return xx;}
			set{xx=value;}
		}
		public double XY
		{
			get{return xy;}
			set{xy=value;}
		}
		public double XZ
		{
			get{return xz;}
			set{xz=value;}
		}
		public double YX
		{
			get{return yx;}
			set{yx=value;}
		}
		public double YY
		{
			get{return yy;}
			set{yy=value;}
		}
		public double YZ
		{
			get{return yz;}
			set{yz=value;}
		}
		public double ZX
		{
			get{return zx;}
			set{zx=value;}
		}
		public double ZY
		{
			get{return zy;}
			set{zy=value;}
		}
		public double ZZ
		{
			get{return zz;}
			set{zz=value;}
		}

        public double StandardDeviationXX
        {
            get { return dxx; }
            set { dxx = value; }
        }
        public double StandardDeviationXY
        {
            get { return dxy; }
            set { dxy = value; }
        }
        public double StandardDeviationXZ
        {
            get { return dxz; }
            set { dxz = value; }
        }
        public double StandardDeviationYX
        {
            get { return dyx; }
            set { dyx = value; }
        }
        public double StandardDeviationYY
        {
            get { return dyy; }
            set { dyy = value; }
        }
        public double StandardDeviationYZ
        {
            get { return dyz; }
            set { dyz = value; }
        }
        public double StandardDeviationZX
        {
            get { return dzx; }
            set { dzx = value; }
        }
        public double StandardDeviationZY
        {
            get { return dzy; }
            set { dzy = value; }
        }
        public double StandardDeviationZZ
        {
            get { return dzz; }
            set { dzz = value; }
        }

		public bool XXVariableQ
		{
			get{return xxv;}
			set{xxv=value;}
		}
		public bool XYVariableQ
		{
			get{return xyv;}
			set{xyv=value;}
		}
		public bool XZVariableQ
		{
			get{return xzv;}
			set{xzv=value;}
		}		
		public bool YXVariableQ
		{
			get{return yxv;}
			set{yxv=value;}
		}
		public bool YYVariableQ
		{
			get{return yyv;}
			set{yyv=value;}
		}
		public bool YZVariableQ
		{
			get{return yzv;}
			set{yzv=value;}
		}		
		public bool ZXVariableQ
		{
			get{return zxv;}
			set{zxv=value;}
		}
		public bool ZYVariableQ
		{
			get{return zyv;}
			set{zyv=value;}
		}
		public bool ZZVariableQ
		{
			get{return zzv;}
			set{zzv=value;}
		}		


		public static Tensor operator +(Tensor a, Tensor b)
		{
			Tensor c = new Tensor();
			c.XX = a.XX + b.XX;
			c.XY = a.XY + b.XY;
			c.XZ = a.XZ + b.XZ;
			c.YX = a.YX + b.YX;
			c.YY = a.YY + b.YY;
			c.YZ = a.YZ + b.YZ;
			c.ZX = a.ZX + b.ZX;
			c.ZY = a.ZY + b.ZY;
			c.ZZ = a.ZZ + b.ZZ;
			return c;
		}
		public static Tensor operator -(Tensor a, Tensor b)
		{
			Tensor c = new Tensor();
			c.XX = a.XX - b.XX;
			c.XY = a.XY - b.XY;
			c.XZ = a.XZ - b.XZ;
			c.YX = a.YX - b.YX;
			c.YY = a.YY - b.YY;
			c.YZ = a.YZ - b.YZ;
			c.ZX = a.ZX - b.ZX;
			c.ZY = a.ZY - b.ZY;
			c.ZZ = a.ZZ - b.ZZ;
			return c;
		}
		public static Tensor operator *(Tensor a, Tensor b)
		{
			Tensor c = new Tensor();
			c.XX = a.XX*b.XX + a.XY*b.YX + a.XZ*b.ZX;
			c.XY = a.XX*b.XY + a.XY*b.YY + a.XZ*b.ZY;
			c.XZ = a.XX*b.XZ + a.XY*b.YZ + a.XZ*b.ZZ;
			c.YX = a.YX*b.XX + a.YY*b.YX + a.YZ*b.ZX;
			c.YY = a.YX*b.XY + a.YY*b.YY + a.YZ*b.ZY;
			c.YZ = a.YX*b.XZ + a.YY*b.YZ + a.YZ*b.ZZ;
			c.ZX = a.ZX*b.XX + a.ZY*b.YX + a.ZZ*b.ZX;
			c.ZY = a.ZX*b.XY + a.ZY*b.YY + a.ZZ*b.ZY;
			c.ZZ = a.ZX*b.XZ + a.ZY*b.YZ + a.ZZ*b.ZZ;
			return c;
		}

		public static Vector operator *(Tensor a, Vector b)
		{
			Vector c = new Vector();
			c.X = a.XX*b.X + a.XY*b.Y + a.XZ*b.Z;
			c.Y = a.YX*b.X + a.YY*b.Y + a.YZ*b.Z;
			c.Z = a.ZX*b.X + a.ZY*b.Y + a.ZZ*b.Z;
			return c;
		}

		public static Vector operator *(Vector a, Tensor b)
		{
			Vector c = new Vector();
			c.X = a.X*b.XX + a.Y*b.YX + a.Z*b.ZX;
			c.Y = a.X*b.XY + a.Y*b.YY + a.Z*b.ZY;
			c.Z = a.X*b.XZ + a.Y*b.YZ + a.Z*b.ZZ;
			return c;
		}
		
		public static Tensor operator *(Scalar a, Tensor b)
		{
			Tensor c = new Tensor();
			c.XX = a.Magnitude * b.XX;
			c.XY = a.Magnitude * b.XY;
			c.XZ = a.Magnitude * b.XZ;
			c.YX = a.Magnitude * b.YX;
			c.YY = a.Magnitude * b.YY;
			c.YZ = a.Magnitude * b.YZ;
			c.ZX = a.Magnitude * b.ZX;
			c.ZY = a.Magnitude * b.ZY;
			c.ZZ = a.Magnitude * b.ZZ;
			return c;
		}		

		public static Tensor operator *(Tensor b, Scalar a)
		{
			Tensor c = new Tensor();
			c.XX = a.Magnitude * b.XX;
			c.XY = a.Magnitude * b.XY;
			c.XZ = a.Magnitude * b.XZ;
			c.YX = a.Magnitude * b.YX;
			c.YY = a.Magnitude * b.YY;
			c.YZ = a.Magnitude * b.YZ;
			c.ZX = a.Magnitude * b.ZX;
			c.ZY = a.Magnitude * b.ZY;
			c.ZZ = a.Magnitude * b.ZZ;
			return c;
		}		
		public static Tensor operator *(double a, Tensor b)
		{
			Tensor c = new Tensor();
			c.XX = a * b.XX;
			c.XY = a * b.XY;
			c.XZ = a * b.XZ;
			c.YX = a * b.YX;
			c.YY = a * b.YY;
			c.YZ = a * b.YZ;
			c.ZX = a * b.ZX;
			c.ZY = a * b.ZY;
			c.ZZ = a * b.ZZ;
			return c;
		}		

		public static Tensor operator *(Tensor b, double a)
		{
			Tensor c = new Tensor();
			c.XX = a * b.XX;
			c.XY = a * b.XY;
			c.XZ = a * b.XZ;
			c.YX = a * b.YX;
			c.YY = a * b.YY;
			c.YZ = a * b.YZ;
			c.ZX = a * b.ZX;
			c.ZY = a * b.ZY;
			c.ZZ = a * b.ZZ;
			return c;
		}		
		public Tensor Inverse
		{
			get
			{				
				Tensor im = new Tensor();
				double[,] a = new double[3,3]
                        {{this.XX,this.XY,this.XZ},
						 {this.YX,this.YY,this.YZ},
						 {this.ZX,this.ZY,this.ZZ}};
                LUdcmp lu = new LUdcmp(a);
				double[,] element = lu.inverse();
                im.XX = element[0, 0];
                im.XY = element[0, 1];
                im.XZ = element[0, 2];
                im.YX = element[1, 0];
                im.YY = element[1, 1];
                im.YZ = element[1, 2];
                im.ZX = element[2, 0];
                im.ZY = element[2, 1];
                im.ZZ = element[2, 2];
          		return im;
			}
		}

        public override string ToString()
        {
            string res = "";
            res += this.XX.ToString() + " +/- " + this.StandardDeviationXX.ToString() + "     ";
            if (this.XY >= 0.0) res += " +";
            res += this.XY.ToString() + " +/- " + this.StandardDeviationXY.ToString() + "     ";
            if (this.XZ >= 0.0) res += " +";
            res += this.XZ.ToString() + " +/- " + this.StandardDeviationXZ.ToString() + " \r\n";
            res += this.YX.ToString() + " +/- " + this.StandardDeviationYX.ToString() + "     ";
            if (this.YY >= 0.0) res += " +";
            res += this.YY.ToString() + " +/- " + this.StandardDeviationYY.ToString() + "     ";
            if (this.YZ >= 0.0) res += " +";
            res += this.YZ.ToString() + " +/- " + this.StandardDeviationYZ.ToString() + " \r\n";
            res += this.ZX.ToString() + " +/- " + this.StandardDeviationZX.ToString() + "     ";
            if (this.ZY >= 0.0) res += " +";
            res += this.ZY.ToString() + " +/- " + this.StandardDeviationZY.ToString() + "     ";
            if (this.ZZ >= 0.0) res += " +";
            res += this.ZZ.ToString() + " +/- " + this.StandardDeviationZZ.ToString() + " \r\n";
            return res;
        }
	}
}