Motion equation of a space engin in low earth orbit

So my goal is to simulate the mouvement/motion of a space capsule from low Earth orbit with initial coditions emuating initial thrust.

Here the code,

Is it right? I know that the next step is to implement Runge Kutt'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 Klaman filter), please help!

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

[James Tursa]

You need to show us the differential equations you are trying to implement, along with your perturbation equations for J2, Moon, and Sun.

[Bjorn Gustavsson]

The terms for the effect of the equatorial bulging doesn't look right - they shouldn't all look the same, right?

[Blob]

@James Tursa here are the differential equations

[Bjorn Gustavsson]

Once you implement these equations and solve them with different RK-methods you should also test that the constants of motion are preserved. That is not guaranteed with standard RK-methods that "knows nothing about mechanics".

[James Tursa]

@Bjorn Gustavsson Maybe do this with all perturbations turned off first so that orbit constants can be checked. Once that gets working and tuned, then add in the perturbations where the orbit constants will not be preserved.

[Bjorn Gustavsson]

Yes, especially the forces from the moon and the sun should be turned off, and first also the J2 term. Then "everything" should be conserved. When J2 is included the precession (and possibly the shift of the perigee) will appear - in a controllable manner.

[Blob]

@James Tursa @Bjorn GustavssonIn order to solve the differential equations I need initial conditions, where do I find them? Do I pu the coordinate to zero? (I am sorry for all the questions, but it i a bit challenging for me!)

[James Tursa]

@Blob No you definitely don't put initial conditions to zero. Your differential equations are singular there, and in fact are not even valid anywhere inside the surface of the Earth. You need to define your initial conditions based on the problem you are solving. What is this "space capsule"? What type of orbit is it in? That will give you the initial conditions.

[Bjorn Gustavsson]

@Blob: In LEO you know that the orbit is close to circular and given an initial circular orbit at an altitude of, let's say, 450 km of altitude you can work out an orbit velocity. If you put the initial position above the equator, you can determine what inclination your initial orbit-plane will have by selecting the direction of the velocity - keeping in mind that it should be tangential to the surface of Earth.

[Blob]

That is my goal: simulating the ideal motion of the space capsule from low Earth orbit with initial conditions emulating the initial thrust, then implement the global nonlinear system accounting for noisy dynamics and noisy observations as well.

Answer 1

[James Tursa]

As Bjorn has pointed out, your J2 terms are all the same but they should be different because the equatorial bulge affects the z-axis differently than it affects the x-y axes.. E.g., your posted equations for Z_2dot have the equivalent of 3-(5*Z^2/r^2) but your code has 1-(5*Z^2/r^2).

Also, I would advise that you rewrite your code using a vector for your states instead of individual variables. E.g., define a state vector as follows:

y(1) = x

y(2) = y

y(3) = z

y(4) = xdot

y(5) = ydot

y(6) = zdot

Then implement a derivative function with the following signature:

function dy = myderiv(t,y, other stuff )
% other code
X_2dot = etc.
Y_2dot = etc.
Z_2dot = etc.
dy = [y(4:6);X_2dot;Y_2dot;Z_2dot];
return
end

The other stuff would be the gravity constants, sun & moon positions, etc.

Your downstream code will be much easier to implement if you have your states in a vector like this.

Seems like the comments for rm should read "Earth to Moon" instead of "Spacecraft to Moon"

You are missing a * multiply in these terms:

mu_m*(X-Xm)
mu_m*(Y-Ym)
mu_m*(Z-Zm)

[Blob]

function sol = orbital_motion()

% Define the function that returns the derivative of the state vector

function dydt = ode(t, y)

% Extract the state variables

X = y(1);

Y = y(2);

Z = y(3);

Xdot = y(4);

Ydot = y(5);

Zdot = y(6);

% Calculate the distances and other parameters

r = sqrt(X^2 + Y^2 + Z^2);

rm = sqrt(Xm^2 + Ym^2 + Zm^2);

delta_m = sqrt((X - Xm)^2 + (Y - Ym)^2 + (Z - Zm)^2);

delta_s = sqrt((X - Xs)^2 + (Y - Ys)^2 + (Z - Zs)^2);

% Calculate the derivatives

X2dot = -((mu_e*Y)/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 - Xs)/delta_s^3) ...

- (mu_s*Xs/r^3);

Y2dot = -((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);

Z2dot = -((mu_e*Z)/r^3) * (1 + (J*(a/r)^2)*(3 - 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);

% Return the derivative vector

dydt = [Xdot; Ydot; Zdot; X2dot; Y2dot; Z2dot];

end

[Blob]

@James Tursa @Bjorn Gustavsson sorry for the insitance, but I wrote this, and I jsut want a verification, is everything alright in my code?