How to build an orbit propagator?

Question

I am trying to build my own satellite orbit prapagator in MATLAB with RK4 integrator and force models added to the propagator to give the perturbed orbit in ECI frame. I am a bit confused and just wanted to clarify if my understanding is right here - Suppose we have to plot the satellite's orbit in a time frame of 6000 seconds. The code first converts orbital elements to the initial position and velocity vector in ECI frame using Kepler equations. Then, based on the integrator chosen (eg: RK4), the integrator propagates the orbit from the initial position and velocity vectors for the time frame of 6000 seconds.

True anomaly is the only orbital element that iterates with time so if true anomaly changes every second, then doesn't the initial position and velocity vectors change every second, making it hard for the integrator as it would propagate a different orbit in each of those 6000 seconds due to different initial positions?

I apologise for the long question as I am trying to get the concepts of an orbit propagator clarified!

Edit:

I think the Python files helped but I wanted to clarify that I wish to add a force model (aPert) which includes drag, SRP, third body and spherical harmonics gravity to my rk4 propagator, like I have done so in the code below. It changes with every second as it is dependant on time so when I use it in the rk4 integration, since the acceleration is already propagated from an initial value, I only need to provide the initial aPert and then it will be propagated using the rk4 propagator right?

function dot = orbitfunc(t,X)
    mu = 398600435507000;    % Gravitational constant (m^3/s^2)
    r = X(1:3);                % Position (m)
    v = X(4:6);                % Velocity (m/s^2)
    dr = v;
    dv = (-mu/(norm(r))^3)*r + aPert; % Newton's law of gravity
    dot = [dr; dv];  
end

function [Xf,K] = rk4orbit(t,X,h)
% Numerical Solution, Runge-Kutta 4th Order
    k1 = orbitfunc(t,X);
    k2 = orbitfunc(t+h/2,X+k1*h/2);
    k3 = orbitfunc(t+h/2,X+k2*h/2);
    k4 = orbitfunc(t+h,X+k3*h);
    K = (k1+2*k2+2*k3+k4)/6; % Acceleration
% Step forward in time
    Xf = X+(h/6)*(k1+2*k2+2*k3+k4);
end

Answer 1

Here's my take on your question "How to build an orbit propagator?"

The code first converts orbital elements to the initial position and velocity vector in ECI frame using Kepler equations. Then, based on the integrator chosen (eg: RK4), the integrator propagates the orbit from the initial position and velocity vectors for the time frame of 6000 seconds.

That's great. If you have the position and velocity in cartesian coordinates (i.e. $(x, y, z, v_x, v_y, v_z)$ at some time say $t=0$ then you have an initial state vector.

At this point, please completely abandon Kepler, his equations and orbits and think only in terms of integrating the orbit in cartesian coordinates. Do not do any further conversions.

True anomaly is the only orbital element that iterates with time so if true anomaly changes every second, then doesn't the initial position and velocity vectors change every second, making it hard for the integrator as it would propagate a different orbit in each of those 6000 seconds due to different initial positions?

Forget Kepler's true anomaly, and just think of time as your independent variable.

I don't know if you are writing your own solver (e.g RK4) or using one that MATLAB provides. It's extremely educational and gratifying to get a manually written solver running, you learn a lot! Once you do, then feel confident to use a library function. Then you can use any ODE solver you like.

I found this Python script with hand-written RK4 and RK4(5) solvers in an old folder on my laptop. I can not speak to it's correctness but I think it gives you a basic idea.

• https://pastebin.com/35P0ERqZ

While @DavidHammen's answer is of course right as always, the world will not end if you use any solver you like to learn how to do it.

You'll find that step sizes can be much bigger than 1 second for good results, maybe 10 seconds or longer can still get you within a few tens of meters of accuracy after one orbit.

What's cool about RK4(5) (Runge–Kutta–Fehlberg method) in that script is that it is variable size. It uses a 5th order to estimate the error and then adjusts the internal step size to try to match the level of acceptable error you specify.

All of this helps you understand the helpfulness of using a variable step size ODE solver, especially when you move from a circular orbit to a highly elliptical orbit where steps can be much larger at apoapsis than at periapsis.

Have fun, and consider switching to python (get the Anaconda installation) or even Go which will be incredibly liberating as you can get help all over the place and it's FREE forever! :-)

Answer 2

RK4 is far from the best numerical integration techniques to use for orbit propagation. This is particularly the case for non-circular orbits. The key problem with using RK4 is that it discards geometry. While any numerical integrator for a first order ordinary differential equation (ODE) can be adapted to an $n^{\text{th}}$ order ODE by creating auxiliary variables, doing so inevitably discards geometry. This is exactly what using RK4 for orbit propagation does. The relevant differential equations are $$\begin{aligned} \frac{d\mathbf x(t)}{dt} &= \mathbf v(t) \\ \frac{d\mathbf v(t)}{dt} &= \mathbf a(\mathbf x(t), \mathbf v(t), t) \end{aligned}$$ where $\mathbf x(t)$ is the instantaneous position vector, $\mathbf v(t)$ is the instantaneous velocity vector, and $\mathbf a(\mathbf x(t), \mathbf v(t), t)$ is the instantaneous acceleration vector expressed as a function of instantaneous position, velocity, and time. That this is inherently a second order ODE is geometry, and treating it as a first order ODE by creating a six vector to be integrated throws geometry under the bus. RK4 does not conserve total energy or angular momentum.

Another issue with RK4 is the fixed step size. I'll start with a circular orbit subject to an inverse square central force problem (e.g., Newtonian point source gravitation). The object doesn't move at all if the step size is ridiculously small and finite precision arithmetic is used. This is because $1+10^{-16}$ is identically equal to one using IEEE double precision (64 bit) arithmetic. As the step size is increased, the computation time decreases as does the error -- to a point. At some point with increased step size, the truncation errors inherent to the integration technique itself begin to dominate over the decrease in finite precision error. Different numerical integration techniques have different inflection points, but they all suffer from these two competing sources of error.

Adaptive techniques regularly adjust the integration step size so as to attempt to stay near that sweet spot where the error is minimized. RK4 is not adaptive. The issue with non-circular orbits is that the ideal step size varies over the course of an orbit. The ideal step size near periapsis is too small for points near apoapsis, resulting in overly large finite precision errors near apoapsis, while the ideal step size near apoapsis is too large for points near periapsis, resulting in overly large truncation errors there.

True anomaly is the only orbital element that iterates with time so if true anomaly changes every second, then doesn't the initial position and velocity vectors change every second, making it hard for the integrator as it would propagate a different orbit in each of those 6000 seconds due to different initial positions?

True anomaly is the only Keplerian orbital element that changes over time if (and only if) the problem in question is an inverse square central force problem such as Newtonian gravitational attraction to a single point source. All of the Keplerian elements are subject to change when there is more than one gravitating body is present, or when the central body does not have a spherical mass distribution, or when non-gravitational perturbing forces are present, or when non-Newtonian concerns (e.g., general relativity) are brought into play. These real-world issues throw Keplerian orbits under the real-world bus, and this is why numerical integration techniques are used for real-world orbit propagation.

With regard to KK4 specifically, there are multiple (an infinite number!) of RK4 integrators, many of which have fourth order accuracy. The Butcher tableau for the canonical RK4 integrator is

$$\begin{array} {c|cccc} 0\\ \frac{1}{2} & \frac{1}{2}\\ \frac{1}{2} &0 &\frac{1}{2} \\ 1& 0& 0& 1\\ \hline & \frac{1}{6} &\frac{1}{3} &\frac{1}{3} &\frac{1}{6} \end{array}$$

What this means is that one first caches the state (position and velocity) at the start of the interval as these will be used repeatedly: $$\begin{aligned} t_0 &\equiv t \\ \mathbf x_0 &\equiv \mathbf x(t) \\ \mathbf v_0 &\equiv \mathbf v(t) \\ \mathbf a_0 &= \mathbf a(\mathbf x_0, \mathbf v_0, t_0) \end{aligned}$$

Note that $\mathbf a_0$ is calculated but not used in this step. It will be used in subsequent steps. Note also that the acceleration may depend on position (e.g., gravitation), velocity (e.g., drag), and time (e.g., Earth orientation). Next one calculates a guess at the state at the midpoint of the step using an Euler step based on state at the start of the interval: $$\begin{aligned} \mathbf t_1 &= t_0 + \frac{\Delta t}2 \\ \mathbf x_1 &= \mathbf x_0 + \frac{\Delta t}2 \,\mathbf v_0 \\ \mathbf v_1 &= \mathbf v_0 + \frac{\Delta t}2 \,\mathbf a_0 \\ \mathbf a_1 &= \mathbf a(\mathbf x_1, \mathbf v_1, t_1) \end{aligned}$$ Next one performs another Euler step, again to the midpoint of the interval but using the state at the first intermediate step to calculate derivatives: $$\begin{aligned} \mathbf t_2 &= t_0 + \frac{\Delta t}2 \\ \mathbf x_2 &= \mathbf x_0 + \frac{\Delta t}2 \,\mathbf v_1 \\ \mathbf v_2 &= \mathbf v_0 + \frac{\Delta t}2 \,\mathbf a_1 \\ \mathbf a_2 &= \mathbf a(\mathbf x_2, \mathbf v_2, t_2) \end{aligned}$$ Note that $t_2 = t_1$, so there's no need to cache this. Next one performs yet another Euler step, this time from the start of the interval to the end of the interval, using the values from the previous step: $$\begin{aligned} \mathbf t_3 &= t_0 + \Delta t \\ \mathbf x_3 &= \mathbf x_0 + \Delta t \,\mathbf v_2 \\ \mathbf v_3 &= \mathbf v_0 + \Delta t \,\mathbf a_2 \\ \mathbf a_3 &= \mathbf a(\mathbf x_3, \mathbf v_3, t_3) \end{aligned}$$

Finally, one calculates a guess at the state at the end of the interval again via an Euler step but this time using weighted averages for the derivatives:

$$\begin{aligned} \mathbf x(t+\Delta t) &= \mathbf x_0 + \Delta t \left(\frac16 \mathbf v_0 + \frac13\mathbf v_1 + \frac13\mathbf v_2 + \frac16\mathbf v_3\right) \\ \mathbf v(t+\Delta t) &= \mathbf v_0 + \Delta t \left(\frac16 \mathbf a_0 + \frac13\mathbf a_1 + \frac13\mathbf a_2 + \frac16\mathbf a_3\right) \end{aligned}$$

An arguably better (in terms of error) but worse (in terms of computation time) fourth order Runge-Kutta integrator has the following Butcher tableau: $$\begin{array} {c|cccc} 0\\ \frac{1}{3} & \frac{1}{3}\\ \frac{2}{3} &\frac{-1}{3} &1 \\ 1& 1& -1& 1\\ \hline & \frac{1}{8} &\frac{3}{8} &\frac{3}{8} &\frac{1}{8} \end{array}$$ This is not widely used. When one says they are using RK4 they almost always mean the canonical RK4 motivated by Simpson's rule.