Simulating orbital dynamics may seem challenging but with a bit of thought, you can easily start with simplified systems and then extend them as needed to get the accuracy that you require. In this article, we will walk through the math that can be used to describe the motion of a simplified system and we will write code that can simulate the motion of a satellite around a planet. Then we will use that simulation to visualize a few types of orbits.
Basic Simulation
Computer simulations of physical systems are used to evaluate things that have not yet been built, as well as situations and conditions in which you might not be able to test an entire system. These simulations are based on states such as position, velocity, orientation,tempertaure etc. These states are used to represent the state of the object being simulated. Time-domain simulations are most common for physical systems as you want to determine how a system will perform over a period of time when exposed to specific initial conditions or environments. These systems can be described by sets of differential equations that show how the system changes for each timestep. As the common phrase goes “All models are wrong, but some are useful”, all simulations are approximate, but some are good enough.
Orbital Dynamics
The most basic simplified orbital dynamic system is a single object orbiting around another object, where the mass difference between the two is so small that the mass of the second object is negligible. This is important because the only force that we will be modeling is gravity. There is no air resistance, and for starters, we will not be modeling any orientation states or propulsion of the satellite.
Point Masses
We can think of our two objects as a mass distribution that is centered around a point. We make the simplification that our two objects are infinitely small points where all of the mass in concentrated at that point. This means that we can plot the location of of our objects as coordinates in 3D space in the form of a vector.
Gravity
The basic gravity model is a function of altitude that is determined by the following equation:
$$F_g=\frac{Gm_em_s}{r^2}$$
where \(r\) is the radius of the earth plus the altitude, \(m_e\) is the mass of the earth \(m_s\) is the mass of the satellite, and \(G\) is the gravitational constant. This model is typically used up to an altitude of 400km, but for our purposes we will use it for altitudes higher than that.
As you can see, the gravity between two objects depends on their mass and the distance between them. In a real system the satellite pulls on the earth with the same force that the earth pulls on the satellite. Here is where our first simplification come in, we will assume that the difference in mass between the earth and our satellite is so large that we don’t need to account for the mass of the satellite. If we say that our satellite has a mass of 1000kg, given the Earth’s mass of 5973600000000000000000000 kg (which we call \(5.9736e24\) kg) then the mass of the earth is 21 orders of magnitude larger than the satellite. Using Newton’s second Law, the equation for a constant mass object is:
$$F=ma$$
We can see that the acceleration of the satellite due to the earth is
$$a_{satellite}=F\frac{1}{1000}$$
And the acceleration of the Earth due to the satellite is
$$a_{earth}=F\frac{1}{5.9736e24}$$
Therefore the acceleration of the earth due to the satellite is around 28 orders of magnitude smaller than the acceleration of the satellite due to the earth.
By using this simplification we can use the point mass of the earth as a reference frame for our simulation this allows us to represent the position as a single vector in 3D space, the position of the satellite with respect to the point mass of the earth.
Differential Equations
For our simulation we want to derive the differential equations of motion for our model. We can represent the motion of the satellite using a position and a velocity vector. These are then discretized into timesteps that we can implement in a program.
Numerical Integration methods are used to solve initial value problems of differential equations. Integration methods make assumptions about how the signal is changing between the discrete points therefore the higher the order of the integration method, the more accurate the results are. Different integration methods may use different numbers of function evaluations to step from a value at one time step to a value at a second time step. Using numerical integration has the potential to introduce instability and errors, but it is usually the only way to solve very complicated systems of equations.
The Runge-Kutta methods are numerical integration methods that use a number of different slopes to calculate the new state. They fall under the category of single-step methods and are the workhorse of numerical integrations. The second and fourth-order schemes are widely used in simulation.
The classic Runge-Kuttta method (RK4) is one of the most popular of the Runge-Kutta family of solvers. It was developed around the year 1900 by the mathematicians Carl Runge and Martin Wilhelm Kutta. The main advantage of the fourth-order Runge-Kutta method is that it is more accurate than the forward Euler method. It has been used for numerical integrations for more than 100 years and has been used in control systems for decades. The RK4 method cancels out errors out to order \(\Delta t^5\) because an n-step RK method has (n+1) order accuracy. The equation is:
$$\begin{matrix}
x_{k+1}=x_k+\frac{\Delta t}{6}(f_1+2f_2+2f_3+f_4)\\ \text{Where} \\
f_1 = f(x_k, t_k) \\
f_2 = f(x_k+\frac{\Delta t}{2}f_1,t_k+\frac{\Delta t}{2}) \\
f_3 = f(x_k+\frac{\Delta t}{2}f_2, t_k+\frac{\Delta t}{2})\\
f_4 = f(x_k+\Delta f_3, t_k+\Delta t)
\end{matrix}
$$
Vector Operations
Vectors are numeric quantities with a magnitude and a direction. The magnitude of a vector is a scalar. Vectors are often underlined or written in bold font for notation and can be represented as matrices in programs.
The vector length (also known as the L2 Norm) is the square root of the sum of the squared components. We use this formula to turn a vector into a unit vector.
A Normalized vector is a vector of length 1 also called a unit vector. It is found by scaling a vector by its length.
To model our system we need a position vector and a velocity vector. Both of these are in 3D space, therefore they each have 3 components. We then also need equations that describe how our system is changing over time. In our model, gravity will always pull towards the center of the Earth. This means that that we can use the negative of the normalized position vector to represent the unit direction of gravity as it has a range between zero and one.
$$\hat{s}=\frac{s}{|s|}$$
Next we have to calculate the acceleration due to gravity at the current position. This is done through the equation that we discussed previously. By multiplying the equation for the acceleration due to gravity by the equation for the unit vector of gravity we get the equation for the gravitational acceleration of our satellite.
$$\vec{a_g}=-\hat{s}\frac{\mu_{uearth}}{|\vec{s}|^2}$$
Acceleration is the change in velocity and velocity is the change in position. This gives us the two differential equations that describe the change in state of our system given the current state. The position vector is changed by the velocity vector, and the velocity vector is changed by the acceleration vector. Our system equations are
$$\begin{matrix}
\dot{\vec{v}}=-\hat{s}\frac{\mu_{Earth}}{|\vec{s}|^2} \\
\dot{\vec{s}}=\vec{v}
\end{matrix}$$
In two dimensions we can visualize this as the velocity vector being rotated by the acceleration due to gravity, and the position vector being rotated by the velocity vector.
Finally we need to have some represntation of the surface of the earth as a boundary. Because our model uses point masses, we could represent an orbit that is 100m away from the center of mass of the earth. This is not actually possible so we need to have a way to check and see if the satellite has passed through the surface and is actually moving underground. We can do this by taking the magnitude of the position vector and comparing it with the radius of the earth. Because our earth model is a perfect sphere, if the magnitude of the position vector of the satellite is less than the surface of the earth then we have hit the planet and we can stop the simulation.
Code
For the code for this simulation we will use the Python language with a few additional packages. The reason to use python is that it has a simple syntax that is easy for people without much experience with programming languages to understand. I will not go into any more detail about the Python language as there are countless resources devoted to that. The reason that we are using the Pytorch library is that as a machine learning library has some useful libraries for matrix operations. It also makes it very easy to accelerate any algorithms using matrices on a GPU, but designing algorithms for hardware acceleration will be left for a future article. The matplotlib library is used for visualization and the numpy library is used for some numerical calcuations.
For the header we need to import the Pytorch library, the Matplotlib library, and the numpy library
We can then define a function for the magnitude (L2 norm) of our vector
Our next function will describe the equations of motion of our system as we defined above.
This takes the position and velocity state vectors as input, and GM is the same as \(u_{earth}\) in the equations above. The torch.mul() function is a scalar multiplication of the unit vector \(-\hat{s}\) by the acceleration due to gravity.
Our ODE solver is an implementation of the 4th order Runge-Kutta method for two states. By separating out the solver from our model, as long as we keep the same inputs and outputs we can easily change our model (or add mutiple simultaneous models) without re-writing most of the simulation.
Visualization
The next functions that we need will be used visualization of the results of our simulation in both 2D and 3D. The first function takes the set of position vectors and plots the orbit in 3D with a reference point representing the center of the earth. It also plots a marker if there is a collision with the Earth.
Our next function plots the satellite orbit as well as the surface of the earth when viewed from above. It will also save an image of the plot with a time stemp.
Main Simulation Loop
Finally, we get to the meat of our simulation. We need to initialize our state vectors to whatever state we want for our simulation. In this case, we will use a 200km orbit starting on the x-axis with an orbital velocity of 7.8km/s in the y-axis direction.
$$\begin{matrix}\vec{s}_0=\begin{bmatrix}
200e3 + \text{Radius of the Earth} \\
0 \\
0 \\
\end{bmatrix} \\
\vec{v}_0=\begin{bmatrix}
0 \\
7.8e3 \\
0 \\
\end{bmatrix}
\end{matrix}$$
We also can set our timestep to 1 minute and the total simulation time to 100 minutes.
For our main simulation loop, for each timestep, we use our ODE solver and system equations to evaluate the new state vectors, check if there has been a surface collision, and save the position vector to an output array.
Finally, we can visualize our two output graphs
As we can see, we have a nice circular 200km orbit. we can use different initial conditions to look at a few other orbits. Here is a retrograde orbit with a starting velocity of -10.8 km/s and a total simulation time increased to 1400 minutes (around 23 hours).
The escape velocity of the earth is 11.2 km/s so therefore if we use our original initial conditions of a 200km orbit but with a velocity of 11.3 km/s we should see a trajectory that never comes back to the Earth. Here is the plot after 100 minutes.
And after 24 hours
We can see that our orbit is asymptotically approaching a line and will never return. Next, we can see what happens if our orbital velocity is too low, around 6.8 km/s
We can see that we have impacted the surface. Another thing that we can do is to add more dimensions to our initial velocity vectors. With an initial x component of 3.9 km/s and the initial y component of 7.8 km/s we get the following orbit.
if we extend the simulation time longer, we can see that eventually, we will impact the surface. This high-arcing suborbital trajectory is similar to the flight path of a ballistic missile.
Conclusion
If you made it all the way to the end good job! This was meant to be a nice simple introduction to physical simulations with a fun application. The code used for this article can be found on GitHub at the following link GitHub; feel free to play around with it on your own! In future articles, we will explore post-processing your simulation data as well as how to write other simulations and maybe even use them to make informational 3D movies.
Sources
03 – Writing a Billiard Simulation in 10 Minutes, 2021. https://www.youtube.com/watch?v=ThhdlMbGT5g.
David. “Vehicle Dynamics: Basic Concepts,” August 25, 2020. https://thef1clan.com/2020/08/25/vehicle-dynamics-basic-concepts/.
Inside 6 DoF Simulations, 2022. https://www.youtube.com/watch?v=30EIDIXvb2Y.
“Numerical Approach to Studying Vehicle Dynamics with a Half-Car Suspension Model,” January 8, 2018. https://nrsyed.com/2018/01/07/numerical-approach-to-studying-vehicle-dynamics-with-a-half-car-suspension-model/.
Rocket Fundamentals (Ideal Rocket Equation Derivation, Specific Impulse) | Rocket Trajectories 1, 2021. https://www.youtube.com/watch?v=bPXjkFAnQio.
“Runge-Kutta Integrator Overview: All Purpose Numerical Integration of Differential Equations – YouTube.” Accessed December 19, 2022. https://www.youtube.com/watch?v=HOWJp8NV5xU.
Sounding Rocket Trajectories | Rocket Trajectories 2, 2021. https://www.youtube.com/watch?v=-4-UfoFgU0A.
Witkin, Andrew. “Physically Based Modeling: Principles and Practice,” n.d.