using System;

namespace Science.Physics.GeneralPhysics
{
	/// <summary>
	/// MomentOfInertia
	/// </summary>
	public class MomentOfInertia : Tensor
	{
		private void SetDim()
		{
			this.DimensionMass = 1;
			this.DimensionLength = 2;
			this.DimensionTime = 0;
		}
		public double kgmSQUARE
		{
			get{return this.ZZ;}
		}
		public MomentOfInertia()
		{
			SetDim();
		}	
		public MomentOfInertia(Mass[] m, Position[] r)
		{
			SetDim();
			double xx = 0.0, xy = 0.0, xz = 0.0;
			double yy = 0.0, yz = 0.0, zz = 0.0;
		    for(int i = 0; i < m.Length; i++)
			{
				xx += m[i].kg*(r[i].Y*r[i].Y+r[i].Z*r[i].Z);
				xy += - m[i].kg*r[i].X*r[i].Y;
				xz += - m[i].kg*r[i].X*r[i].Z;
				yy += m[i].kg*(r[i].X*r[i].X+r[i].Z*r[i].Z);
				yz += - m[i].kg*r[i].Y*r[i].Z;
				zz += m[i].kg*(r[i].X*r[i].X+r[i].Y*r[i].Y);
			}
			this.XX = xx; 
			this.XY = xy;
			this.XZ = xz;
			this.YX = xy; 
			this.YY = yy;
			this.YZ = yz;
			this.ZX = xz; 
			this.ZY = yz;
			this.ZZ = zz;
		}

		public static void ParallelAxisTheorem(MomentOfInertia I, 
			Vector axisUnitVector, Mass m, Length d)  // after set axis
		{
			if(axisUnitVector.X > 0.9) I.XX = I.XX + m.kg*d.m*d.m;
			else if(axisUnitVector.Y > 0.9) I.YY = I.YY + m.kg*d.m*d.m;
			else I.ZZ = I.ZZ + m.kg*d.m*d.m;
		}

        private Scalar.FunctionOfPosition fp;
		public MomentOfInertia(Density rho, Volume v)
		{
			SetDim();
		    fp = rho.ScalarFunctionOfPosition;

            Scalar s = new Scalar();

            s.ScalarFunctionOfPosition = new Scalar.FunctionOfPosition(XXrho);
            Scalar ss1 = VectorCalculus.VolumeIntegral(s, v);			
			this.XX = ss1.Magnitude;
            this.StandardDeviationXX = ss1.StandardDeviation;

            s.ScalarFunctionOfPosition = new Scalar.FunctionOfPosition(XYrho);
            Scalar ss2 = VectorCalculus.VolumeIntegral(s, v);
            this.XY = ss2.Magnitude;
            this.StandardDeviationXY = ss2.StandardDeviation;

            s.ScalarFunctionOfPosition = new Scalar.FunctionOfPosition(XZrho);
            Scalar ss3 = VectorCalculus.VolumeIntegral(s, v);
            this.XZ = ss3.Magnitude;
            this.StandardDeviationXZ = ss3.StandardDeviation;

            s.ScalarFunctionOfPosition = new Scalar.FunctionOfPosition(YYrho);
            Scalar ss4 = VectorCalculus.VolumeIntegral(s, v);
            this.YY = ss4.Magnitude;
            this.StandardDeviationYY = ss4.StandardDeviation;

            s.ScalarFunctionOfPosition = new Scalar.FunctionOfPosition(YZrho);
            Scalar ss5 = VectorCalculus.VolumeIntegral(s, v);
            this.YZ = ss5.Magnitude;
            this.StandardDeviationYZ = ss5.StandardDeviation;

            s.ScalarFunctionOfPosition = new Scalar.FunctionOfPosition(ZZrho);
            Scalar ss6 = VectorCalculus.VolumeIntegral(s, v);
            this.ZZ = ss6.Magnitude;
            this.StandardDeviationZZ = ss6.StandardDeviation;
            
			this.YX = this.XY;
			this.ZX = this.XZ;
			this.ZY = this.YZ;

            this.StandardDeviationYX = this.StandardDeviationXY;
            this.StandardDeviationZX = this.StandardDeviationXZ;
            this.StandardDeviationZY = this.StandardDeviationYZ;
		}
	
		private Scalar XXrho(Position x)
		{
            Scalar s = new Scalar();
            s.Magnitude = (x.Y*x.Y + x.Z*x.Z)*fp(x).Magnitude;
            return s;
		}
        private Scalar XYrho(Position x)
		{
            Scalar s = new Scalar();
            s.Magnitude = -x.X * x.Y * fp(x).Magnitude;
            return s;
		}
        private Scalar XZrho(Position x)
		{
            Scalar s = new Scalar();
            s.Magnitude = -x.X * x.Z * fp(x).Magnitude;
            return s;
		}
        private Scalar YYrho(Position x)
		{
            Scalar s = new Scalar();
            s.Magnitude = (x.X * x.X + x.Z * x.Z) * fp(x).Magnitude;
            return s;
		}
        private Scalar YZrho(Position x)
		{
            Scalar s = new Scalar();
            s.Magnitude = -x.Y * x.Z * fp(x).Magnitude;
            return s;
		}
        private Scalar ZZrho(Position x)
		{
            Scalar s = new Scalar();
            s.Magnitude = (x.X * x.X + x.Y * x.Y) * fp(x).Magnitude;
            return s;
		}

        public override string ToString()
        {
            return base.ToString() + "(kgm^2)";
        }
	}
}