Borrar filtros
Borrar filtros

Motion equation of a space engin in low earth orbit

3 visualizaciones (últimos 30 días)
Blob
Blob el 3 de Mzo. de 2023
Comentada: Blob el 30 de Abr. de 2023
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);
  10 comentarios
Bjorn Gustavsson
Bjorn Gustavsson el 4 de Abr. de 2023
@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
Blob el 4 de Abr. de 2023
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.

Iniciar sesión para comentar.

Respuestas (1)

James Tursa
James Tursa el 31 de Mzo. de 2023
Editada: James Tursa el 31 de Mzo. de 2023
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)
  2 comentarios
Blob
Blob el 30 de Abr. de 2023
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
Blob el 30 de Abr. de 2023
@James Tursa @Bjorn Gustavsson sorry for the insitance, but I wrote this, and I jsut want a verification, is everything alright in my code?

Iniciar sesión para comentar.

Categorías

Más información sobre Gravitation, Cosmology & Astrophysics en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by