
Math & Physics Mathematical and physical problems of space flight and astronomy. 

Thread Tools 
07092018, 01:42 PM  #1 
Orbinaut

Putting it all together  a sixthorder integrator
This note is going to be about building a highfidelity integration engine for Orbiter that takes into account the gravitational influence of other bodies (tidal forces) and the 'nonspherical 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 Earthcentred 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 nonrotating reference frame centred on the Earth's. The ECI reference frame (or more generally, the Bodycentred Inertial (BCI) reference frame) is noninertial. 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 'nonspherical gravity sources'. Typically, these are the assorted spherical harmonics whose coefficients are given as , 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 terms.) This gravitational model for the Earth can also be generalised to any gravitating body. These nonspherical gravity contributions are most conveniently expressed in the Earthcentred Earthfixed (ECEF)  or more generally the Bodycentred Bodyfixed (BCBF) reference frame  which is a reference frame that corotates 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 nonspherical 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 highfidelity 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 linebyline 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.67259e11 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.6269e6, 2.51e6, 1.60e6, 0.15e6], "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[n1] + (1n) * P[n2]; P[n] /= n DP[n] = (2*n  1) * P[n1] + (2*n  1) * xi * DP[n1] + (1n) * DP[n2]; 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 EarthMoon 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 EarthSun 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 (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 Next, we have: Code:
# the universal gravitatinal constant (as per Orbiter documentation) G = 6.67259e11 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.6269e6, 2.51e6, 1.60e6, 0.15e6], "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] } Next, we define a number of Python 'dictionaries' (or lookup tables): * 'mu' the standard gravitational parameters for the Earth, the Moon and the Sun; * 'R' and 'J' the reference radius and nonspherical 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 Earthcentred 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 ): 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): 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 EarthMoon 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 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] 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) 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 EarthMoon system A_mon_BCI = (1 + mu_dict["Moon"] / mu_dict["Earth"]) * A_mon_BCI 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") Code:
return A_mon_BCI 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 EarthSun 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 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 * 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 secoondorder 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 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 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) 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 Code:
# repack the position vector and update time Q1 = [Q1_mon, Q1_sun, Q1_veh] mjd = mjd + dt / 86400.0 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 result of this sequence of calculations is a secondorder (selfadjoint 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 Now we have a fourthorder 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 As one might expect, there is a whole series of similar general methods that upgrade a sixthorder integrator into an eigthorder integrator; and an eigthorder integrator into a tenthorder 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 timesteps of 30 seconds, the sixthorder integrator integrates the trajectory of a satellite with near machine precision and using progressively higherorder 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 nonspherical 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 sixthorder integrator to solve a number of rendezvous problems in Orbiter. Last edited by MontBlanc2012; 07092018 at 11:38 PM. 

Thread Tools  


Quick Links  Need Help? 