small additions

Author: elv1s42

RK/GravitationalMomentEquations.cs

@@ -11,18 +11,18 @@ public GravitationalMomentEquations(double epsilon, double delta) : base(7)
             _delta = delta;
         }
 
-        public override void Equations(double t, double[] y, ref double[] fy)
+        public override void Equations(double t, double[] y, ref double[] yDot)
         {
             var p = y[0] * y[0] + y[1] * y[1] + y[2] * y[2] + y[3] * y[3] - 1;
-            fy[0] = -0.5 * y[1] * y[4] - 0.5 * y[2] * y[5] - 0.5 * y[3] * y[6] + 0.5 * y[2] - 0.5 * y[0] * p;
-            fy[1] = 0.5 * y[0] * y[4] + 0.5 * y[2] * y[6] - 0.5 * y[3] * y[5] - 0.5 * y[3] - 0.5 * y[1] * p;
-            fy[2] = 0.5 * y[0] * y[5] + 0.5 * y[3] * y[4] - 0.5 * y[1] * y[6] - 0.5 * y[0] - 0.5 * y[2] * p;
-            fy[3] = 0.5 * y[0] * y[6] + 0.5 * y[1] * y[5] - 0.5 * y[2] * y[4] + 0.5 * y[1] - 0.5 * y[3] * p;
-            fy[4] = (_epsilon - _delta) * (-y[5] * y[6] + 3 * (2 * y[2] * y[3] + 2 * y[0] * y[1]) 
+            yDot[0] = -0.5 * y[1] * y[4] - 0.5 * y[2] * y[5] - 0.5 * y[3] * y[6] + 0.5 * y[2] - 0.5 * y[0] * p;
+            yDot[1] = 0.5 * y[0] * y[4] + 0.5 * y[2] * y[6] - 0.5 * y[3] * y[5] - 0.5 * y[3] - 0.5 * y[1] * p;
+            yDot[2] = 0.5 * y[0] * y[5] + 0.5 * y[3] * y[4] - 0.5 * y[1] * y[6] - 0.5 * y[0] - 0.5 * y[2] * p;
+            yDot[3] = 0.5 * y[0] * y[6] + 0.5 * y[1] * y[5] - 0.5 * y[2] * y[4] + 0.5 * y[1] - 0.5 * y[3] * p;
+            yDot[4] = (_epsilon - _delta) * (-y[5] * y[6] + 3 * (2 * y[2] * y[3] + 2 * y[0] * y[1]) 
                 * (y[0] * y[0] - y[1] * y[1] - y[2] * y[2] + y[3] * y[3]));
-            fy[5] = (1 - _epsilon) / _delta * (-y[6] * y[4] + 3 * (2 * y[1] * y[3] - 2 * y[0] * y[2]) 
+            yDot[5] = (1 - _epsilon) / _delta * (-y[6] * y[4] + 3 * (2 * y[1] * y[3] - 2 * y[0] * y[2]) 
                 * (y[0] * y[0] - y[1] * y[1] - y[2] * y[2] + y[3] * y[3]));
-            fy[6] = (_delta - 1) / _epsilon * (-y[4] * y[5] + 3 * (2 * y[1] * y[3] - 2 * y[0] * y[2]) 
+            yDot[6] = (_delta - 1) / _epsilon * (-y[4] * y[5] + 3 * (2 * y[1] * y[3] - 2 * y[0] * y[2]) 
                 * (2 * y[2] * y[3] + 2 * y[0] * y[1]));
         }
     }

RK/RungeKuttaBase.cs

@@ -5,13 +5,12 @@ namespace RK
     public class RungeKuttaBase
     {
         public int N;
-        double _t; // текущее время 
-        public double[] Y; // искомое решение Y[0] - само решение, Y[i] - i-тая производная решения
+        double _t; // current time 
+        public double[] Y; //solutions
 
-        readonly double[] _yy; 
-        double[] _y1, _y2, _y3, _y4;
+        private double[] _y1, _y2, _y3, _y4, _yy;
 
-        public RungeKuttaBase(int aN) // aN - размерность системы 
+        private void BaseSetup(int aN)
         {
             N = aN; // сохранить размерность системы
 
@@ -29,6 +28,17 @@ public RungeKuttaBase(int aN) // aN - размерность системы
             _y4 = new double[N];
         }
 
+        public RungeKuttaBase(int aN) // aN - размерность системы 
+        {
+            BaseSetup(aN);
+        }
+
+        public RungeKuttaBase(int aN, double t0, double[] y0)
+        {
+            BaseSetup(aN);
+            SetInit(t0, y0);
+        }
+
         public void SetInit(double t0, double[] y0) // установить начальные условия.
         {                                           // t0 - начальное время, Y0 - начальное условие
             _t = t0;
@@ -39,12 +49,12 @@ public void SetInit(double t0, double[] y0) // установить началь
             }
         }
 
-        public double GetTime() // вернуть текущее время
+        public double GetTime()
         {
             return _t;
         }
 
-        public virtual void Equations(double t, double[] y, ref double[] fy)
+        public virtual void Equations(double t, double[] y, ref double[] yDot)
         {
             throw new NotImplementedException();
         }

RK/Solver.cs

@@ -6,31 +6,57 @@ namespace RK
 {
     public static class Solver
     {
-        static double[] SolveOnePoint(double tInit, double[] yInit, 
+        public const double Pi = Math.PI;
+
+        public static double[] InitialConditions(double phi0, double theta0, double psi0)
+        {
+            var sPh0 = Math.Sin(phi0 * Pi / 2);
+            var sPs0 = Math.Sin(psi0 * Pi / 2);
+            var sTh0 = Math.Sin(theta0 * Pi / 2);
+            var cPh0 = Math.Cos(phi0 * Pi / 2);
+            var cPs0 = Math.Cos(psi0 * Pi / 2);
+            var cTh0 = Math.Cos(theta0 * Pi / 2);
+
+            var lambda0 = cPh0*cPs0*cTh0 + sPh0*sPs0*sTh0;
+            var lambda1 = sPh0*cPs0*cTh0 - cPh0*sPs0*sTh0;
+            var lambda2 = cPh0*cPs0*sTh0 + sPh0*sPs0*cTh0;
+            var lambda3 = cPh0*sPs0*cTh0 - sPh0*cPs0*sTh0;
+            return new[] { lambda0, lambda1, lambda2, lambda3, 0, 1, 0 };
+            
+        }
+
+        private static double Psi(this double[] y)
+        {
+            var alpha1 = y[0] * y[0] + y[1] * y[1] - y[2] * y[2] - y[3] * y[3];
+            var beta1 = 2 * (y[1] * y[2] + y[0] * y[3]);
+            return Math.Abs(Math.Atan(beta1 / alpha1));
+        }
+
+        private static double Phi(this double[] y)
+        {
+            var gamma2 = 2 * (y[2] * y[3] + y[0] * y[1]);
+            var gamma3 = y[0] * y[0] - y[1] * y[1] - y[2] * y[2] + y[3] * y[3];
+            return Math.Abs(Math.Atan(gamma2 / gamma3));
+        }
+
+        private static double Theta(this double[] y)
+        {
+            var gamma1 = 2 * (y[1] * y[3] - y[0] * y[2]);
+            return Math.Abs(-Math.Asin(gamma1));
+        }
+
+        static void SolveOnePoint(double tInit, double[] yInit, 
             double epsilon, double delta, 
             double timePeriod, double timeStep)
         {
             var rk4 = new GravitationalMomentEquations(epsilon, delta);
             rk4.SetInit(tInit, yInit);
-            var write = new StreamWriter("Test.txt");
-
+            
             while (rk4.GetTime() < timePeriod)
             {
-                write.WriteLine("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}", 
-                    Convert.ToString(rk4.GetTime(), CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[0], CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[1], CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[2], CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[3], CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[4], CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[5], CultureInfo.InvariantCulture), 
-                    Convert.ToString(rk4.Y[6], CultureInfo.InvariantCulture));
-
+                //rk4.Y
                 rk4.NextStep(timeStep);
             }
-            
-            write.Close();
-            return rk4.Y;
         }
     }
 }