compute an elliptical spacraft orbit in matlab

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

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

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