Numerical model of the Solar System

A Solar system numerical model

A law of motion for entities in a Solar system numerical model

The Solar system is a planetary system that may be modeled numerically using both Newton's Law of Gravitation and Newton's Second Law of Motion.

Newton's Universal Law of Gravitation states that if a test mass m resides at a coordinate vector ${\displaystyle {\vec {r}}}$ , then the gravitational force acting on the test mass, due to another mass that is present is:

${\displaystyle {\vec {F}}=G{\frac {mM_{i}}{|{\vec {r}}_{i}-{\vec {r}}|^{3}}}({\vec {r}}_{i}-{\vec {r}})}$ 

where:

${\displaystyle {\vec {F}}}$  is the gravitational force acting on mass m.

${\displaystyle G}$  is the gravitation constant.

${\displaystyle M_{i}}$  is the mass of an entity i. For purposes of modeling a planetary system, it is assumed to be a mass that is large enough to make a difference in the motion of test mass m. It could be a star, a planet, a moon, a large asteroid, or maybe even a large comet. A small rock or a spacecraft would not be massive enough.

${\displaystyle {\vec {r}}_{i}}$  is the coordinate vector of entity i.

In a planetary system, there will be many entities i. Therefore, the total gravitational force acting on test mass m will be the sum of the forces for each i, over all n forces that are present:

${\displaystyle {\vec {F}}=\sum _{i}^{n}G{\frac {mM_{i}}{|{\vec {r}}_{i}-{\vec {r}}|^{3}}}({\vec {r}}_{i}-{\vec {r}})}$ 

Here we have changed the definition of ${\displaystyle {\vec {F}}}$  to mean the total gravitational force due to all entities i.

${\displaystyle {\vec {F}}}$  is the total gravitational force acting on mass m.

Newton's Second Law of Motion states that if a net force of ${\displaystyle {\vec {F}}}$  acts on a test mass m, that the test mass will undergo an acceleration ${\displaystyle {\vec {a}}}$  such that all three quantities are related by:

${\displaystyle {\vec {F}}=m{\vec {a}}}$ 

These two laws have the common left-hand-side ${\displaystyle {\vec {F}}}$ . They may therefore be combined. When this is done and when the resulting equation is divided through by the common factor m, an equation of motion for test mass m results:

${\displaystyle {\vec {a}}=\sum _{i}^{n}G{\frac {M_{i}}{|{\vec {r}}_{i}-{\vec {r}}|^{3}}}({\vec {r}}_{i}-{\vec {r}})}$ 

At this point the concept of a test mass m is dispensed with. In order to model a solar system, each of the entities j of mass ${\displaystyle M_{j}}$ , in turn, take the place of the test mass. The coordinate vector ${\displaystyle {\vec {r}}}$  of the test mass therefore becomes the coordinate vector ${\displaystyle {\vec {r}}_{j}}$  of entity j, and the summation is done over all other entities i except for the case where i equals j itself.

${\displaystyle {\vec {a}}_{j}=\sum _{i\neq j}^{n}G{\frac {M_{i}}{|{\vec {r}}_{i}-{\vec {r}}_{j}|^{3}}}({\vec {r}}_{i}-{\vec {r}}_{j})}$ 

where obviously both i and j run from 1 to n.

The above equation is written in vector form. In order to have a set of formulas that is useful, the vector equation must be split into its three vector component equations.

${\displaystyle (a_{j})_{x}=\sum _{i\neq j}^{n}G{\frac {M_{i}}{d^{3}}}(x_{i}-x_{j})}$ 

${\displaystyle (a_{j})_{y}=\sum _{i\neq j}^{n}G{\frac {M_{i}}{d^{3}}}(y_{i}-y_{j})}$ 

${\displaystyle (a_{j})_{z}=\sum _{i\neq j}^{n}G{\frac {M_{i}}{d^{3}}}(z_{i}-z_{j})}$ 

where

${\displaystyle d=((x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}+(z_{i}-z_{j})^{2})^{1/2}}$ 

At this point we have an equation of motion for all entities j in the panetary system to be modeled.

Reduced units

The mks system of units consists of three basic units for the three basic mechanical scalars. Length is in meters [m], mass is in kilograms [kg], and time is in seconds [s].

This would be a terribly awkward system of units for modeling a planetary system. A far more convenient system of units would be: length in astronomical units [ AU ], mass in Terran masses [mTerran], and time in days [d].

An astronomical unit [AU] is the mean distance from the center of the Earth to the center of the Sun, over the course of any given recent year. In this system of units, the Earth could be given an initial coordinate of, say, (x = 1.0, y = 0.0). A planet that is toward the outer edge of the Solar system, such as Pluto, could be given an initial coordinate of, say, (x = 40.0, y = 0.0). The magnitude of the distances that occur in a Solar system model is quite manageable in this system of units.

A Terran mass [mTerran] is simply the mass of the Earth. In this system of units, the Earth would have a mass of (m = 1.0). A large planet, (the gas giant Jupiter being an example of such), could have a mass of, say (m = 300.0). A star, (our Sun being an example of such), could have a mass of, say, (m = 300,000.0) The magnitude of masses that occur in a solar system model is quite intuitive in this system of units. It is fairly easy to picture that another planet could have, say 10 times, or 1000 times the mass of the Earth, and that a star could have, say 1 million times the mass of the Earth.

A mean solar day of Earth [d] is a familiar unit that is easy to picture. It takes our Earth a year made up of 364.24 means solar days to orbit 360 degrees around our Sun. If a planet at the outer edge of the solar system took, say, 200.0 years to orbit around the Sun, then this would be an interval of time that could be easily pictured, and easily converted to days. The magnitude of the time intervals that occur in a Solar system model is both manageable and intuitive in this system of units.

For reference, the conversion between the MKS system, and our AU-mTerran-d system, is given below:

 1.0 AU      = 149.5985 x 10^9 m
 1.0 mTerran = 5.977 x 10^24 kg
 1.0 d       = 60. * 60. * 24. s = 86,400.0 s

The Universal Constant of Gravitation, G is, in MKS units:

 G = 6.673 x 10^-11 m^3 kg^-1 s^-2

If we are going to go ahead and use our AU-mTerran-d system of units, then we must convert this constant to that system of units and use the converted constant in our numerical model of the Solar system. Let us call this converted constant 'g'.

 g = 6.673 x 10^-11 (m^3 kg^-1 s^-2) (149.5985 x 10^9 m/AU)^-3 (5.977 x 10^24 kg/mTerran)^+1 (86,400.0 s/d)^+2

   = (6.673 x 10^-11) (149.5985 x 10^9)^-3 (5.977 x 10^24)^+1 (86,400.0)^-2 AU^3 mTerran^-1 d^-2

   = 1.595 861 x 10^-29 AU^3 mTerran^-1 d^-2

It is this converted constant g that will be used in the code. The converted constant g may appear to be a less convenient number than the original constant G. However, the more convenient value for G appears in one of the input files, and the awkward conversion to g is done automatically in the code.

The initial velocities of the Solar system entities are given in units of [km/s]. However, users who have become used to our AU-mTerran-d system may wish to assign velocities in units of [AU/d]. The conversion factors are:

 1.0 AU / d = 1731.464 120 km / s
 1.0 km / s = 0.000 577 459 AU / d

The code expects the initial velocities to be in units of [km/s]. The code was written this way for the reason that users may wish to generalize it so that the trajectories of spacecraft through the Solar system may be calculated. Spacecraft using conventional chemical rocket fuels function by applying an engine burn in order to attain a change in velocity (a delta-v) in any given direction. The distance over which the spacecraft travels while a given engine burn is applied, is negligible compared to the distance scale of the Solar system through which it will travel. Therefore a given delta-v may be considered to occur instantaneously. A delta-v is typically in units of [km/s]. Therefore while distances are often given in [AU], and time intervals in [d], velocities are often given in [km/s].

Now that we have both a set of equations of motion and a convenient system of units, we next need an algorithm for computing the coordinates of the Solar system entities, as time passes.

A time-stepping algorithm

The next step is to come up with an algorithm in which the coordinates ${\displaystyle {\vec {r}}_{j}}$  may be calculated as time passes.

We start with the definitions for velocity and acceleration:

${\displaystyle {\vec {a}}={\frac {d{\vec {v}}}{dt}}}$ 

${\displaystyle {\vec {v}}={\frac {d{\vec {r}}}{dt}}}$ 

If instead of an infinitesimal time step ${\displaystyle dt}$  we consider a finite time step ${\displaystyle \delta t}$ , then we obtain definitions for the average velocity and acceleration over the time interval ${\displaystyle \delta t}$ 

${\displaystyle {\vec {a}}_{ave}={\frac {\Delta {\vec {v}}}{\Delta t}}}$ 

${\displaystyle {\vec {v}}_{ave}={\frac {\Delta {\vec {r}}}{\Delta t}}}$ 

or

${\displaystyle \Delta {\vec {v}}={\vec {a}}_{ave}\Delta t}$ 

${\displaystyle \Delta {\vec {x}}={\vec {v}}_{ave}\Delta t}$ 

Splitting these vector equations into the vector component equations:

${\displaystyle (a_{ave})_{x}={\frac {\Delta v_{x}}{\Delta t}}}$ 

${\displaystyle (a_{ave})_{y}={\frac {\Delta v_{y}}{\Delta t}}}$ 

${\displaystyle (a_{ave})_{z}={\frac {\Delta v_{z}}{\Delta t}}}$ 

${\displaystyle (v_{ave})_{x}={\frac {\Delta x}{\Delta t}}}$ 

${\displaystyle (v_{ave})_{y}={\frac {\Delta y}{\Delta t}}}$ 

${\displaystyle (v_{ave})_{z}={\frac {\Delta z}{\Delta t}}}$ 

or:

${\displaystyle \Delta v_{x}=(a_{ave})_{x}\Delta t}$ 

${\displaystyle \Delta v_{y}=(a_{ave})_{y}\Delta t}$ 

${\displaystyle \Delta v_{z}=(a_{ave})_{z}\Delta t}$ 

${\displaystyle \Delta x=(v_{ave})_{x}\Delta t}$ 

${\displaystyle \Delta y=(v_{ave})_{y}\Delta t}$ 

${\displaystyle \Delta z=(v_{ave})_{z}\Delta t}$ 

Next we need an algorithm for stepping the coordinates (x, y, z)_j through time. Let us consider just coordinate x for just one entity, and let us consider the following algorithm:

 input t_f, (delta)t         (We need input values for the time step size (delta)t
                              and the final time t_f)
 t = 0.                      (We start the time at zero)
 x = x_0                     (We set the coordinate x to its initial value)
 v = v_0                     (We set the velocity v to its initial value)
 a = f(x)                    (The acceleration is a function of x and we calculate it here) (3.)
 output t, x
 loop until (t = t_f)
 |   t = t + (delta)t        (We advance the time)
 |   x = x + v * (delta)t    (We advance the coordinate) (1.)
 |   v = v + a * (delta)t    (We advance the velocity)   (2.)
 |   a = f(x)                (We calculate the new acceleration as a function of x) (3.)
 |   output t, x

This algorithm is called the Euler Algorithm. It is the most basic time stepping algorithm that can be used to solve a point particle classical mechanics problem. The Euler Algorithm, however, has problems. The formulas (1.) and (2.) are based on the definitions given earlier of average velocity and average acceleration over a finite time interval, and yet it uses velocities and accelerations that occur at the beginning of the time interval.

A better algorithm is to use velocities and accelerations that occur in the middle of the time intervals. This is quite simply accomplished. We start by advancing the velocity by only a half step rather than a full step. From that point on, we advance the coordinate and the velocity by a full step, except that in this algorithm, the velocity that is used to advance the coordinate, occurs in the middle of the coordinate time step, and the acceleration that is used to advance the velocity, occurs in the middle of the velocity time step.

 input t_f, (delta)t         (We need input values for the time step size (delta)t
                              and the final time t_f)
 t = 0.                      (We start the time at zero)
 x = x_0                     (We set the coordinate x to its initial value)
 a = f(x)                    (The acceleration is a function of x and we calculate it here) (3.)
 v = v_0                     (We set the velocity v to its initial value)
 output t, x
 
 v = v + a * (delta)t/2.     (We advance the velocity by half a time step)
 
 loop until (t = t_f)
 |   t = t + (delta)t        (We advance the time)
 |   x = x + v * (delta)t    (We advance the coordinate by a full time step) (1.)
 |   a = f(x)                (We calculate the new acceleration as a function of x) (3.)
 |   v = v + a * (delta)t    (We advance the velocity by a full time step)   (2.)
 |   output t, x

This algorithm is called the Feynman Algorithm, or the Verlet Algorithm.

We now have all the tools we need to construct a numerical model of the solar system. We have an equation of motion for all the solar system entities j, which provides us with the acceleration of all the solar system entities j at any given instant. We also have a time stepping algorithm which uses this acceleration information for a given entity to step its coordinates and velocity through time.

A working code

The following code constitutes such a numerical model. The code is written in fortran 77. It reads an input file containing start-up information for a solar system, and writes an output file that contains the coordinates of each solar system entity, at each time step in the simulation, from an initial time of zero, up until the final time.

      program solarV2f
      implicit none
c
      call readConsts()
      call readSystem()
      call moveSystem()
c
      stop
      end
c
      subroutine readConsts()
      implicit none
c
      double precision rAu, mTerran, g, deltaT, tF
      double precision tDay
      common /fundConsts/ rAu, mTerran, g, deltaT, tF, tDay
c
      integer nin, nout, nlog
      parameter(nin=1, nout=2, nlog=3)
c
      open(nin, file='const.dat', status='old')
      read(nin, *) rAu
      read(nin, *) mTerran
      read(nin, *) g
      read(nin, *) deltaT
      read(nin, *) tF
      close(nin)
c
      tDay = 60. * 60. * 24.
      g = g * mTerran * tDay * tDay / (rAu * rAu * rAu)
c
      return
      end
c
      subroutine readSystem()
      implicit none
c
      double precision PI
      parameter(PI=3.141592)
c
      double precision rAu, mTerran, g, deltaT, tF
      double precision tDay
      common /fundConsts/ rAu, mTerran, g, deltaT, tF, tDay
c
      integer nMax, n
      parameter(nMax = 20)
      character*20 name(nMax)
      double precision mass(nMax)
      double precision rx(nMax), ry(nMax), rz(nMax)
      double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)
      double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)
      common /size/ n
      common /system/ name, mass, rx, ry, rz,
     &  vOldx, vOldy, vOldz, vAvex, vAvey, vAvez
c
      double precision r, theta, phi
      double precision vr, vtheta, vphi
      double precision vx(nMax), vy(nMax), vz(nMax)
      integer i
c
      integer nin, nout, nlog
      parameter(nin=1, nout=2, nlog=3)
c
      open(nin, file='solar.dat', status='old')
      read(nin, *) n
      do i = 1, n
        read(nin, 10) name(i)
        read(nin, *) mass(i)
        read(nin, *) r, theta, phi
        theta = theta * (PI/180.)
        phi = phi * (PI/180.)
        rx(i) = r * sin(theta) * cos(phi)
        ry(i) = r * sin(theta) * sin(phi)
        rz(i) = r * cos(theta)
c
        read(nin, *) vr, vtheta, vphi
        vx(i) = vr * sin(theta) * cos(phi)
     &        + vtheta * cos(theta) * cos(phi)
     &        - vphi * sin(phi)
        vy(i) = vr * sin(theta) * sin(phi)
     &        + vtheta * cos(theta) * sin(phi)
     &        + vphi * cos(phi)
        vz(i) = vr * cos(theta)
     &        - vtheta * sin(theta)
c
        vx(i) = (tDay * 1000. / rAu) * vx(i)
        vy(i) = (tDay * 1000. / rAu) * vy(i)
        vz(i) = (tDay * 1000. / rAu) * vz(i)
c
        vOldx(i) = vx(i)
        vOldy(i) = vy(i)
        vOldz(i) = vz(i)
c
        vAvex(i) = vx(i)
        vAvey(i) = vy(i)
        vAvez(i) = vz(i)
c
      enddo
      close(nin)
      return
   10 format(a20)
      end
c
      subroutine moveSystem
      implicit none
c
      double precision rAu, mTerran, g, deltaT, tF
      double precision tDay
      common /fundConsts/ rAu, mTerran, g, deltaT, tF, tDay
c
      integer nMax, n
      parameter(nMax = 20)
      character*20 name(nMax)
      double precision mass(nMax)
      double precision rx(nMax), ry(nMax), rz(nMax)
      double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)
      double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)
      common /size/ n
      common /system/ name, mass, rx, ry, rz,
     &  vOldx, vOldy, vOldz, vAvex, vAvey, vAvez
c
      integer i
      double precision t
      double precision vNewx, vNewy, vNewz
      double precision ax, ay, az
c
      integer nin, nout, nlog
      parameter(nin=1, nout=2, nlog=3)
c
      open(nout, file='coords.dat', status='unknown')
c
c     Start at time zero
c
      t = 0.
      call writeSystem(t)
c
c     Perform initial half-step
c
      do i = 1, n
        call accelMass(i, g, ax, ay, az)
        vOldx(i) = vOldx(i) + (deltaT * ax) / 2.
        vOldy(i) = vOldy(i) + (deltaT * ay) / 2.
        vOldz(i) = vOldz(i) + (deltaT * az) / 2.
      enddo
c
c     Perform time steps
c
      do while (t < tF)
        t = t + deltaT
        do i = 1, n
          rx(i) = rx(i) + (deltaT * vOldx(i))
          ry(i) = ry(i) + (deltaT * vOldy(i))
          rz(i) = rz(i) + (deltaT * vOldz(i))
        enddo
        do i = 1, n
          call accelMass(i, g, ax, ay, az)
          vNewx = vOldx(i) + (deltaT * ax)
          vNewy = vOldy(i) + (deltaT * ay)
          vNewz = vOldz(i) + (deltaT * az)
          vAvex(i) = (vOldx(i) + vNewx) / 2.
          vAvey(i) = (vOldy(i) + vNewy) / 2.
          vAvez(i) = (vOldz(i) + vNewz) / 2.
          vOldx(i) = vNewx
          vOldy(i) = vNewy
          vOldz(i) = vNewz
        enddo
        call writeSystem(t)
      enddo
c
c     Calculations are done
c
      close(nout)
      return
      end
c
      subroutine accelMass(i, g, ax, ay, az)
      implicit none
c
      integer i
      double precision g
      double precision ax, ay, az
c
      integer nMax, n
      parameter(nMax = 20)
      character*20 name(nMax)
      double precision mass(nMax)
      double precision rx(nMax), ry(nMax), rz(nMax)
      double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)
      double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)
      common /size/ n
      common /system/ name, mass, rx, ry, rz,
     &  vOldx, vOldy, vOldz, vAvex, vAvey, vAvez
c
      integer j
      double precision rsx, rsy, rsz
      double precision rox, roy, roz
      double precision rosx, rosy, rosz
      double precision ms
      double precision r2, r, r3
c
c     Initialize the acceleration
c
      ax = 0.
      ay = 0.
      az = 0.
c
c     Set the coordinates of the observer mass
c
      rox = rx(i)
      roy = ry(i)
      roz = rz(i)
c
c     Loop over all other masses
c
      do j = 1, n
        if (j .ne. i) then
          ms = mass(j)
c
          rsx = rx(j)
          rsy = ry(j)
          rsz = rz(j)
c
          rosx = rsx - rox
          rosy = rsy - roy
          rosz = rsz - roz
c
          r2 = rosx*rosx + rosy*rosy + rosz*rosz
          r = dsqrt(r2)
          r3 = r * r * r
c
          ax = ax + g * ms * rosx / r3
          ay = ay + g * ms * rosy / r3
          az = az + g * ms * rosz / r3
        endif
      enddo
c
      return
      end
c
      subroutine writeSystem(t)
      implicit none
c
      double precision t
c
      integer nMax, n
      parameter(nMax = 20)
      character*20 name(nMax)
      double precision mass(nMax)
      double precision rx(nMax), ry(nMax), rz(nMax)
      double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)
      double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)
      common /size/ n
      common /system/ name, mass, rx, ry, rz,
     &  vOldx, vOldy, vOldz, vAvex, vAvey, vAvez
c
      integer i
c
      integer nin, nout, nlog
      parameter(nin=1, nout=2, nlog=3)
c
c     Write the time
c
      write(nout,200) t
c
c     Write the number of masses in this time frame
c
      write(nout,100) n
c
c     Loop over all the masses and write the coordinates of each one
c
      do i = 1, n
        write(nout,300) rx(i), ry(i), rz(i)
      enddo
c
      return
  100 format(I4)
  200 format(F12.6)
  300 format(F12.6,4X,F12.6,4X,F12.6)
      end
c

The following file contains the physical constants that the numerical model requires, as well as the time information.

149.5985e09
5.977e24
6.673e-11
1.
1000.

The file is in the form of:

 [1.0 AU, in m]
 [1.0 mTerran, in kg]
 [G, in m^3 kg^-1 s^-2]
 [Time step, in d]
 [Total time interval of simulation, in d]

The next file contains contains the information for a simplified Solar system.

3
Sun
332943
0.031254    90.         180.
0.          0.          0.079
Jupiter
2000
5.203       88.694      0.
0.          0.          13.07
Earth
1.
1.          90.         0.
0.          0.          29.78

The file is in the form of:

 [Number of Solar system entities]
 [Name of Solar system entity 1]
 [Mass in mTerran]
 [Initial Radius from center, in AU]              [Initial Theta, in degrees]                    [Initial Phi, in degrees]
 [Initial Radius-direction velocity, in km/s]     [Initial Theta-direction velocity, in km/s]    [Initial Phi-direction velocity, in km/s]

The time info in the first file, and the solar system info in the second file, may be modified.

This article is intended to explain the physical theory that is required to numerically simulate Solar system dynamics using Newton's Law of Gravitation, and Newton's Second Law of Motion, and to provide a simple code that serves as an example of how the theory may be coded. No attempt has been made here to graphically display or animate the motion of the Solar system entities.

Users may wish to develop more elaborate codes that will do this. They may wish to use the code provided here as the 'physics engine' in such a code.

Readers may note that the code performs some complicated coordinate system conversion calculations, particularly in the calculation of the initial velocity components. These calculations are described in a separate article, on a transformation from spherical coordinates to rectangular coordinates.

See also

• astronomical algorithm