ホーム 道しるべ 憩いの広場 濃緑空間 濃緑研の日記

RungeKutta法プログラム
ホーム ] 上へ ] オイラ−法プログラム ] オイラー法(3DNow!) ] [ RungeKutta法プログラム ] オイラー法FPUアセンブリ ]

 

            R = sat.P.X*sat.P.X+sat.P.Y*sat.P.Y+sat.P.Z*sat.P.Z;
            w = float(GMdt/(sqrt(R)*R));
            RK[1].aX = w*sat.P.X;
            RK[1].vX = sat.V.X*dt;
            RK[1].aY = w*sat.P.Y;
            RK[1].vY = sat.V.Y*dt;
            RK[1].aZ = w*sat.P.Z;
            RK[1].vZ = sat.V.Z*dt;
            RK[1].x = sat.P.X + RK[1].vX/2;
            RK[1].y = sat.P.Y + RK[1].vY/2;
            RK[1].z = sat.P.Z + RK[1].vZ/2;

            R = RK[1].x*RK[1].x + RK[1].y*RK[1].y + RK[1].z*RK[1].z;
            w = float(GMdt/(sqrt(R)*R));
            RK[2].aX = w*RK[1].x;
            RK[2].vX = (sat.V.X+RK[1].aX/2)*dt;
            RK[2].aY = w*RK[1].y;
            RK[2].vY = (sat.V.Y+RK[1].aY/2)*dt;
            RK[2].aZ = w*RK[1].z;
            RK[2].vZ = (sat.V.Z+RK[1].aZ/2)*dt;
            RK[2].x = sat.P.X + RK[2].vX/2;
            RK[2].y = sat.P.Y + RK[2].vY/2;
            RK[2].z = sat.P.Z + RK[2].vZ/2;
           
            R = RK[2].x*RK[2].x + RK[2].y*RK[2].y + RK[2].z*RK[2].z;
            w = float(GMdt/(sqrt(R)*R));
            RK[3].aX = w*RK[2].x;
            RK[3].vX = (sat.V.X+RK[2].aX/2)*dt;
            RK[3].aY = w*RK[2].y;
            RK[3].vY = (sat.V.Y+RK[2].aY/2)*dt;
            RK[3].aZ = w*RK[2].z;
            RK[3].vZ = (sat.V.Z+RK[2].aZ/2)*dt;
            RK[3].x = sat.P.X + RK[3].vX/2;
            RK[3].y = sat.P.Y + RK[3].vY/2;
            RK[3].z = sat.P.Z + RK[3].vZ/2;

            R = RK[3].x*RK[3].x + RK[3].y*RK[3].y + RK[3].z*RK[3].z;
            w = float(GMdt/(sqrt(R)*R));
            RK[4].aX = w*RK[3].x;
            RK[4].vX = (sat.V.X+RK[3].aX)*dt;
            RK[4].aY = w*RK[3].y;
            RK[4].vY = (sat.V.Y+RK[3].aY)*dt;
            RK[4].aZ = w*RK[3].z;
            RK[4].vZ = (sat.V.Z+RK[3].aZ)*dt;

            sat.V.X += (RK[1].aX+2*(RK[2].aX+RK[3].aX)+RK[4].aX)/6;
            sat.P.X += (RK[1].vX+2*(RK[2].vX+RK[3].vX)+RK[4].vX)/6;
            sat.V.Y += (RK[1].aY+2*(RK[2].aY+RK[3].aY)+RK[4].aY)/6;
            sat.P.Y += (RK[1].vY+2*(RK[2].vY+RK[3].vY)+RK[4].vY)/6;
            sat.V.Z += (RK[1].aZ+2*(RK[2].aZ+RK[3].aZ)+RK[4].aZ)/6;
            sat.P.Z += (RK[1].vZ+2*(RK[2].vZ+RK[3].vZ)+RK[4].vZ)/6;