|
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;
|