# Putting it all together - a sixth-order integrator

#### MontBlanc2012

##### Member
This note is going to be about building a high-fidelity integration engine for Orbiter that takes into account the gravitational influence of other bodies (tidal forces) and the 'non-spherical gravity sources' that Orbiter attributes to each body. In a subsequent post, I will show how this integrator can be used to solve the orbital rendezvous problem.

To build this integration engine, I will make use of ideas presented in three recent posts in Orbiter’s Maths & Physics section. These three posts are:

* "The ECI reference frame and tidal forces"

* "J4 - J5 perturbations"

* "From BCI to BCBF and back again"

The first of these posts shows how to write the equations of motion of a satellite in the Earth-centred Inertial (ECI) reference frame. This is a natural reference frame to use when considering the motion of satellites within its sphere of influence. It is a non-rotating reference frame centred on the Earth's. The ECI reference frame (or more generally, the Body-centred Inertial (BCI) reference frame) is non-inertial. Because of this we need to add tidal forces attributable to other bodies to the equations of motion. For the Earth, it is common to include tidal forces attributable to the Moon and the Sun.

The second post shows how to calculate the force contributions from the 'non-spherical gravity sources'. Typically, these are the assorted spherical harmonics whose coefficients are given as [MATH]J_2[/MATH], [MATH]J_3[/MATH] and so on. (It is possible to build more complicated gravitational models for the Earth as discussed in that post but for Orbiter we need only consider the axially symmetrically contributions, i.e., the [MATH]J_n[/MATH] terms.) This gravitational model for the Earth can also be generalised to any gravitating body. These non-spherical gravity contributions are most conveniently expressed in the Earth-centred Earth-fixed (ECEF) - or more generally the Body-centred Body-fixed (BCBF) reference frame - which is a reference frame that co-rotates with the Earth (or, more generally, some other gravitating body).

This brings us on, then, to the third post which discusses the algorithms needed to convert from the ECI (or BCI) reference frame in which the equations of motion of a satellite are most conveniently expressed; and the ECEF (or BCBF) reference frame in which the contributions to gravitational forces from the non-spherical gravity sources are calculated.

These force calculations and coordinate transformations are then 'wrapped up' into a 6th order numerical integrator which, when run, produces a high-fidelity integration of the trajectory of a satellite moving in proximity to the Earth (or, more generally, some other gravitating body). Having built this integrator, I will then use this integrator to carry out some precision trajectory design work.

The integrator - a detailed description

To keep things as straight forward as possible, I'm simply going to present the Python code for the integrator. Then, I'm going to go through the code line-by-line to describe what the code is doing. So, without further ado, here is the integrator:

Code:
import numpy as np

# the universal gravitatinal constant (as per Orbiter documentation)
G = 6.67259e-11

mu_dict = { "Earth"  : G * 5.973698968e+24,
"Moon"   : G * 7.347664e+22,
"Sun"    : G * 1.9889194444e+30 }

R_dict  = { "Earth"  : 6.37101e6,
"Moon"   : 1.738e6,
"Sun"    : 6.96e8 }

J_dict  = { "Earth"  : [0.0, 0.0, 1082.6269e-6, -2.51e-6, -1.60e-6, -0.15e-6],
"Moon"   : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
"Sun"    : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] }

# this routine taks a vector 'Q' in the ECEF reference frame and returns the gravitational acceleration in the BCBF reference frame
def g(Q, body):

nmax  = 5
mu    = mu_dict[body]
R     = R_dict[body]
J     = J_dict[body]

x     = Q[0]
y     = Q[1]
z     = Q[2]

r     = np.sqrt(x*x + y*y + z*z)
xi    = z / r
phi   = R / r
k     = r * r; k *= k; k = mu / k

P     = np.zeros(nmax + 1)
DP    = np.zeros(nmax + 1)
f     = np.zeros(3)

P[0]  = 1.0; DP[0] = 0.0
P[1]  = xi; DP[1] = 1.0

for n in range(2, nmax + 1):
P[n]  = (2*n - 1) * xi * P[n-1] + (1-n) * P[n-2]; P[n] /= n
DP[n] = (2*n - 1) * P[n-1] + (2*n - 1) * xi * DP[n-1] + (1-n) * DP[n-2]; DP[n] /= n

q = phi
temp = mu / r / r / r;
f[0] = - temp * x
f[1] = - temp * y
f[2] = - temp * z
for n in range(2, nmax + 1):
q *= phi
f[0] += k * J[n] * q * ((n+1) * P[n] * x * r + DP[n] * x * z)
f[1] += k * J[n] * q * ((n+1) * P[n] * y * r + DP[n] * y * z)
f[2] += k * J[n] * q * ((n+1) * P[n] * z * r - DP[n] * (x * x + y * y))

return f

def BCBF2BCI( X, mjd ):

# Earth parameters
pLan = 0.0            # Precession LAN
pObl = 0.0            # Precession obliquity
Tp   = -9413040.4     # Precession period
LAN  = 0.0            # LAN
mjdL = 51544.5        # LAN MJD
Ts   = 86164.10132    # Sidereal rotation period
eps  = 0.4090928023   # Obliquity
psi0 = 4.88948754     # Sidereal rotation offset

deld = mjd - mjdL
tau  = LAN + 2 * np.pi * deld / Tp
psi  = 2 * np.pi * 86400.0 / Ts * deld - 2 * np.pi * np.cos(eps) / Tp * deld + psi0

x    = np.empty(6)
y    = np.empty(6)
z    = np.empty(6)

x[5] = X[0]
y[5] = X[1]
z[5] = X[2]

c    = np.cos(psi)
s    = np.sin(psi)
x[4] = +x[5] * c - y[5] * s
y[4] = +x[5] * s + y[5] * c
z[4] = +z[5]

c    = np.cos(eps)
s    = np.sin(eps)
x[3] = +x[4]
y[3] = +y[4] * c + z[5] * s
z[3] = -y[4] * s + z[5] * c

c    = np.cos(tau)
s    = np.sin(tau)
x[2] = +x[3] * c - y[3] * s
y[2] = +x[3] * s + y[3] * c
z[2] = +z[3]

c    = np.cos(pObl)
s    = np.sin(pObl)
x[1] = +x[2]
y[1] = +y[2] * c + z[2] * s
z[1] = -y[2] * s + z[2] * c

c    = np.cos(pLan)
s    = np.sin(pLan)
x[0] = +x[1] * c - y[1] * s
y[0] = +x[1] * s + y[1] * c
z[0] = +z[1]

return np.array([x[0],y[0],z[0]])

def BCI2BCBF( X, mjd ):

# Earth parameters
pLan = 0.0            # Precession LAN
pObl = 0.0            # Precession obliquity
Tp   = -9413040.4     # Precession period
LAN  = 0.0            # LAN
mjdL = 51544.5        # LAN MJD
Ts   = 86164.10132    # Sidereal rotation period
eps  = 0.4090928023   # Obliquity
psi0 = 4.88948754     # Sidereal rotation offset

deld = mjd - mjdL
tau  = LAN + 2 * np.pi * deld / Tp
psi  = 2 * np.pi * 86400.0 / Ts * deld - 2 * np.pi * np.cos(eps) / Tp * deld + psi0

x    = np.empty(6)
y    = np.empty(6)
z    = np.empty(6)

x[0] = X[0]
y[0] = X[1]
z[0] = X[2]

c    = np.cos(pLan)
s    = np.sin(pLan)
x[1] = +x[0] * c + y[0] * s
y[1] = -x[0] * s + y[0] * c
z[1] = +z[0]

c    = np.cos(pObl)
s    = np.sin(pObl)
x[2] = +x[1]
y[2] = +y[1] * c - z[1] * s
z[2] = +y[1] * s + z[1] * c

c    = np.cos(tau)
s    = np.sin(tau)
x[3] = +x[2] * c + y[2] * s
y[3] = -x[2] * s + y[2] * c
z[3] = +z[2]

c    = np.cos(eps)
s    = np.sin(eps)
x[4] = +x[3]
y[4] = +y[3] * c - z[3] * s
z[4] = +y[3] * s + z[3] * c

c    = np.cos(psi)
s    = np.sin(psi)
x[5] = +x[4] * c + y[4] * s
y[5] = -x[4] * s + y[4] * c
z[5] = +z[4]

return np.array([x[5],y[5],z[5]])

def a_mon(Q_BCI, mjd):

# unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

# gravitational acceleration of the Moon due to the Earth
Q_mon_BCBF = BCI2BCBF(Q_mon_BCI, mjd)
A_mon_BCBF = g( Q_mon_BCBF, "Earth")
A_mon_BCI  = BCBF2BCI(A_mon_BCBF, mjd)

# correct for reduced mass of the Earth-Moon system
A_mon_BCI  = (1 + mu_dict["Moon"] / mu_dict["Earth"]) * A_mon_BCI

# add the tidal term due to the Sun's gravitational influence
A_mon_BCI += g( Q_sun_BCI, "Sun") - g( Q_sun_BCI - Q_mon_BCI, "Sun")

return A_mon_BCI

def a_sun(Q_BCI, mjd):

# unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

# gravitation acceleration of the Sun due to the Earth
Q_sun_BCBF = BCI2BCBF(Q_sun_BCI, mjd)
A_sun_BCBF = g( Q_sun_BCBF, "Earth")
A_sun_BCI  = BCBF2BCI(A_sun_BCBF, mjd)

# correct for reduced mass of the Earth-Sun system
A_sun_BCI  = (1 + mu_dict["Sun"] / mu_dict["Earth"]) * A_sun_BCI

# add the tidal term due to the Moon's gravitational influence
A_sun_BCI += g( Q_mon_BCI, "Moon") - g( Q_mon_BCI - Q_sun_BCI, "Moon")

return A_sun_BCI

def a_veh(Q_BCI, mjd):

# unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

# gravitation acceleration of the Vehicle due to the Earth
Q_veh_BCBF = BCI2BCBF(Q_veh_BCI, mjd)
A_veh_BCBF = g(Q_veh_BCBF, "Earth")
A_veh_BCI  = BCBF2BCI(A_veh_BCBF, mjd)

# add the tidal terms due to the gravitational influece of the Moon and the Sun
A_veh_BCI += g( Q_mon_BCI, "Moon") - g( Q_mon_BCI - Q_veh_BCI, "Moon")
A_veh_BCI += g( Q_sun_BCI, "Sun" ) - g( Q_sun_BCI - Q_veh_BCI, "Sun" )

return A_veh_BCI

def Integrator2(Q0, V0, mjd, dt):

# unpack the state vectors
[Q0_mon, Q0_sun, Q0_veh]  = Q0
[V0_mon, V0_sun, V0_veh]  = V0

# first 'kick' step
Va_mon  = V0_mon + 0.5 * dt * a_mon(Q0, mjd)
Va_sun  = V0_sun + 0.5 * dt * a_sun(Q0, mjd)
Va_veh  = V0_veh + 0.5 * dt * a_veh(Q0, mjd)

# drift step
Q1_mon  = Q0_mon + dt * Va_mon
Q1_sun  = Q0_sun + dt * Va_sun
Q1_veh  = Q0_veh + dt * Va_veh

# repack the position vector and update time
Q1  = [Q1_mon, Q1_sun, Q1_veh]
mjd = mjd + dt / 86400.0

# second 'kick' step
V1_mon  = Va_mon + 0.5 * dt * a_mon(Q1, mjd)
V1_sun  = Va_sun + 0.5 * dt * a_sun(Q1, mjd)
V1_veh  = Va_veh + 0.5 * dt * a_veh(Q1, mjd)

# repack the velocity vector
V1 = [V1_mon, V1_sun, V1_veh]

return Q1, V1, mjd

def Integrator4(Q0, V0, mjd, dt):

w0 = +1.35120719195965763
w1 = -1.70241438391931527

Qa, Va, mjd = Integrator2(Q0, V0, mjd, dt * w0)
Qb, Vb, mjd = Integrator2(Qa, Va, mjd, dt * w1)
Q1, V1, mjd = Integrator2(Qb, Vb, mjd, dt * w0)

return Q1, V1, mjd

def Integrator6(Q0, V0, mjd, dt):

w0 = +1.17467175808936338
w1 = -1.34934351617872677

Qa, Va, mjd = Integrator4(Q0, V0, mjd, dt * w0)
Qb, Vb, mjd = Integrator4(Qa, Va, mjd, dt * w1)
Q1, V1, mjd = Integrator4(Qb, Vb, mjd, dt * w0)

return Q1, V1, mjd

Even a cursory look through this code suggests that quite a lot is going on in this integrator. But it is all necessary if one wants to build a precision tool.

(I should also mention that this integrator has been tailored to integrate the equations of motion for one satellite in the vicinity of Earth in Orbiter 2016. However, this code can be generalised quite easily to other gravitating bodies and other ephemeris systems without much effort.)

Let's start going though this Python code, so that we can begin to see what's going on and why. To begin with, we have:

Code:
import numpy as np

This code snippet imports a Python package called 'Numpy' which contains a number of useful tools for manipulating vectors and matrices that aren't available in the standard Python installation.

Next, we have:

Code:
# the universal gravitatinal constant (as per Orbiter documentation)
G = 6.67259e-11

mu_dict = { "Earth"  : G * 5.973698968e+24,
"Moon"   : G * 7.347664e+22,
"Sun"    : G * 1.9889194444e+30 }

R_dict  = { "Earth"  : 6.37101e6,
"Moon"   : 1.738e6,
"Sun"    : 6.96e8 }

J_dict  = { "Earth"  : [0.0, 0.0, 1082.6269e-6, -2.51e-6, -1.60e-6, -0.15e-6],
"Moon"   : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
"Sun"    : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] }

This snippet sets some global constants that are relevant for integrating trajectories in the vicinity of the Earth. The first variable that is set is 'G', the universal gravitation constant. The value that it is set to is the value given to it in Orbiter 2016, rather than the current generally accepted best estimate.

Next, we define a number of Python 'dictionaries' (or look-up tables):

* 'mu' the standard gravitational parameters for the Earth, the Moon and the Sun;

* 'R' and 'J' the reference radius and non-spherical gravity coefficients used by Orbiter for each body.

All of these values have been extracted from Orbiter's planetary configuration files. Note that although we are integrating the equations of motion in the Earth-centred Inertial reference frame, we include parameter values for the Moon and for the Sun because we want to include the tidal forces from those bodies in the equations of motion.

The code then goes on to define three functions:

Code:
def g(Q, body):

def BCBF2BCI( X, mjd ):

def BCI2BCBF( X, mjd ):

The first of these is basically the same routine as set out (and explained) in the earlier post "J4 - J5 perturbations". It calculates the gravitational force of a body at an x-y-z point Q in the BCBF reference frame. The second and third are routines that were presented in the post "From BCI to BCBF and back again". As the names suggest the BCBF2BCI routine transforms a vector 'X' in the BCBF reference frame to a vector in the BCI reference frame on a given date ('mjd'). The inverse routine to take the vector back to the BCI reference frame is BCI2BCBF. At the moment, these routines are hard-wired to calculate the transformations for the ECI and ECEF reference frames but it is easy to generalise these to other gravitating bodies other than the Earth.

Next we have three routines that define the gravitational forces acting on the Moon, the Sun and the satellite in the ECI reference frames:

Code:
def a_mon(Q_BCI, mjd):

def a_sun(Q_BCI, mjd):

def a_veh(Q_BCI, mjd):

Now, I am calculating the gravitational forces acting on the Moon and the Sun in this algorithm because I want to calculate the tidal forces acting on a satellite attributable to these two bodes. To do this, I need to know where the Moon and the Sun are. And rather than use an external ephemeris for these bodies, I have elected to co-integrate the motions of the Sun and the Moon in the ECI reference frame.

Let's look at the first of these functions. The full routine is given as:

Code:
def a_mon(Q_BCI, mjd):

# unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

# gravitational acceleration of the Moon due to the Earth
Q_mon_BCBF = BCI2BCBF(Q_mon_BCI, mjd)
A_mon_BCBF = g( Q_mon_BCBF, "Earth")
A_mon_BCI  = BCBF2BCI(A_mon_BCBF, mjd)

# correct for reduced mass of the Earth-Moon system
A_mon_BCI  = (1 + mu_dict["Moon"] / mu_dict["Earth"]) * A_mon_BCI

# add the tidal term due to the Sun's gravitational influence
A_mon_BCI += g( Q_sun_BCI, "Sun") - g( Q_sun_BCI - Q_mon_BCI, "Sun")

return A_mon_BCI

This function takes a list of postion vectors 'Q_BCI' of the Moon, the Sun and the satellite in the ECI reference frame on a given date, 'mjd'. It returns the gravitational acceleration of the Moon for that date. Let's go through this code step-by-step.

First, we have:

Code:
    # unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

All this does is unpack the list of position vectors contained in 'Q_BCI'. The first entry in the list is the position vector of the Moon in the ECI reference frame and we assign it to 'Q_mon_BCI'. The second entry is the position vector of the Sun in the ECI reference frame and we assign it to 'Q_sun_BCI'. And the third entry is the position vector of the satellite in the ECI reference frame and we assign it to 'Q_veh_BCI'.

Next we have:

Code:
  # gravitational acceleration of the Moon due to the Earth
Q_mon_BCBF = BCI2BCBF(Q_mon_BCI, mjd)
A_mon_BCBF = g( Q_mon_BCBF, "Earth")
A_mon_BCI  = BCBF2BCI(A_mon_BCBF, mjd)

The first line of this code snippet converts the position vector of the Moon in the ECI reference frame to the corresponding position vector in the ECEF reference frame. Once in the ECEF reference frame, we calculate the acceleration of the Moon due to the gravitational attraction of the Earth (including the acceleration components due to the Earth's non-spherical gravitational perturbations). The result is an acceleration vector in the ECEF reference frame. So, finally, we need to convert this back to the ECI reference frame so that we can use it in our trajectory integration calculations.

To calculate the acceleration of the Moon in the ECI reference frame, we have to make an adjustment for the fact that the Earth and the Moon are rotating around the Earth Moon Barycentre (EMB). This adjustment is captured in the code snippet:

Code:
  # correct for reduced mass of the Earth-Moon system
A_mon_BCI  = (1 + mu_dict["Moon"] / mu_dict["Earth"]) * A_mon_BCI

The derivation of this adjustment was outlined in the post "The ECI reference frame and tidal forces".

In addition, to making this barycentre adjustment, we also need to take into account the tidal force of the Sun acting on the Moon. We do this in the code snippet:

Code:
  # add the tidal term due to the Sun's gravitational influence
A_mon_BCI += g( Q_sun_BCI, "Sun") - g( Q_sun_BCI - Q_mon_BCI, "Sun")

This tidal force adjustment follows the adjustment as also set out in "The ECI reference frame and tidal forces". The astute reader may observe that we are using position vectors in the Sun's BCI reference frame here rather than in the Sun's BCBF reference frame. We can get away with this here because in Orbiter, the Sun's gravitational field doesn't include non-spherical terms so it makes no difference whether or not we work in the BCBF or BCI reference frames. If in Orbiter, the Sun were modeled with non-spherical gravity terms, then we would have to work a little harder to transform to the appropriate BCBF reference frame before calculating the tidal contributions. We then return the gravitational acceleration of the Moon in the ECI reference frame:

Code:
    return A_mon_BCI

Having calculated the acceleration of the Moon in the ECI reference frame in this way, we then do much the same to calculate the acceleration of the Sun in the ECI reference frame. The code for doing this parallels that of the calculations for the Moon, except that we now need to add the tidal contribution due to the Moon rather than that of the Sun.

Code:
def a_sun(Q_BCI, mjd):

# unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

# gravitation acceleration of the Sun due to the Earth
Q_sun_BCBF = BCI2BCBF(Q_sun_BCI, mjd)
A_sun_BCBF = g( Q_sun_BCBF, "Earth")
A_sun_BCI  = BCBF2BCI(A_sun_BCBF, mjd)

# correct for reduced mass of the Earth-Sun system
A_sun_BCI  = (1 + mu_dict["Sun"] / mu_dict["Earth"]) * A_sun_BCI

# add the tidal term due to the Moon's gravitational influence
A_sun_BCI += g( Q_mon_BCI, "Moon") - g( Q_mon_BCI - Q_sun_BCI, "Moon")

return A_sun_BCI

Again, in calculating the tidal terms acting on the Sun due to the Moon, we can take advantage of the fact that Orbiter has no non-spherical gravitational terms for the Moon and so we can elect to calculate the gravitational acceleration directly in the Moon's BCI reference frame rather that the Moon's BCBF reference frame. If this were not the case, we weould have to include a few additional coordinate transformations at this juncture.

Finally, we can calculate the gravitational acceleration of the satellite in the ECI reference frame:

Code:
def a_veh(Q_BCI, mjd):

# unpack the state vector
Q_mon_BCI = Q_BCI[0]
Q_sun_BCI = Q_BCI[1]
Q_veh_BCI = Q_BCI[2]

# gravitation acceleration of the Vehicle due to the Earth
Q_veh_BCBF = BCI2BCBF(Q_veh_BCI, mjd)
A_veh_BCBF = g(Q_veh_BCBF, "Earth")
A_veh_BCI  = BCBF2BCI(A_veh_BCBF, mjd)

# add the tidal terms due to the gravitational influece of the Moon and the Sun
A_veh_BCI += g( Q_mon_BCI, "Moon") - g( Q_mon_BCI - Q_veh_BCI, "Moon")
A_veh_BCI += g( Q_sun_BCI, "Sun" ) - g( Q_sun_BCI - Q_veh_BCI, "Sun" )

return A_veh_BCI

Again, the procedure for doing this is the same except that:

* we do not need to make a barycentre adjustment; and

* we need to include the tidal forces due to both the Moon and the Sun.

Again, we take advantage of the spherical symmetry of Orbiter's assumed gravitational fields for the Sun and the Moon to simplify the calculation of the tidal terms.

Now we come onto the heart of the integrator: building the core secoond-order integrator:

Code:
def Integrator2(Q0, V0, mjd, dt):

# unpack the state vectors
[Q0_mon, Q0_sun, Q0_veh]  = Q0
[V0_mon, V0_sun, V0_veh]  = V0

# first 'kick' step
Va_mon  = V0_mon + 0.5 * dt * a_mon(Q0, mjd)
Va_sun  = V0_sun + 0.5 * dt * a_sun(Q0, mjd)
Va_veh  = V0_veh + 0.5 * dt * a_veh(Q0, mjd)

# drift step
Q1_mon  = Q0_mon + dt * Va_mon
Q1_sun  = Q0_sun + dt * Va_sun
Q1_veh  = Q0_veh + dt * Va_veh

# repack the position vector and update time
Q1  = [Q1_mon, Q1_sun, Q1_veh]
mjd = mjd + dt / 86400.0

# second 'kick' step
V1_mon  = Va_mon + 0.5 * dt * a_mon(Q1, mjd)
V1_sun  = Va_sun + 0.5 * dt * a_sun(Q1, mjd)
V1_veh  = Va_veh + 0.5 * dt * a_veh(Q1, mjd)

# repack the velocity vector
V1 = [V1_mon, V1_sun, V1_veh]

return Q1, V1, mjd

This type of integrator is a 'kick-drift-kick' integrator. It consists of three steps. The first is a momentum 'kick' that changes the speed of the bodies; the second is a drift at constant speed after having received the momentum kick; and the third is a second momentum 'kick' that makes a final adjustment to the velocity of the bodies.

In this integrator, the motion of three bodies are modeled: the Moon, the Sun and the satellite. The position and velocities of all three of these bodies are stored in two vector lists, Q0 and V0. The first action we take is to unpack these lists:

Code:
    # unpack the state vectors
[Q0_mon, Q0_sun, Q0_veh]  = Q0
[V0_mon, V0_sun, V0_veh]  = V0

Now we have explicit names for each of the state vectors of each of the bodies that we are modeling. The next action is to perform the firs momentum 'kick' for each body:

Code:
    # first 'kick' step
Va_mon  = V0_mon + 0.5 * dt * a_mon(Q0, mjd)
Va_sun  = V0_sun + 0.5 * dt * a_sun(Q0, mjd)
Va_veh  = V0_veh + 0.5 * dt * a_veh(Q0, mjd)

Here we update the velocity vectors of the Moon, the Sun and the satellite using the force calculation routines described earlier. The factor of 0.5 appears because we want the 'kick' to apply for only half a time-step. The other half a time-step will come from the second momentum 'kick' performed at the end of the integration time-step.

Next we perform the 'drift' calculation:

Code:
    # drift step
Q1_mon  = Q0_mon + dt * Va_mon
Q1_sun  = Q0_sun + dt * Va_sun
Q1_veh  = Q0_veh + dt * Va_veh

At this point, we take the opportunity to repack the final position vectors of the bodies and to update the MJD:

Code:
    # repack the position vector and update time
Q1  = [Q1_mon, Q1_sun, Q1_veh]
mjd = mjd + dt / 86400.0

And, finally, we perform the second momentum 'kick' using the previously updated positions and velocities:

Code:
    # second 'kick' step
V1_mon  = Va_mon + 0.5 * dt * a_mon(Q1, mjd)
V1_sun  = Va_sun + 0.5 * dt * a_sun(Q1, mjd)
V1_veh  = Va_veh + 0.5 * dt * a_veh(Q1, mjd)

The velocity vector is re-packed into its list form and this is returned from the integration routine along with the updated MJD.

The result of this sequence of calculations is a second-order (self-adjoint and symplectic) integration scheme. That's all fine, but we can do better. And that's where the next bit of the code comes in:

Code:
def Integrator4(Q0, V0, mjd, dt):

w0 = +1.35120719195965763
w1 = -1.70241438391931527

Qa, Va, mjd = Integrator2(Q0, V0, mjd, dt * w0)
Qb, Vb, mjd = Integrator2(Qa, Va, mjd, dt * w1)
Q1, V1, mjd = Integrator2(Qb, Vb, mjd, dt * w0)

return Q1, V1, mjd

It turns that any self-adjoint and symplectic integrator accurate to second-order can be upgraded to a fourth-order (self-adjoint and symplectic) integrator by this algorith,. The coefficients 'w0' and 'w1' are special magic numbers that divide the time-step into three pieces of unequal length. In each of these three-time step the second-order integrator is called. In total, the functions used to calculate the accelerations of the Moon, the Sun and the satellite are each called six times.

Now we have a fourth-order integrate. This too is good, but we can do better. And that's where the final bit of the code comes in:

Code:
def Integrator6(Q0, V0, mjd, dt):

w0 = +1.17467175808936338
w1 = -1.34934351617872677

Qa, Va, mjd = Integrator4(Q0, V0, mjd, dt * w0)
Qb, Vb, mjd = Integrator4(Qa, Va, mjd, dt * w1)
Q1, V1, mjd = Integrator4(Qb, Vb, mjd, dt * w0)

return Q1, V1, mjd

Just as there is a general method of upgrading any self-adjoint and symplectic second-order integrator into a fourth-order integrator, this code snippet provides a general method for upgrading a fourth-order integrator into a sixth-order integrator. Again, the coefficients are 'magic' numbers that divide the the time-step into three non-equal parts; and in each of these three time-steps, the fourth-order integrator is called. So now the functions used to calculate the body accelerations are called a total of eighteen times (3 x 6) in each time-step.

As one might expect, there is a whole series of similar general methods that upgrade a sixth-order integrator into an eigth-order integrator; and an eigth-order integrator into a tenth-order integrator; and so on. But the computational cost of these order upgrades is exponential and, for the problem in hand, gives rapidly diminishing returns. Even with time-steps of 30 seconds, the sixth-order integrator integrates the trajectory of a satellite with near machine precision and using progressively higher-order gives almost no further benefit. So, we stop the ugrade process at this point.

The final integration routine that we call to integrate the equations of motion of a satellite in the ECI reference frame is 'Integrator6', then. In turn, 'Integrator6' calls 'Integrator4' which, in turn calls 'Integrator2'. 'Integrator2' uses 'a_mon', 'a_sun' and 'a_veh' to calculate the gravitational accelerations of the Moon, the Sun and the satellite. In each of those functions, we transform from the ECI reference frame, calculate the non-spherical contributions; and transform back to the ECI reference frame. By calling 'Integrator6', we do a lot of computational work. But that computational effort yields accurate results.

In the next post, I will start using the integrator to test its accuracy and then use the sixth-order integrator to solve a number of rendezvous problems in Orbiter.

Last edited: