<p>As an engineering student, I am doing a project about Kalman filtering and the mission Apollo 11. Now I tried to simulate the trajectory of the moon and sun for 3 days (the duration in which the spacecraft flight from the earth to the moon). But to my surprise, when I solved the motion equations of the spacecraft, I did not get good results.</p>
<p>How does this trajectory look like?</p>
<p>Here are blocks of the code:</p>
<pre><code>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);

%% date du départ 16 juillet 1969

mois = 7;

jours = [16 17 18];

ann = 1969;

%d = 367*Y - (7*(Y + ((M+9)/12)))/4 + (275*M)/9 + D - 730530

d1 = 367*ann - (7*(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(1) - 
730530;

d2 = 367*ann - (7*(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(2) - 730530;

d3 = 367*ann - (7*(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(3) - 730530;

%% sun position

if t &lt;= 24*60
    d = d1;
elseif t&gt;24*60 &amp; t &lt;= 48*60
    d = d1 + (t - 24*60) / (48*60 - 24*60) * (d2 - d1);
else
    d = d2 + (t - 48*60) / (72*60 - 48*60) * (d3 - d2);
end

w = 282.9404 + 4.70935e-5* d ;%   (longitude of perihelion)

a = 1.000000    ;   %   (mean distance, a.u.)

e = 0.016709 - 1.15e-9* d ; %%  (eccentricity)

M = 356.0470 + 0.9856002585* d ; %  (mean anomaly)

oblecl = 23.4393 - 3.563e-7 * d;

L = w + M ;

E = M + (180/pi) * e * sin(M) * (1 + e * cos(M));

x = cos(E) - e;

y = sin(E) * sqrt(1 - e*e);

r = sqrt(x*x + y*y);

v = atan2( y, x );

lon = v + w;

x = r * cos(lon);

y = r * sin(lon);

z = 0.0;

xequat = 0.881048;

yequat = 0.482098 * cos(oblecl) - 0.0 * sin(oblecl);

zequat = 0.482098 * sin(oblecl) + 0.0 * cos(oblecl);

Xs = xequat*1.496e+11;

Ys = yequat*1.496e+11;

Zs = zequat*1.496e+11;


%% lune position     !!!!

N = 125.1228 - 0.0529538083 * d;   % (Long asc. node)

inc =   5.1454;                          % (Inclination)

w = 318.0634 + 0.1643573223 * d;   % (Arg. of perigee)

a =  60.2666  ;                             % (Mean distance)

e = 0.054900  ;                             % (Eccentricity)

M = 115.3654 + 13.0649929509 * d;    %(Mean anomaly)

% normalisation

M = M + 129*360;

w = w + 360;

E0 = M + (180/pi) * e * sin(M) * (1 + e * cos(M));

E1 = E0 - (E0 - (180/pi) * e * sin(E0) - M) / (1 - e * cos(E0));

E = E1;


x = a * (cos(E) - e);

y = a * sqrt(1 - e*e) * sin(E);

r = sqrt( x*x + y*y );

v = atan2( y, x );

xeclip = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(inc) );

yeclip = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(inc) );

zeclip = r * sin(v+w) * sin(inc);


xeq = xeclip;

yeq = yeclip * cos(oblecl) - zeclip * sin(oblecl);

zeq = yeclip * sin(oblecl) + zeclip * cos(oblecl);

Xm = xeq*1.496e+11 ;

Ym = yeq*1.496e+11;

Zm = zeq*1.496e+11;

%% constantes

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.

% Calculate the distances and other parameters

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

rs = sqrt(Xs^2 + Ys^2 + Zs^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);

Y_s = Ys;

Z_s = Zs;


% 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/rs^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/rs^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/rs^3);


% Return the derivative vector

%X2dot = interp1(Xdot,X2dot,t);

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

end

Then we have the ODE45 to resolve the equation

close all
clear all; clc;

%% Set the initial conditions

    
    X0 =100 + 6.37826e6;

    Y0 = 100 + 6.37826e6;

    Z0 = 100 + 6.37826e6; 

     Xdot0 =  -2.0225e+14;    %Xm/72*3600, repère géocentrique %8.9919e+10  ;

     Ydot0 =  -1.8023e+14;%-1.0908e+11;

     Zdot0 =  3.4945e+14; %-4.8960e+10;

    y0 = [X0; Y0; Z0; Xdot0; Ydot0; Zdot0];
    
    tspan = 0:1:72*60*60;
    
[t, y] = ode45(@ode, tspan, y0);
    
  

plot3(y(:,1),y(:,2),y(:,3))

title('Solution of the Motion Equation with ODE45');
</code></pre>

CC BY-SA 4.0
2023-05-06T15:49:42.700Z
A few comments.  Firstly, are you aware of [space.se] it is about human exploration of space and may be more appropriate.  Secondly we probably can't help in debugging code.  You'll need to break down the code into testable sections and test each part.  However, There are folks here who can discuss the mathematics of orbital trajectories. Thirdly, you didn't attach the code.
James K
CC BY-SA 4.0
2023-05-06T15:49:57.127Z
You can compare your results to JPL Horizons: https://ssd.jpl.nasa.gov/horizons/app.html#/
Greg Miller
CC BY-SA 4.0
2023-05-06T18:19:15.913Z
The code works! it is the calculation of the initial conditions that I am not sure of
Blobmou
CC BY-SA 4.0
2023-05-06T18:21:47.433Z
@GregMiller I calculated the Moon and Sun positions using the equations in this website https://stjarnhimlen.se/comp/tutorial.html#7
Blobmou
CC BY-SA 4.0
2023-05-06T22:15:27.203Z
As Uhoh said "please debug" doesn't often get answered. The linked page for calculating the position of the moon looks sound, (though remember you are really finding the position of the Earth, not the position of the sun) But your code is one big unstructured mass of assignments.  If it isn't doing what you expect it to do, you need to split it up into testable functions, and test them
James K
CC BY-SA 4.0
2023-05-06T22:24:37.410Z
Ok, but how do I verify of my calculations are correct (initial positions, sun and moon position)?
Blobmou
CC BY-SA 4.0
2023-05-06T22:31:19.777Z
See the comment above by Greg, you can compare your positions to those calculated by the high accuracy formulae used by the Horizons system
James K
CC BY-SA 4.0
2023-05-07T08:47:00.273Z
Thanks for the edit! Hopefully the drive-by downvoting will stop :-)
uhoh
CC BY-SA 4.0
2023-05-07T17:10:06.100Z
@uhoh I think I found what is wrong, the initial velocity of the spacecraft, does anyone know how do I calculate it to reach the moon?
Blobmou
CC BY-SA 4.0
2023-05-07T19:59:14.020Z
@Blobmou based on [Greg Miller's comment](https://astronomy.stackexchange.com/questions/53623/difficulty-simulating-trajectory-of-moon-and-sun-for-apollo-11-mission?noredirect=1#comment119816_53623) I looked in [Horizons](https://ssd.jpl.nasa.gov/horizons/app.html#/) at the trajectory of the first four days of Apollo 12's trajectory for example, see https://i.sstatic.net/n2hTv.png and https://i.sstatic.net/fZccK.png where state vectors (x, y, z, vx, vy, vz) can be obtained.
uhoh
CC BY-SA 4.0
2023-05-09T19:18:56.333Z
You can use my script at the end of [this answer](https://astronomy.stackexchange.com/a/49823/16685) to view these trajectories in 3D. But you won't see much detail if you include the Sun. The ID for Apollo 11 is [-399110](https://ssd.jpl.nasa.gov/api/horizons.api?format=text&MAKE_EPHEM=NO&COMMAND=-399110)
PM 2Ring
<orbit><n-body-simulations><modeling>
Difficulty simulating trajectory of moon and sun for Apollo 11 mission
Blobmou
# Difficulty simulating trajectory of moon and sun for Apollo 11 mission

(Asked by: username_0 on 2023-05-06)

As an engineering student, I am doing a project about Kalman filtering and the mission Apollo 11. Now I tried to simulate the trajectory of the moon and sun for 3 days (the duration in which the spacecraft flight from the earth to the moon). But to my surprise, when I solved the motion equations of the spacecraft, I did not get good results.
How does this trajectory look like?
Here are blocks of the code:
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);

%% date du départ 16 juillet 1969

mois = 7;

jours = [16 17 18];

ann = 1969;

%d = 367*Y - (7*(Y + ((M+9)/12)))/4 + (275*M)/9 + D - 730530

d1 = 367*ann - (7*(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(1) - 
730530;

d2 = 367*ann - (7*(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(2) - 730530;

d3 = 367*ann - (7*(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(3) - 730530;

%% sun position

if t <= 24*60
    d = d1;
elseif t>24*60 & t <= 48*60
    d = d1 + (t - 24*60) / (48*60 - 24*60) * (d2 - d1);
else
    d = d2 + (t - 48*60) / (72*60 - 48*60) * (d3 - d2);
end

w = 282.9404 + 4.70935e-5* d ;%   (longitude of perihelion)

a = 1.000000    ;   %   (mean distance, a.u.)

e = 0.016709 - 1.15e-9* d ; %%  (eccentricity)

M = 356.0470 + 0.9856002585* d ; %  (mean anomaly)

oblecl = 23.4393 - 3.563e-7 * d;

L = w + M ;

E = M + (180/pi) * e * sin(M) * (1 + e * cos(M));

x = cos(E) - e;

y = sin(E) * sqrt(1 - e*e);

r = sqrt(x*x + y*y);

v = atan2( y, x );

lon = v + w;

x = r * cos(lon);

y = r * sin(lon);

z = 0.0;

xequat = 0.881048;

yequat = 0.482098 * cos(oblecl) - 0.0 * sin(oblecl);

zequat = 0.482098 * sin(oblecl) + 0.0 * cos(oblecl);

Xs = xequat*1.496e+11;

Ys = yequat*1.496e+11;

Zs = zequat*1.496e+11;


%% lune position     !!!!

N = 125.1228 - 0.0529538083 * d;   % (Long asc. node)

inc =   5.1454;                          % (Inclination)

w = 318.0634 + 0.1643573223 * d;   % (Arg. of perigee)

a =  60.2666  ;                             % (Mean distance)

e = 0.054900  ;                             % (Eccentricity)

M = 115.3654 + 13.0649929509 * d;    %(Mean anomaly)

% normalisation

M = M + 129*360;

w = w + 360;

E0 = M + (180/pi) * e * sin(M) * (1 + e * cos(M));

E1 = E0 - (E0 - (180/pi) * e * sin(E0) - M) / (1 - e * cos(E0));

E = E1;


x = a * (cos(E) - e);

y = a * sqrt(1 - e*e) * sin(E);

r = sqrt( x*x + y*y );

v = atan2( y, x );

xeclip = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(inc) );

yeclip = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(inc) );

zeclip = r * sin(v+w) * sin(inc);


xeq = xeclip;

yeq = yeclip * cos(oblecl) - zeclip * sin(oblecl);

zeq = yeclip * sin(oblecl) + zeclip * cos(oblecl);

Xm = xeq*1.496e+11 ;

Ym = yeq*1.496e+11;

Zm = zeq*1.496e+11;

%% constantes

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.

% Calculate the distances and other parameters

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

rs = sqrt(Xs^2 + Ys^2 + Zs^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);

Y_s = Ys;

Z_s = Zs;


% 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/rs^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/rs^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/rs^3);


% Return the derivative vector

%X2dot = interp1(Xdot,X2dot,t);

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

end

Then we have the ODE45 to resolve the equation

close all
clear all; clc;

%% Set the initial conditions

    
    X0 =100 + 6.37826e6;

    Y0 = 100 + 6.37826e6;

    Z0 = 100 + 6.37826e6; 

     Xdot0 =  -2.0225e+14;    %Xm/72*3600, repère géocentrique %8.9919e+10  ;

     Ydot0 =  -1.8023e+14;%-1.0908e+11;

     Zdot0 =  3.4945e+14; %-4.8960e+10;

    y0 = [X0; Y0; Z0; Xdot0; Ydot0; Zdot0];
    
    tspan = 0:1:72*60*60;
    
[t, y] = ode45(@ode, tspan, y0);
    
  

plot3(y(:,1),y(:,2),y(:,3))

title('Solution of the Motion Equation with ODE45');

**Tags:** orbit, n-body-simulations, modeling

> **Comments:**
> **username_1**: A few comments.  Firstly, are you aware of [space.se] it is about human exploration of space and may be more appropriate.  Secondly we probably can't help in debugging code.  You'll need to break down the code into testable sections and test each part.  However, There are folks here who can discuss the mathematics of orbital trajectories. Thirdly, you didn't attach the code.  
> **username_2**: You can compare your results to JPL Horizons: https://ssd.jpl.nasa.gov/horizons/app.html#/  
> **username_0**: The code works! it is the calculation of the initial conditions that I am not sure of  
> **username_0**: @GregMiller I calculated the Moon and Sun positions using the equations in this website https://stjarnhimlen.se/comp/tutorial.html#7  
> **username_1**: As Uhoh said "please debug" doesn't often get answered. The linked page for calculating the position of the moon looks sound, (though remember you are really finding the position of the Earth, not the position of the sun) But your code is one big unstructured mass of assignments.  If it isn't doing what you expect it to do, you need to split it up into testable functions, and test them  
> **username_0**: Ok, but how do I verify of my calculations are correct (initial positions, sun and moon position)?  
> **username_1**: See the comment above by Greg, you can compare your positions to those calculated by the high accuracy formulae used by the Horizons system  
> **username_3**: Thanks for the edit! Hopefully the drive-by downvoting will stop :-)  
> **username_0**: @username_3 I think I found what is wrong, the initial velocity of the spacecraft, does anyone know how do I calculate it to reach the moon?  
> **username_3**: @username_0 based on [username_2's comment](https://astronomy.stackexchange.com/questions/53623/difficulty-simulating-trajectory-of-moon-and-sun-for-apollo-11-mission?noredirect=1#comment119816_53623) I looked in [Horizons](https://ssd.jpl.nasa.gov/horizons/app.html#/) at the trajectory of the first four days of Apollo 12's trajectory for example, see https://i.sstatic.net/n2hTv.png and https://i.sstatic.net/fZccK.png where state vectors (x, y, z, vx, vy, vz) can be obtained.  
> **username_4**: You can use my script at the end of [this answer](https://astronomy.stackexchange.com/a/49823/16685) to view these trajectories in 3D. But you won't see much detail if you include the Sun. The ID for Apollo 11 is [-399110](https://ssd.jpl.nasa.gov/api/horizons.api?format=text&MAKE_EPHEM=NO&COMMAND=-399110)  

---

## 0 answers

