Motion equation of a space engine in low-earth orbit

Question

My goal is to simulate the movement/motion of a space capsule from low-Earth orbit with initial conditions emulating an initial thrust. Here is the code, Is it right? I know that the next step is to implement a Runge Kutta algorithm (to solve the ordinary differential equation), but shouldn't I go to the state-space representation (knowing that the end goal is to implement a Kalman filter)?

r = sqrt (X.^2 + Y.^2 + Z^.2); % calculates the distance from the spacecraft to the origin                                                                                       (which could be the center of the Earth or some other reference point).
rm = sqrt (Xm.^2 + Ym.^2 + Zm^.2); %calculates the distance from the spacecraft to the Moon.
delta_m = sqrt ((X-Xm).^2 + (Y-Ym).^2 + (Z-Zm)^.2); %calculates the distance between the spacecraft and the Moon.
delta_s = sqrt ((X-Xs).^2 + (Y-Ys).^2 + (Z-Zs)^.2); %calculates the distance between the spacecraft and the Sun.
mu_e = 3.986135e14; %sets the gravitational parameter for the Earth.
mu_m = 4.89820e12;%sets the gravitational parameter for the Moon.
mu_s = 1.3253e20;%sets the gravitational parameter for the Sun.
a = 6.37826e6; %the equatorial radius of the Earth.
J = 1.6246e-3;%sets the Earth's second dynamic form factor.

X_2dot = -((mu_e*X)/r^3 ) * [1 + (J(a/r)^2)*(1-(5*Z^2/r^2))] - (mu_m(X-Xm)/delta_m^3) - (mu_m*Xm/rm^3) - (mu_s*(X-X_s)/delta_s^3) - (mu_s*Xs/r^3);
Y_2dot = -((mu_e*Y)/r^3 ) * [1 + (J(a/r)^2)*(1-(5*Z^2/r^2))] - (mu_m(Y-Ym)/delta_m^3) - (mu_m*Ym/rm^3) - (mu_s*(Y-Y_s)/delta_s^3) - (mu_s*Ys/r^3);
Z_2dot = -((mu_e*Z)/r^3 ) * [1 + (J(a/r)^2)*(1-(5*Z^2/r^2))] - (mu_m(Z-Zm)/delta_m^3) - (mu_m*Zm/rm^3) - (mu_s*(Z-Z_s)/delta_s^3) - (mu_s*Zs/r^3);

Answer 1

Here the code, Is it right?

"Please check my code" questions are not really appropriate in any Stack Exchange site.

Instead, you might first show the equations you're trying to implement along with citing your source for them, and then at least ask if you've coded them properly.

I know that the next step is to implement Runge Kutta's algorithm (to solve the ODE), but shouldn't I go to the state space representation (knowing that the ned goal is to implement a Kalman filter)?

Yes! Your ODE solver will solve an IVP or initial value problem. In classical mechanics (including orbital mechanics) that means you'll need both an initial position and an initial velocity.

Your initial state vector will then be $\mathbf{X} = [x, y, z, v_x, v_y, v_z]$

Your RK ODE solver will require a function that takes your state vector $\mathbf{X}$ and returns it's first time derivative $\text{d} \mathbf{X}/\text{d}t$ to integrate.

For the first three values, you just return the three velocities you were passed. For the second three, you return the accelerations that you are calculating in your question. $[v_x, v_y, v_y, a_x, a_y, a_z]$

If you are using Python with the state vector as a numpy array, that might look like this (from here) for a central force where GMe is the standard gravitational parameter of Earth:

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -GMe * x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

If you want to include J2, mine then looks like this (from here):

def deriv(X, t):
    x, v = X.reshape(2, -1)
    accs = []
    for xx in x.reshape(-1, 3):
        acc = -xx * GMe * ((xx**2).sum())**-1.5
        if use_J2:
            x, y, z = xx
            xsq, ysq, zsq = xx**2
            rm7 = (xsq + ysq + zsq)**-3.5
            accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))
            acc -= J2 * np.hstack((accJ2x, accJ2y, accJ2z))
        accs.append(acc)
    return np.hstack((v, np.hstack((accs))))