Can anyone please explain me how to plot for the attached algorithm?

The attached text is in a orbital mechanics problem. I got the answer in two methods one is using the universal formulation (based on iterations) and the other one is using ODE45 method. I got graph for the ODE45 method. However, I don't how to plot for the universal formulation and arrive with the same graph I got for the ODE45 method. Can anyone kindly help me on how to get the same graph with the universal formulation algorithm? Thanks

The below attached graphs are the graphs I got from the ODE45 method

(figure a) Plot the trajectory of the spacecraft from the initial position and for the next 3 hours. In the plot, indicate clearly the initial and final positions (both solutions).

(figure b) Plot the magnitude of the radius as a function of time (both solutions).

Answer 1

[Altaïr]

Hey @Zev M,

Great job implementing both orbital propagation methods! The custom universal formulation using Kepler's equation with Stumpff functions and the numerical ODE45 approach seem to successfully calculate the final positions. The difficulty in generating trajectory and radius plots for the universal method likely stems from how the code currently computes only the endpoint state (at t=3 hours), whereas plotting requires intermediate states throughout the propagation period. Unlike ODE45—which automatically returns states at multiple time steps—the analytical approach needs explicit time-stepping to capture the full trajectory.

To resolve this, consider modifying the universal formulation code with a time-stepping loop that computes positions at regular intervals. Here's a suggested adjustment:

% Universal Formulation Plot Generator
global mu
mu = 398600; 
R0 = [7000; -12124; 3000]; 
V0 = [2.6679; 4.6210; 1];  
t_total = 3600*3;  % 3 hours

% Time-stepping parameters
dt = 60;  % 60-second intervals (adjust as needed)
t_steps = 0:dt:t_total;
positions = zeros(length(t_steps), 3);

% Compute states at each time step
for i = 1:length(t_steps)
    [R, ~] = rv_from_r0v0(R0, V0, t_steps(i));
    positions(i, :) = R';
end

% Calculate radii and plot (use same plotting code as ODE45 method)
radii = vecnorm(positions, 2, 2);

Alternatively, for MATLAB R2024b or later, the propagateOrbit function offers built-in trajectory generation. Control the propagation model via the PropModel attribute:

https://www.mathworks.com/help/aerotbx/ug/propagateorbit.html#mw_9467fd42-d4fa-4e03-ac74-703b12a0207b

For earlier MATLAB versions, an alternative is creating a Satellite object and configuring its OrbitPropagator property. Use the following command to read more about the Satellite object:

web(fullfile(docroot, 'aerotbx/ug/satellite.html'))