Satellites are very important part of modern civilization, they are used for many purposes and each of us is using or at least has used some of their services. I won’t be however talking about the purposes of the satellites today. I would rather spend some time writing about how do they stay on the sky, what needs to be fulfilled and what is their trajectory, which is not a circular orbit as it may sound. We will get there there soon. First of all, lets take a look on the derivation of the circular orbit. We will need the Newton’s law of gravitation and the description of the centrifugal force. We can say, that the circular orbit is a situation, where all forces are in balance eg: Eq1

Where κ is the universal gravitation constant, which is equal to 6.67 10-11m3kg-1s-2 and Mg is the mass of the earth, which is approximately 6.1024kg. One can easily obtain the formula needed to calculate the requested velocity based on the distance from the center of the earth: (Substituing r = 6378 000 + 100 000 | Radius of the Earth + Altitude of 100km | yields velocity of 7.86 km/s. The higher the distance of the 2 object,the slower the object moves. This is in fact one of the Kepler’s laws.)


For example, GPS Satellites are located at the MEO (Medium Orbit, which is approximately 22 000 km above the surface. Their speed is “only 3.8 km/s”). This circular orbit speed is also the first cosmic speed. Giving the object lower speed (at the same altitude) would result in the eventual crash of the object on the Earth surface. The question then is,if we increase the speed,what can happen? There are 2 possibilities. The first and in case of most the satellites is that the orbit’s shape changes from circular into the elliptical one. Satellites with speed higher than the first cosmic speed moves on eliptical orbits. Even this has however some limitations. The second cosmic speed to be correct. Its derivation can be found elsewhere and is not extra hard,the formula says:


As you can see,the formula is very similar and in our case of altitude of 100km above the earth surface is approximately 11km/s. This is the speed that needs to be given to an object to escape from the gravitational field of the earth. Objects with speeds equal to or higher than the second cosmic speed no longer have “orbital trajectories” instead, they just fly away on a parabolic trajectory. Of course there is another cosmic speed and no surprise its the third cosmic speed. this is the speed that need to be reached in order to escape from the gravitational field of the sun. I will however now concentrate on a different and general topic,which includes the Newton’s law of gravitation along with Newton’s second law of motion:


We have just subsititued the general force with the gravitational force. This is still nothing worth noticing except for the minus sign, that implies the fact that the moves towards the source of the force (meaning just direction in other words), so we take a look at a more generalized form of the equation by taking the definition of the acceleration a as the second derivate of the location and further considering a 2-D space:


We can see, that the situation is much more interesting and not as easy to understand. First of all, the distance is defined by both spacial variables x and y. Then at second in (3) we see the x – acceleration defined as the second derivative of the x – position. This is identical to situation in (4) with the y spacial variable. Furthermore, we can see some strange sin() and cos() functions appearing in the Equations. These functions however only represents the projection of the total force acting onto the object into both spatial dimensions x and y. Users with basic knowledge of solving differential equations should know, that in order to solve this system of 2 second order differential equations, we would need 4 initial conditions, that represents the initial position of the object in 2D space along with its velocity in both dimensions. Although these equations can be probably solved analytically, I will concentrate on the numerical solutions, which includes discretization in time and the 2 spacial dimensions x and y.

The hardest thing I had to overcome was the fact, that these are the second order differential equations,which I have never been solving before. Therefore I will outline the basic concept of solving them,which is that you really cant solve them at once. Instead, the best option is to divide each of the second order equations into 2 first order differential equations thus increasing the total number of ODEs (Ordinary differential equations) to 4. These are however easy to solve and their implementations can be found in any of the matlab functions provided (EulerCompute / RK2Compute / RK4Compute / R6Compute). All of these functions are some numerical methods I have used to solve this problem. I suggest to start with the Euler’s method implementation, because it is the easiest to understand. RK stands for Runge-Kutta and the number represents the order of the method. Higher orders are more accurate at the expense of more operations.


I call the s1 and s2 the “help” variables,although more accurate would be the interpretation of the velocity – please continue for details. To illustrate Euler’s method, we take a look at (5) for example. We would like to know x, not its derivative (Which is the velocity … you know dx/dt = v ) … In order to obtain x, we just integrate the whole equation thus obtaining x. Its as simple as it sounds. 


I have created a script demonstrating the movement of an object with the gravitational pull of the “Earth” – It doesn’t have to be necessary Earth, but because I have used the Earth’s mass for calculations, I call it the Earth simulation. This parameter can be adjusted in the script as well as the initial conditions and reference position of the Earth in 2D. The condition to stop the simulation is defined as crash to the Earth’s surface. I always do use the base units to make thing easier to understand,which on the other hand results in long numbers … One of the most important parameters is the time step however,which by default equals to 1 second. Increasing the time step speeds-up the simulation while greatly increasing the numerical errors. These can be suppressed by higher-order methods and or adaptive methods (Adaptivness of the algorithms is not implemented ) , but I would not recommend to increase it much. Just for demonstration purposes.

Matlab 2D demo: (E)MotionEQ
Matlab 3D demo: (E)MotionEQ

3D Version added 4.7.2016
It is implemented only as the 4th order Runge-Kutta method.
Due to 3D Graphics, the demo is slower than the 2D version.
Several improvements exist in the 3D. To speedup the computation,
angular corrections can be skipped via a variable in the main script.
The simulation core is almost identical to the 2D problem,
additional info can be found within the demo or in RK4Compute function.