Viaje de la Tierra a la Luna

Sistema de ecuaciones diferenciales de segundo orden

Resolver el sistema de ecuaciones diferenciales de segundo orden

$\begin{array}{l} \frac{d^{2} x_{R}}{d t^{2}} - 2 \cdot 0.231 \frac{d y_{R}}{d t} - {0.231}^{2} x_{R} = - 11519.568 \frac{(x_{R} + 0.731)}{d_{T}^{3}} - 141.394 \frac{(x_{R} - 59.552)}{d_{L}^{3}} \\ \frac{d^{2} y_{R}}{d t^{2}} + 2 \cdot 0.231 \frac{d x_{R}}{d t} - {0.231}^{2} y_{R} = - 11519.568 \frac{y_{R}}{d_{T}^{3}} - 141.394 \frac{y_{R}}{d_{L}^{3}} \\ d_{T} = \sqrt{{(x_{R} + 0.731)}^{2} + y_{R}^{2}} d_{L} = \sqrt{{(x_{R} - 59.552)}^{2} + y_{R}^{2}} \end{array}$

La constante J se escribe

$J = \frac{1}{2} ({\dot{x}}_{R}^{2} + {\dot{y}}_{R}^{2}) - \frac{1}{2} {0.231}^{2} (x_{R}^{2} + y_{R}^{2}) - (\frac{11519.568}{d_{T}} + \frac{141.394}{d_{L}})$

En el instante t=0, la nave espacial parte de la posición

x₀=-r_T+r·cosθ
y₀= r·sinθ

Con velocidad inicial

v_0x=-(v+Δv)·sinθ
v_0y=(v+Δv)·cosθ

Donde Δv es el incremento de velocidad proporcionado por los cohetes de la nave de forma casi instantánea

public class Sistema extends RungeKutta{
    final double  rT=0.73095;
    final double rL=59.55162;
    final double w=0.23072;  // wTierra-Luna*un día
    final double cL=141.39403;
    final double cT=11519.56835;
    pulic Sistema(double h){
      super(h);
    }
    public double f(double x, double y, double vx, double vy, double t){
        double dL=Math.sqrt((x-rL)*(x-rL)+y*y);
        double dT=Math.sqrt((x+rT)*(x+rT)+y*y);
        double z=2*w*vy+w*w*x-cL*(x-rL)/(dL*dL*dL)-cT*(x+rT)/(dT*dT*dT);
         return z;
    }
    public double g(double x, double y, double vx, double vy, double t){
        double dL=Math.sqrt((x-rL)*(x-rL)+y*y);
        double dT=Math.sqrt((x+rT)*(x+rT)+y*y);
        double z=-2*w*vx+w*w*y-cL*y/(dL*dL*dL)-cT*y/(dT*dT*dT);
        return z;
    }
    public double energia(double x, double y, double vx, double vy){
        double dL=Math.sqrt((x-rL)*(x-rL)+y*y);
        double dT=Math.sqrt((x+rT)*(x+rT)+y*y);
       double z=(vx*vx+vy*vy)/2-(cL/dL+cT/dT)-w*w*(x*x+y*y)/2;
       return z;
    }
}
//Objeto de la clase Sistema
 void setNuevo(double v0, double ang){
    double x0=-rT+rNave*Math.cos(ang);
    double y0=rNave*Math.sin(ang);
    double v0x=-(rNave*wNave+v0*3600*24/6.37e6)*Math.sin(ang);
    double v0y=(rNave*wNave+v0*3600*24/6.37e6)*Math.cos(ang);
    estado=new Estado(0.0, x0, y0, v0x, v0y);
    sis=new Sistema(calculaPaso());
    eInicial=oscilador.energia(x0, y0, v0x, v0y);
 }
  double calculaPaso(){
     double dL=Math.sqrt((estado.x-rL)*(estado.x-rL)+estado.y*estado.y);
     double dT=Math.sqrt((estado.x+rT)*(estado.x+rT)+estado.y*estado.y);
     double distMed=Math.sqrt(dT*dL);
     double dt=0.0005*Math.pow((distMed/(60.282*0.1)), 1.5);  //para un día
     return dt;
 }
 void mover(){
        sis.resolver(estado);
        sis.setPaso(calculaPaso());
        double energia=oscilador.energia(estado.x, estado.y, estado.vx, estado.vy);
        error=Math.abs((energia-eInicial)*100/eInicial);
        if(error>1.0){
         //se para ;
        }
   }
 }