compute an elliptical spacraft orbit in matlab

To compute an elliptical spacecraft orbit in MATLAB, we need to follow these steps:

  1. Define the gravitational constant mu of the central body around which the spacecraft is orbiting.
  2. Define the semi-major axis a and the eccentricity e of the spacecraft's orbit.
  3. Calculate the semi-minor axis b of the spacecraft's orbit using the following formula: b = a * sqrt(1 - e^2).
  4. Define the inclination i, the longitude of ascending node Omega, and the argument of periapsis omega of the spacecraft's orbit.
  5. Define the mean anomaly M of the spacecraft's orbit.
  6. Compute the eccentric anomaly E of the spacecraft's orbit using the following iterative formula: E = M + e * sin(E).
  7. Compute the true anomaly nu of the spacecraft's orbit using the following formula: nu = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2)).
  8. Compute the position and velocity vectors of the spacecraft in the perifocal frame using the following formulas:
main.m
r_pqw = [a * (cos(E) - e); a * sqrt(1 - e^2) * sin(E); 0];
v_pqw = [-(sqrt(mu / a) / (1 - e * cos(E))) * sin(E); sqrt(mu / a) / (1 - e * cos(E)) * (sqrt(1 - e^2) * cos(E)); 0];
177 chars
3 lines
  1. Transform the position and velocity vectors from the perifocal frame to the geocentric equatorial frame using the following transformation matrix:
main.m
Q_pqw_to_eci = [cos(omega) * cos(Omega) - sin(omega) * sin(Omega) * cos(i), -cos(omega) * sin(Omega) - sin(omega) * cos(Omega) * cos(i), sin(omega) * sin(i); sin(omega) * cos(Omega) + cos(omega) * sin(Omega) * cos(i), -sin(omega) * sin(Omega) + cos(omega) * cos(Omega) * cos(i), -cos(omega) * sin(i); sin(i) * sin(Omega), sin(i) * cos(Omega), cos(i)];
r_eci = Q_pqw_to_eci * r_pqw;
v_eci = Q_pqw_to_eci * v_pqw;
412 chars
4 lines

Here's the MATLAB code that implements these steps:

main.m
% Define constants and parameters
mu = 3.986e5;                 % Gravitational constant of Earth [km^3/s^2]
a = 7000;                     % Semi-major axis of the spacecraft's orbit [km]
e = 0.2;                      % Eccentricity of the spacecraft's orbit
i = deg2rad(30);              % Inclination of the spacecraft's orbit [rad]
Omega = deg2rad(45);          % Longitude of ascending node of the spacecraft's orbit [rad]
omega = deg2rad(60);          % Argument of periapsis of the spacecraft's orbit [rad]
M = deg2rad(90);              % Mean anomaly of the spacecraft's orbit [rad]

% Compute the semi-minor axis of the spacecraft's orbit
b = a * sqrt(1 - e^2);

% Compute the eccentric anomaly of the spacecraft's orbit
E = M;
max_iterations = 1000;        % Maximum number of iterations for convergence
tolerance = 1e-8;             % Convergence tolerance
for i = 1:max_iterations
    f = E - e * sin(E) - M;
    df = 1 - e * cos(E);
    dE = -f / df;
    E = E + dE;
    if abs(dE) < tolerance
        break;
    end
end

% Compute the true anomaly of the spacecraft's orbit
nu = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2));

% Compute the position and velocity vectors of the spacecraft in the perifocal frame
r_pqw = [a * (cos(E) - e); b * sin(E); 0];
v_pqw = [-(sqrt(mu / a) / (1 - e * cos(E))) * sin(E); sqrt(mu / a) / (1 - e * cos(E)) * (sqrt(1 - e^2) * cos(E)); 0];

% Compute the transformation matrix from the perifocal frame to the geocentric equatorial frame
Q_pqw_to_eci = [cos(omega) * cos(Omega) - sin(omega) * sin(Omega) * cos(i), -cos(omega) * sin(Omega) - sin(omega) * cos(Omega) * cos(i), sin(omega) * sin(i); sin(omega) * cos(Omega) + cos(omega) * sin(Omega) * cos(i), -sin(omega) * sin(Omega) + cos(omega) * cos(Omega) * cos(i), -cos(omega) * sin(i); sin(i) * sin(Omega), sin(i) * cos(Omega), cos(i)];

% Transform the position and velocity vectors from the perifocal frame to the geocentric equatorial frame
r_eci = Q_pqw_to_eci * r_pqw;
v_eci = Q_pqw_to_eci * v_pqw;
2003 chars
40 lines

Note that this code assumes a circular central body, which is a good approximation for orbits around Earth. If the central body is not spherical or has significant oblateness, more complex calculations may be necessary.

gistlibby LogSnag