Construct a Halo orbit around L2 of the Earth-Moon system

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
OK, it's been some months since my last substantial post on Lissajous and Halo orbits; and, to be frank, it's also been some months since I last thought seriously about this topic. However, when I last wrote a post on Halo Orbits (way back on 30 October, 2018), I ended up by saying that (now that I had ironed out the basic maths of Lissajous and Halo orbits that:

I'm going to work through a number of examples and attach them to the relevant posts - so, even though the arcane maths may not be entirely clear, the mechanics of constructing the state vectors of a Halo orbit should be.

After that, I'm going to show how to 'lock' a vessel onto a Halo orbit, starting from the initial state vectors calculated using this approach. This will serve as a proxy for active station keeping and provide Halo orbit targets for Orbinauts to rendezvous with, etc.

Well, even though three months have passes, these still seems to me to be a worth goals, so in this post, I'm going to start the first of these - namely working through a number of specific examples of constructing Halo orbits - for the Earth-Moon system; for the Sun-Earth system and (just for kicks) for the the Sun-Mars system. In each case, the result of this exercise will be the calculation of Orbiter-friendly state vectors (position and velocity) that will allow one to place a spacecraft in a prescribed Halo orbit using Orbiter's Scenario Editor (or, if one wants to automate the process, via Orbiter's api commands).

As I've repeatedly mentioned in previous posts on Lissajous orbits, even though we can calculate the initial position and velocity of a spacecraft in a prescribed Halo orbit, we do not expect that spacecraft to remain on that Halo orbit because (unlike the usual two-body Keplerian dynamics), dynamics in the vicinity of these Halo orbits is unstable and so some form of active station-keeping is needed to keep them on the prescribed Halo orbit. Once I've worked through the various examples of setting up the initial conditions for Halo orbits, I'll later present a simple Lua script for performing a hot-wired version of this active station-keeping that will, in effect, lock the spacecraft onto the prescribed Halo (or more generally, Lissajous) orbit and keep it there indefinitely


General process
In each case, the general process for construction of a Halo orbit initial conditions is exactly the same.

1. Set up some basic information: For a given MJD, capture relevant information about the elliptical orbit of the primary gravitating body around the secondary gravitating body

1. Choose the scale parameters of the Haol orbt: A Halo orbit is a specific case of a general Lissajous orbit where the frequency of oscillation in the in-plane 'x-y' direction is the same as the frequency of oscillation in the out-of-plane 'z' direction. In turn, these frequencies were a function of the amplitude of oscillations in the 'x-y' and 'z' directions, so that by choosing a scale in, say, the 'z' direction; the spreadsheet widget in the the spreadsheet of the post A digression on Halo orbits gives the required amplitude of oscillation in the 'x-y' direction. By choosing a different scale parameter, we choose a different Halo orbit. In fact, there are infinitely many such orbits that we can choose from (just as there are infinitely many circular orbits in two-body Keplerian orbital dynamics). Ultimately, the choice of scale depends on the operational objectives of the Halo orbit.

2. Having determined the scale parameters of the Halo orbit using the aforementioned spreadsheet widget, we then set the values of the six parameters that we have used to parameterise our Halo orbit. These six parameters are:

- [MATH]d\alpha_1[/MATH] - the scale parameter in the 'x-y' direction
- [MATH]d\alpha_2[/MATH] - the scale parameter in the 'z' direction
- [MATH]e[/MATH] - the orbital eccentricity of the secondary gravitating body around the primary (e.g., the eccentricity of the Moon around the Earth)
- [MATH]\nu[/MATH] - the true anomaly of the secondary gravitating body around the primary (e.g., the true anomaly of the Moon around the Earth
- [MATH]\phi_1[/MATH] - a phase angle in the 'x-y' direction
- [MATH]\phi_2[/MATH] - a phase angle in the 'z' direction

So, let's go through this parameter list and see where we stand. By using the Halo spreadsheet widget, we've set the scale parameter in the 'z' direction and from this determined the scale parameter in the 'x-y' direction. The orbital eccentricity and true anomaly and orbital eccentricity have been calculated as part of the first step. So, that's four parameters done. And, finally we set both phase angles to 0 (for a 'northern' Halo orbit); and the 'x-y' phase angle to 0 and the 'z' phase angle to $\pi$ (180 degrees) for a 'southern' Halo orbit). Easy!

3. Having determined the orbital parameters, we then have all the information that we need to uniquely define a Halo orbit. Then, we use the datasets set out in previous posts (e.g., Dialing-up arbitrary Lagrange Point orbits to calculate the state vectors of the spacecraft in the somewhat byzantine rotating-pulsating coordinate system of the ER3BP (Elliptical Restricted Three-Body Problem).

4. Then, we transform the state vectors from the ER3BP rotating-pulsating coordinate system to the Orbiter's standard global coordinate system.

5. Finally, we set the spacecraft's state vectors to these calculated values using Orbiter's Scenario Editor; or by using Orbiter's api.

Although this process sounds quite complicated, it is in fact a highly mechanical process and can be automated easily. But we'll work through it step-by-step for the Earth Moon system below.


Example 1 - a Halo orbit for Earth-Moon L2

OK, then, time to do some specific examples - starting with the Earth-Mon system. The most commonly discussed Halo orbits in the Earth-Moon system are those around the L2 Lagrange point. So, the first thing to do is to choose an MJD and determine the orbital eccentricity and the true anomaly of the Moon with respect to the Earth on that date. Here's some Lua code that interrogates Orbiter's system to calculate these quantities:


Code:
-- set the system MJD
mjd = 58511.000000
oapi.set_simmjd(mjd)

-- set the gravitational constant
G       = 6.67259e-11
primary = oapi.get_objhandle("Earth")
secondary = oapi.get_objhandle("Moon")

-- get the current state vectors and gravitational parameter of the primary gravitating body
q_pri   = oapi.get_globalpos(primary)
p_pri   = oapi.get_globalvel(primary)
m_pri   = oapi.get_mass(primary)
mu_pri  = G * m_pri
	
-- get the current state vectors and gravitational parameter of the secondary gravitating body
q_sec   = oapi.get_globalpos(secondary)
p_sec   = oapi.get_globalvel(secondary)
m_sec   = oapi.get_mass(secondary)
mu_sec  = G * m_sec

-- calculate the orbital eccentricity vector
mu      = mu_pri + mu_sec
X       = vec.sub(q_sec, q_pri)
V       = vec.sub(p_sec, p_pri)
evec    = vec.mul(vec.dotp(V, V) / mu - 1 / vec.length(X), X)
evec    = vec.sub(evec, vec.mul(vec.dotp(V, X) / mu, V))

-- calculate the true anaomaly
nu      = math.acos(vec.dotp(evec, X) / vec.length(evec) / vec.length(X))
if vec.dotp(X, V) < 0 then nu = 2 * math.pi - nu end

-- print out the results
file = io.open("testmb2012.txt", "w")
io.output(file)
file:write(mjd,"\n")
file:write("evec = \{", evec.x, ", ", evec.y, ", ", evec.z, "\}\n")
file:write("ecc  = ", vec.length(evec), "\n")
file:write("nu   = ", nu, "\n")
file:write("q_pri:  ", q_pri.x, "  ", q_pri.z, "  ", q_pri.y, "\n")
file:write("p_pri:  ", p_pri.x, "  ", p_pri.z, "  ", p_pri.y, "\n")
file:write("q_sec:  ", q_sec.x, "  ", q_sec.z, "  ", q_sec.y, "\n")
file:write("p_sec:  ", p_sec.x, "  ", p_sec.z, "  ", p_sec.y, "\n")
file:write("mu_pri: ", mu_pri, "\n")
file:write("mu_sec: ", mu_sec, "\n")
file:close()

(In this script, we have set the the primary gravitating body to the Earth; and the secondary gravitating body to the Moon. Of course for other systems - e.g. the Sun-Mars system, we would set the primary to the Sun and the secondary to Mars.)

If one runs this Lua script in Orbiter, and then opens the designated output file 'testmb2012.txt', one should find the following output:

Code:
58511
evec = {-0.02985654146866, 0.00016241939735106, 0.056132311490403}
ecc  = 0.063578894625387
nu   = 1.7602449698559
q_pri:  -89818288265.986  117996567403.37  -13177279.81454
p_pri:  -24140.309588757  -18232.438540029  1.4286445780124
q_sec:  -90114909520.995  117757041901.07  21314097.759869
p_sec:  -23550.70413747  -19067.747671548  -13.045629242125
mu_pri: 3.9860043996887e+014
mu_sec: 4902794932976

This output file lists the MJD that we have chosen (58511.0000); the orbital eccentricity of the Moon relative to the Earth (0.063578895875671); the true anomaly of the Moon relative to the Earth (1.7602450571985 radians); the state vectors (position and velocity) of the primary and secondary bodies (Earth and Moon in this instance in Orbiter's global coordinate system (albeit with the 'y' and 'z' components transposed to convert to a more conventional right-handed coordinate system); and the gravitational parameters for the primary and secondary gravitating bodies in the ER3BP.

Time now to calculate the scale parameters [MATH]d\alpha_1[/MATH] and [MATH]d\alpha_2[/MATH]. these are determined using the spreadsheet widget as discussed in A digression on Halo orbits. To do this open the Halo.xlsx spreadsheet and move to the Excel tab labelled 'EM L2' (because we want to calculate the parameters of a Halo orbit around the Earth-Moon L2 Lagrange point). Next, enter the value for the eccentricity calculated above, i.e., [MATH]e = 0.063578894625387[/MATH] into cell B18. Then, set the cell B19 to a suitable value for the 'z' direction scale parameter, [MATH]d\alpha_2[/MATH]. Because the Halo orbit that we are going to construct has no particular mission objectives, there is a degree of arbitrariness to this. However, we don't want to make the value too large and a value of around 0.03 is probably a convenient size. (So, what does a value of 0.03 mean? Well, it's a length measured as a multiple of the Earth-Moon distance. Since the Earth-Moon distance is around 385,000 km, a value of 0.03 is roughy 0.03 x 385,000 or about 11,550 km. Which means that the vertical 'z' direction size of the Halo obit is about twice that - or around 23,000 km. This isn't a small orbit.) From the spreadsheet we can now read of the scale parameter in the 'x-y' direction. This should be [MATH]\alpha_1 = 0.03782520[/math]. This means that the end-to-end 'x-y' length of the Halo orbit is around 65,000 km. Again, this isn't small.

Now, then, we have four of the six parameters that are needed to describe the halo orbit - [MATH]e[/MATH], [MATH]\nu[/MATH], [MATH]d\alpha_1[/MATH] and [MATH]d\alpha_2[/MATH]. Just two more to go: the two phase parameters, [MATH]\phi_1[/MATH] and [MATH]\phi_2[/MATH] - but these are very easy. Because we have no particular mission objectives in mind for this Halo orbit, we can set [MATH]\phi_1 = 0[/MATH] and supposing that we want a 'northern' Halo orbit, we set [MATH]\phi_2 = 0[/MATH] as well. (If we wanted a 'southern' Halo orbit, we would set [MATH]\phi_2 = \pi[/MATH].

Now we are in a position to calculate the state vectors of the Halo orbit (in the ER3BP's rotating=pulsating coordinates) for MJD = 58511.0000 for the L2 Lagrange point of the Earth-Moon system. To do this, we run a number of Python scripts posted earlier in Dialing-up arbitrary Lagrange Point orbits. The Python script used to read the Lagrange point data sets and set up the function to calculate the state vectors is:

Code:
import openpyxl
import math
import numpy as np

wb = openpyxl.load_workbook('path to file/EM L2 - Lissajous - ER3BP.xlsx')

# extract the data from the spreadsheet
wsgeneral = list(wb['general'    ].rows)
wsomega   = list(wb['omega-data' ].rows)
wslambda  = list(wb['lambda-data'].rows)
wsx       = list(wb['x-data'     ].rows)
wsy       = list(wb['y-data'     ].rows)
wsz       = list(wb['z-data'     ].rows)

mu    = wsgeneral[0][1].value
alpha = wsgeneral[1][1].value


def getOmega(e, da1, da2):
    
    omeg = 0
    for idx in range(1, len(wsomega)):
        x     = wsomega[idx]
        omeg += x[4].value * e**x[1].value * da1**x[2].value * da2**x[3].value
    
    return omeg


def getLambda(e, da1, da2):
    
    lamb = 0
    for idx in range(1, len(wslambda)):
        x     = wslambda[idx]
        lamb += x[4].value * e**x[1].value * da1**x[2].value * da2**x[3].value
        
    return lamb


def pos(e, da1, da2, nu, phi1, phi2):
    
    omeg = getOmega (e, da1, da2)
    lamb = getLambda(e, da1, da2)    
   
    theta1 = omeg * nu + phi1
    theta2 = lamb * nu + phi2
    
    eta = 0
    for idx in range(1, len(wsx)):
        x      = wsx[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        eta   += x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.cos(angle)
        
    kappa = 0
    for idx in range(1, len(wsy)):
        x      = wsy[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        kappa += x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.sin(angle)
        
    sigma = 0
    for idx in range(1, len(wsz)):
        x      = wsz[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        sigma += x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.cos(angle)
        
    deta = 0
    for idx in range(1, len(wsx)):
        x      = wsx[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        w      = x[4].value + omeg * x[5].value + lamb * x[6].value
        deta  += - w * x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.sin(angle)
        
    dkappa = 0
    for idx in range(1, len(wsy)):
        x      = wsy[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        w      = x[4].value + omeg * x[5].value + lamb * x[6].value
        dkappa+= + w * x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.cos(angle)
        
    dsigma = 0
    for idx in range(1, len(wsz)):
        x      = wsz[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        w      = x[4].value + omeg * x[5].value + lamb * x[6].value
        dsigma+= -w * x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.sin(angle)
        
    return [[eta + alpha, kappa, sigma], [deta, dkappa, dsigma]]

Here, we the file "EM L2 - Lissajous - ER3BP.xlsx" can be downloaded from the .zip file in Dialing-up arbitrary Lagrange Point orbits. If we run this script by invoking:


Code:
## use the parameter values calculated earlier
e    = 0.063578894625387
nu   = 1.760244969
da1  = 0.03782520
da2  = 0.03
phi1 = 0.0
phi2 = 0.0

p = pos(e, da1, da2, nu, phi1, phi2)
print(p)

we should find that this routine returns:


Code:
[[1.1102069240218368, 0.005087416690776363, -0.051590690909498756],
 [0.0025792635799946413, 0.19776761282432434, 0.006101524551592115]]

This is a list consisting of two 3-vectors. The first is [MATH][1.1102069240218368, 0.005087416690776363, -0.051590690909498756][/MATH] which is the position of the spacecraft on the prescribed Halo orbit at MJD = 58511.000000 in the ER3BP rotating-pulsating coordinate system; and the second is [MATH][0.0025792635799946413, 0.19776761282432434, 0.006101524551592115][/MATH] which is the velocity of the spacecraft in the ER3BP rotating-pulsating coordinate system. However, we can't directly use this information to construct the state vectors of a vessel in the prescribed Halo orbit in Orbiter. Instead, we have to apply a final script to map the coordinates back into Orbiter's coordinate system. To do this, we have to invoke one further Python script:

Code:
# input the epoch MJD
mjd = 58511.0000

# input the state vectors of the primary gravitating body - right-handed coordinate system
Qpri  = np.array([-89818288265.957,  117996567403.43,  -13177279.813415])
Ppri  = np.array([-24140.309588727,  -18232.43854004,  1.4286445767117])

# input the state vectors of the secondary gravitating body - right-handed coordinate system
Qsec  = np.array([-90114909520.951,  117757041901.14,  21314097.759466])
Psec  = np.array([-23550.70413747,  -19067.747671548,  -13.045629242125])

# input the gravitational constants
mupri = 3.9860043996887e+014
musec = 4902794932976.0


mu    = mupri + musec

def getECIcoords(a, e, nu, pos, mu):
    
    [[eta, kappa, sigma], [deta, dkappa, dsigma]] = pos
    cs = math.cos(nu)
    sn = math.sin(nu)

    f  = a * (1 - e*e) / (1 + e * cs)
    x  = f * (eta * cs - kappa * sn)
    y  = f * (eta * sn + kappa * cs)
    z  = f * sigma
    
    g  = math.sqrt(mu / a / (1-e*e))
    vx = g * (-sn*eta   - (e + cs)*kappa + (1 + e*cs)*(cs*deta   - sn*dkappa))
    vy = g * (-sn*kappa + (e + cs)*eta   + (1 + e*cs)*(cs*dkappa + sn*deta  ))
    vz = g * ( e*sn*sigma + (1 + e*cs)*dsigma)
    
    return [[x, y, z], [vx, vy, vz]]


# clculate the centre of mass and the centre of velocity
com = (mupri * Qpri + musec * Qsec) / mu
cov = (mupri * Ppri + musec * Psec) / mu

r   = Qsec - Qpri
v   = Psec - Ppri
R   = np.linalg.norm(r)

# calculate a, e, nu
evec= r * np.dot(v,v)/mu - v * np.dot(r,v)/mu - r/R
e   = np.linalg.norm(evec)
a   = mu / (2*mu / R - np.dot(v,v))
nu  = math.acos(np.dot(evec,r)/e/R)
if np.dot(r,v)<0:
    nu = 2 * math.pi - nu
    
xhat= evec / e
zhat= np.cross(r, v)
zhat= zhat / np.linalg.norm(zhat)
yhat= np.cross(zhat,xhat)

# now calculate the ECI state vectors relative to the primary - right-handed coordinate system
[[Qx, Qy, Qz], [Px, Py, Pz]] = getECIcoords(a, e, nu, p, mu)
Q = Qx * xhat + Qy * yhat + Qz * zhat + com - Qpri
P = Px * xhat + Py * yhat + Pz * zhat + cov - Ppri

# print the MJD
print(mjd)

# print the position vector relative to the primary - (Orbiter-freindly left-handed coordinate system, ecliptic, cartesian)
print([Q[0], Q[2], Q[1]])

# print the velocity vector relative to the primary - (Orbiter-friendly left-handed coordinate system, ecliptic, cartesian)
print([P[0], P[2], P[1]])


In this script, we have to enter the MJD (58511.0000), the state vectors of the primary and the secondary (the Earth and the Moon) in Orbiter's global coordinate system (with the 'y' and 'z' components swapped to convert to a right-handed coordinate system); and the gravitational parameters of the primary and the secondary. This information, together with the state vectors calculated earlier in the ER3BP reference frame provides sufficient information to calculate the state vectors in Orbiter's coordinate system. When run, we should obtain as output:

Code:
58511.0
[-333326298.9295349, 19007567.817135274, -271161377.07902527]
[790.6644000361339, -17.607214696246114, -1093.2669313970146]

This records the MJD for which the state vectors are valid; the vector [MATH][-333326298.9295349, 19007567.817135274, -271161377.07902527][/MATH] is the position vector (in metres) in an Earth-centric, ecliptic, cartesian coordinate system; and the vector [MATH][790.6644000361339, -17.607214696246114, -1093.2669313970146][/MATH] is the velocity vector (metres / second) again in an Earth-centric, ecliptic, cartesian coordinate system. These values can now be inserted directly into Orbiter's using the Scenario Editor or, if automated, by using relevant api commands.

For convenience, here is a simple Orbiter scenario that sets up a Delta Glider with the state vectors as calculated.

Code:
BEGIN_ENVIRONMENT
  System Sol
  Date MJD 58511.00000
END_ENVIRONMENT

BEGIN_FOCUS
  Ship GL-NT
END_FOCUS

BEGIN_SHIPS
GL-NT:DeltaGlider
  STATUS Orbiting Earth
  RPOS -333326298.9295349 19007567.817135274 -271161377.07902527
  RVEL 790.6644000361339 -17.607214696246114 -1093.2669313970146
  AROT -179.993 -38.127 -179.973
  RCSMODE 2
  AFCMODE 7
  PRPLEVEL 0:1.000000 1:0.988429
  NAVFREQ 0 0 0 0
  XPDR 0
  HOVERHOLD 0 1 0.0000e+000 0.0000e+000
  AAP 0:0 0:0 0:0
END
END_SHIPS


In summary

In summary, in this note, we've walked through the process of placing a vessel on a prescribed Halo orbit around the L2 Lagrange point in the Earth-Moon system. As always with Halo (and, more generally, Lissajous orbits), these orbits are unstable and small errors in the underlying perturbation approximations as well as solar effects will 'kick' the vessel of the Halo orbit. Consequently, active station-keeping is needed to keep it there - but writing a script that does this is beyond the scope of this post, although it's something I'll come back to later.

In the next post, and just to demonstrate that approach is quite general I'll go through the mechanics of this a second time for the Sun-Mars system and construct a Halo orbit around L1 of that system.
 
Top