# Numerical calculation of the Earth-Moon L2 halo orbit family

#### MontBlanc2012

##### Member
As many of the regular readers of this Maths & Physics forum will know, over the last six months I've been working on an exercise of 'mapping' the centre-manifold orbits (planar and vertical Lyapunov, Lissajous and halo orbits) of the L1 and L2 Lagrange points - particularly those of the Earth-Moon and the Sun-Earth systems.

So far I've posted only maps using a form of perturbation theory (known as Lindstedt-Poincaré theory). Perturbation theory only takes one so far, however in mapping some of these orbit families. In particular, the work so far posted only maps 'small' Halo orbits and does not map the whole of the Halo orbit families.

Over the weekend, I took a step further in this mapping exercise by using numerical integration (bootstrapped with the existing perturbation theory) to map the whole of the L1 and L2 northern and southern halo orbit families. The result is a set of Excel-based interpolation tables that allow one to calculate (for Orbiter 2016) the state vectors (i.e., position and velocity) of any point on any halo orbit for any time. This seems to me to be quite a useful thing so I'm going to share these results - together with some Python code that allow one to decode the Excel spreadsheets to calculate the orbit state vectors. - But I'll actually post this in a second instalment below this one in this thread.

The L2 southern halo orbit family
To whet the appetite, when for example all of the L2 southern halo orbits are pasted together, one developed the following picture of what the southern halo orbit family looks like: The orange surface is the halo orbit manifold (i.e., just the ensemble of all of the halo orbits in the Circular Restricted Three Body Problem - CR3BP). The white sphere depicts the position of the Moon and is drawn to scale. The black lines on the surface of the orange manifold show some specific examples of halo orbits that lie on this manifold.

The halo orbit marked in green is of particular interest. This is an example of something called a Near Rectilinear Halo Orbit (NRHO). Although it looks like a Keplerian elliptical orbit, it isn't - it's actually a halo orbit. But it is on a NRHO very much like this one that NASA propose to park LOP-G (the Lunar Orbital Platform Gateway station). Consequently, the spreadsheet interpolation tables provide a representation of NRHO orbits for use in Orbiter 2016 - and I'll show how to do this in a subsequent post in this thread. As you can see, though these NRHOs have a perilune that pass close to the Moon and there is an interesting problem to be investigated about how to launch from a point on the Moon and rendezvous with a LOP-G-like station parked in such an NRHO; and to complete the reverse transfer back to a base on the surface of the Moon.

The halo orbits bifurcate from the planar Lyapunov family
If we can say that the natural end point of the halo orbit family terminates in the NRHOs, what about their starting point. If we look at the same halo orbit manifold families from the other side, we can identify the originating planar Lyapunov orbit. The highlighted red 'halo orbit' is, in fact, the planar Lyapunov orbit (lying in the orbital plane of the Earth and the Moon) from which the whole of the L2 halo orbit family originates. In this image, the Moon is 'hidden' behind the halo orbit manifold, but the blue sphere depicts the location of the Earth (also drawn to scale in this image).

In a subsequent post in this thread, I'll return to posting the interpolation tables, assorted decoding Python scripts, and demonstrate how the information can be used to calculate the state vectors of these orbits for use in Orbiter 2016.

(N.B. the images in this post were prepared using Mathematica)

#### MontBlanc2012

##### Member
Well, in the preceding post in this thread I wrote:

In a subsequent post in this thread, I'll return to posting the interpolation tables, assorted decoding Python scripts, and demonstrate how the information can be used to calculate the state vectors of these orbits for use in Orbiter 2016.

So, that's exactly what I'm now going to do:

interpolation tables
The complete picture of the shape of the entire (southern) L2 halo orbit family is built up from a series of ~1300 individual snapshot high precision halo orbit - each of which provides an interpolation 'point' from which an arbitrary L2 halo orbit can be determined. These interpolation 'points' are stored in a 7 MB .xlsx spreadsheet file. This file (even when compressed) exceeds this forum's files size upload limits, so the spreadsheet can be downloaded from Dropbox here:

Raw halo orbit interpolation data file.

This spreadsheet consists of eight tabs. The first three are the interpolation points of the 'x', 'y' and 'z' position coordinates in the rotating coordinate system of the Circular Restricted Three Body Problem (CR3BP); the next three are the interpolation points of the 'x', 'y' and 'z' velocity coordinates in the same CR3BP reference frame. The seventh tab is an interpolation tab for orbital period. And the last tab is the set of 70 interpolation points of the underlying Gauss-Lobatto style variational integrator used to construct the 1300 or so individual halo orbits.

On the first six sheet, the Gauss-Lobatto interpolation points are given across the sheet and correspond to a single halo orbit. Vertically on each of these sheets data for the 1300 separately calculated halo orbits are listed.

This interpolation data format isn't excessively elegant but, meh, it works. And it provides a very accurate representation of solutions to the CR3BP model.

Python scripts to decode the raw interpolation data
So assuming that you want to (and have) download the interpolation file and stored it in a known location on your computer, the next step is to decode this interpolation data and construct a specific state vectors for a specific halo orbit of your choosing. Below is a description of a Python script that does that. Although one could automate this as an MFD and save you the trouble of implementing this code myself, my real interest in using Orbiter is in understanding the underlying maths and physics of some of these orbits. So I'm just going to provide a basic mechanism for decoding the interpolation data - and I'll leave it to others to develop better ways of accessing this data.

Let's break the Python script down into a number of chunks:

First, the snippet to load the required Python packages:

Code:
import math
import pandas as pd
import numpy  as np
from   scipy  import interpolate

The 'math' package is loaded so that Python knows something about sine, cosines and pi. the 'pandas' package is loaded because it provides a rather elegant interface to Excel files and provides a comprehensive set of data manipulation tools. 'numpy' is loaded because this is my 'go to' package for all things involving vectors and matrices. And the 'interpolate' module of the 'scripy' package is loaded because we need to carry out a spline interpolation to connect the 1300 or so separate halo orbits encoded in the data file.

Next, we load the data file into Python using 'pandas' and convert that data to numpy data arrays using the 'numpy' package:

Code:
file = 'EM L2 - halo - southern - 70.xlsx'

xl   = pd.ExcelFile(file)
x    = np.transpose(np.array((xl.parse('x'          )).values))
y    = np.transpose(np.array((xl.parse('y'          )).values))
z    = np.transpose(np.array((xl.parse('z'          )).values))
vx   = np.transpose(np.array((xl.parse('vx'         )).values))
vy   = np.transpose(np.array((xl.parse('vy'         )).values))
vz   = np.transpose(np.array((xl.parse('vz'         )).values))
T    = np.transpose(np.array((xl.parse('Period'     )).values))
itp  = np.transpose(np.array((xl.parse('Interp. pts')).values))

We now need to define two functions. The first of these is:

Code:
def getStateVectors(nu, period):

def lbp(i, xval, x):

prod = 1.0
for idx in range(len(x)):
if idx != i:
prod *= (xval - x[idx]) / (x[i] - x[idx])

return prod

if period >= np.min(T) and period <= np.max(T) and nu >= 0.0 and nu <= 1.0:

tval  = list(map(lambda idx: lbp(idx, nu, itp), range(70)))
xval  = list(map(lambda t: 1 * interpolate.interp1d(T, t, kind = 'cubic')(period), x ))
yval  = list(map(lambda t: 1 * interpolate.interp1d(T, t, kind = 'cubic')(period), y ))
zval  = list(map(lambda t: 1 * interpolate.interp1d(T, t, kind = 'cubic')(period), z ))
vxval = list(map(lambda t: 1 * interpolate.interp1d(T, t, kind = 'cubic')(period), vx))
vyval = list(map(lambda t: 1 * interpolate.interp1d(T, t, kind = 'cubic')(period), vy))
vzval = list(map(lambda t: 1 * interpolate.interp1d(T, t, kind = 'cubic')(period), vz))

Xvec  = list(np.dot([ xval,  yval,  zval], tval))
Vvec  = list(np.dot([vxval, vyval, vzval], tval))

return [Xvec, Vvec]

else:

return 'Input values out of bounds'

This function 'getStateVectors' takes two inputs - 'nu' and 'period' and converts to a set of state vectors in the rotating coordinate system of the CR3BP using a simple spline interpolation scheme, followed by a Legendre polynomial interpolation scheme to get to the final result. It is here that the all of the interpolation data points are converted into one seamless, smooth halo orbit manifold.

Two inputs into this function are required because the halo orbit manifold is a curved two-dimensional structure in three-dimensional space (think a curved piece of paper hanging in space). Any point on the surface is defined by two numbers, just as you would use to define the location of a point on a piece of paper.

So, what are these input parameters? The first is the orbital period. The maximum orbital period of a halo orbit on the manifold is 14.851 days; and the minimum orbital period of a halo orbit on the manifold is 7.683 days. The 14.851 day orbit corresponds to the planar Lyapunov orbit (given as the red orbit in the diagrams of the preceding post); and the 7.683 day orbit corresponds to the NRHO orbit (given as the green orbit in the diagrams of the preceding post. You can select a halo orbit with any orbital period that you like between these two values.

The second input parameter is a value called 'nu'. This identifies 'where' you are on the orbit in much the same away as the mean anomaly describes your location on a Keplerian orbit. 'nu' can take any value between 0 and 1. If 'nu' is set to 0 (or 1), you select a point at perilune (the point of closest approach to the Moon); if you set 'nu' to 0.5, you select apolune (the point of farthest approach to the Moon). If 'nu' is set to a value between 0.0 and 0.5, one is departing perilune; and if 'nu' is set to a value between 0.5 and 1.0, one is approaching perilune. For example, setting 'nu' to 0.95 would describe a position on the orbit that is close to (but not yet at) perilune.

A second function also needs to be defined - one that converts from the rotating coordinate system of the CR3BP to the normal rectilinear coordinate system used by Orbiter 2016. The function that does this is given below:

Code:
def getMCIcoords(frac, period):

[[eta, kappa, sigma], [deta, dkappa, dsigma]] = getStateVectors(frac, period)
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(G / 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]]

Using this function requires that a few things - the Moon's orbital eccentricity; its semi-major axis, and its true anomaly in its orbit around Earth be available as global quantities.

OK, so now to tie all of this together and calculate Orbiter 2106 compatible state vectors for a chosen halo orbit, we must first provide Python with some data about the current time and position of the Earth and the Moon as well as the gravitational parameters of the Earth and the Moon.

Code:
# input the epoch MJD
mjd = 52013.754909351

# input the state vectors of the Earth - right-handed coordinate system
Qe  = np.array([-136757075937.97, -63826187339.785, 20730124.687839])
Pe  = np.array([12032.791283836, -27150.845638622, 0.8388150790952])

# input the state vectors of the Moon - right-handed coordinate system
Qm  = np.array([-136653091584.32, -64212987439.536, 17225090.561136])
Pm  = np.array([12980.539253522, -26932.733069022, -85.194965428945])

# input the gravitational constants for the Earth and the Moon
GE  = 3.986004418e14
GM  = 4.9048695e12

All of this information can be gathered from Orbiter by interrogating Orbiter's integrator using the API commands for either Lua or C. Having loaded this basic information into Orbiter, we now calculate 'a', 'e', and true anomaly for the Moon:

Code:
# calcuale a, e, and nu
G   = GE + GM
muE = GE / (GE + GM)
muM = GM / (GE + GM)

com = muE * Qe + muM * Qm
cov = muE * Pe + muM * Pm

r   = Qm - Qe
v   = Pm - Pe
R   = np.linalg.norm(r)

evec= r * np.dot(v,v)/G - v * np.dot(r,v)/G - r /R
e   = np.linalg.norm(evec)
a   = G / (2*G / 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

These are standard calculations for calculating these quantities and are described in the orbiter.pdf file that comes standard with Orbiter.

Now, finally, we can use all of this information to calculate the state vectors of our halo orbit. So, first, we have two define out two parameters that determine our halo orbit and point on the halo orbit:

Code:
# define the halo orbit
frac = 0.967
period = 7.7

This says that we want to select a halo orbit with a period of 7.7 days; and that we want to be close to perilune. (An orbit with a period of this length is a NRHO orbit.) In addition, because we choose a 'nu' value of 0.967 (here called 'frac'), we choose a point on the orbit that is just a few hours away from perilune.

So, we've provided basic information about the date, location of the Earth and Moon and, of course, the point on our chosen halo orbit, it's time to put all of this together and calculate the state vectors of the vessel at that point on the chosen halo orbit:

Code:
# now calculate the MCI state vectors - right-handed coordinate system
period = period * 2 * math.pi / 27.32
[[Qx, Qy, Qz], [Px, Py, Pz]] = getMCIcoords(frac, period)

xhat= evec / e
zhat= np.cross(r,v)
zhat= zhat / np.linalg.norm(zhat)
yhat= np.cross(zhat,xhat)

Q = Qx * xhat + Qy * yhat + Qz * zhat + com - Qm
P = Px * xhat + Py * yhat + Pz * zhat + cov - Pm

# print the position vector - left-handed coordinate system
print('Moon-centric, ecliptic, Orbiter 2016 left-handed frame')
print('-----------------------------------------------')
print()
print('MJD:   ', mjd)
print()
print('Position x-coordinate (m):    ', Q)
print('Position y-coordinate (m):    ', Q)
print('Position z-coordinate (m):    ', Q)
print()

# print the velocity vector - left-handed coordinate system
print('Velocity x-coordinate (m/s):  ', P)
print('Velocity y-coordinate (m/s):  ', P)
print('Velocity z-coordinate (m/s):  ', P)
print()

This small script now does all the grunt work to perform the calculations and produces information that can be inserted directly into Orbiter.

example of use of scripts to construct NRHO orbit state vectors
When you run the above script, you should see the following information produced:

Code:
Moon-centric, ecliptic, Orbiter 2016 left-handed frame
-----------------------------------------------

MJD:    52013.754909351

Position x-coordinate (m):     -16115050.988265991
Position y-coordinate (m):     -3050794.458271753
Position z-coordinate (m):     -5524352.302307129

Velocity x-coordinate (m/s):   384.45187900043493
Velocity y-coordinate (m/s):   534.7812752345075
Velocity z-coordinate (m/s):   146.61220797426722

This is all the information that you need to insert a vessel on the prescribed point on the prescribed halo orbit at the prescribed time. The simplest way of doing this is to use the Scenario Editor:

1. Open a scenario and select the vessel that you want placed on the halo orbit.

2. Pause the Orbiter program and open the Scenario Editor.

3. Select the date option and paste in the MJD value given in the output and the hit 'Apply' 4. Select 'Done', then 'Edit', then 'State Vectors'.

5. Select the Moon as the 'orbit reference' (make sure that the frame is the 'ecliptic' and the coordinates are 'cartesian').

6. Now paste in the 6 position and velocity coordinate values. At the end of this, you should see something that looks like this: 7. Now hit 'Apply' and unpause Orbiter. Congratulations, your selected vessel is now in an NRHO orbit.

• francisdrake

#### ShootDrag

##### New member
Hey, hope you are doing well.

I was searching on Earth-Moon L2 Halo orbits and came across this thread.
Could you please direct me to the source, weblink or book where you learned math for the L2 halo Orbits?

Best,

#### MontBlanc2012

##### Member
Hey, hope you are doing well.

I was searching on Earth-Moon L2 Halo orbits and came across this thread.
Could you please direct me to the source, weblink or book where you learned math for the L2 halo Orbits?

Best,

Hi, ShootDrag

Not sure what your level of mathematical background is - but the basic scheme I used as a starting point was the paper:

"High-order solutions of invariant manifolds associated with libration point orbits in the elliptic restricted three-body system".

I've attached a copy of this for your reference.

This paper basically sets out a complicated scheme for calculating the Lindstedt–Poincaré (L-P) expansion of orbits around L1 and L2 Lagrange points. In principle, the scheme works for the Elliptical Restricted Three-Body Problem (ER3BP) - but the halo orbit family is really only well-defined for the Circular Restricted Three-Body Problem (CR3BP), which is a subset of the former (i.e., where one effectively assumes that the Moon orbits the Earth in a circular orbit). So, the expansions that I've used are based on that. You can use the L-P expansion scheme to generate accurate halo orbits 'close to' L1 or L2. For the near-rectilinear halo orbits, which are a 'long way' from L1 and L2, eventually one has to use more brute force numerical approaches. But if your interest is in halo orbits close to L2, this approach is good enough

Overall, this isn't a simple paper to work through - nor is it easy to apply in practice. I would think some background in Classical Mechanics at University level would probably be necessary to make sense of the paper.

I have developed some Mathematica code that solves the high-order L-P expansion problem for the super-set ER3BP (from which CR3BP high-order expansions are directly available), if that is of any use to you. I had plans to convert this code to a more friendly Python format but, alas, these ambitions were overrun by macro-economic chaos and COVID-19.

Regards

M.B.

#### Attachments

• High-order solutions of invariant manifolds associated with libration point orbits in the elli...pdf
1.9 MB · Views: 15

#### ShootDrag

##### New member
Hi, ShootDrag

Not sure what your level of mathematical background is - but the basic scheme I used as a starting point was the paper:

"High-order solutions of invariant manifolds associated with libration point orbits in the elliptic restricted three-body system".

I've attached a copy of this for your reference.

This paper basically sets out a complicated scheme for calculating the Lindstedt–Poincaré (L-P) expansion of orbits around L1 and L2 Lagrange points. In principle, the scheme works for the Elliptical Restricted Three-Body Problem (ER3BP) - but the halo orbit family is really only well-defined for the Circular Restricted Three-Body Problem (CR3BP), which is a subset of the former (i.e., where one effectively assumes that the Moon orbits the Earth in a circular orbit). So, the expansions that I've used are based on that. You can use the L-P expansion scheme to generate accurate halo orbits 'close to' L1 or L2. For the near-rectilinear halo orbits, which are a 'long way' from L1 and L2, eventually one has to use more brute force numerical approaches. But if your interest is in halo orbits close to L2, this approach is good enough

Overall, this isn't a simple paper to work through - nor is it easy to apply in practice. I would think some background in Classical Mechanics at University level would probably be necessary to make sense of the paper.

I have developed some Mathematica code that solves the high-order L-P expansion problem for the super-set ER3BP (from which CR3BP high-order expansions are directly available), if that is of any use to you. I had plans to convert this code to a more friendly Python format but, alas, these ambitions were overrun by macro-economic chaos and COVID-19.

Regards

M.B.

Thank you, MontBlanc2012.

I will go through the paper; I believe it shouldn't be that difficult.
Will ask if I encounter any problems. I hope that's okay.

Yeah, this Covid19 has been a real bummer. Hope you and your loved ones are healthy and doing good. Anyway, python would have been a little difficult as I mostly work in MATLAB, and seldom I use Mathematica. So yeah, Mathematica file, if available, should serve as a good understanding tool besides the paper itself.

Once again, thank you for taking out time to help.

Best,
S.D.

#### MontBlanc2012

##### Member
Thank you, MontBlanc2012.

I will go through the paper; I believe it shouldn't be that difficult.
Will ask if I encounter any problems. I hope that's okay.

Yeah, this Covid19 has been a real bummer. Hope you and your loved ones are healthy and doing good. Anyway, python would have been a little difficult as I mostly work in MATLAB, and seldom I use Mathematica. So yeah, Mathematica file, if available, should serve as a good understanding tool besides the paper itself.

Once again, thank you for taking out time to help.

Best,
S.D.

Hi, ShootDrag

I've attached my Mathematica notebook. This has been added with a '.txt' suffix rather that a '.nb' notebook suffix in order to 'fool' this forum's file uploading service. I suggest that after downloading, you change the file extension back to '.nb' from '.txt'.

This notebook employs a simplified version of the expansion scheme set out in the paper f, but just restricts the expansion to the centre manifold solutions (planar and vertical Lyapunov, Lissajous and Halo orbits). The stable and unstable manifold solutions can be examined in a similar fashion, but at the time I developed this notebook, I was only focusing on the centre manifold solutions so this notebook ignores other solutions.

You probably have access to Mathematica, but in case anyone else reading this thread doesn't, a Mathematica Player can be downloaded from https://www.wolfram.com/player/. This can be used to view the file, at least.

You should note that I've broadly checked that the algorithm works by doing a series of spot checks. But don't take it as guaranteed that the code necessarily works properly in all cases and all the time: this work was developed for my own interest - so, caveat emptor.

Apparently, '.nb' files

#### Attachments

• ER3BP - Lissajous - L1 & L2.txt
188.4 KB · Views: 14
• ShootDrag

#### ShootDrag

##### New member
Hey, MontBlanc.

I recently stumbled upon your youtube videos on similar topics, such as stable manifold trajectory & orbit insertion using stable manifold etc.
Since those videos only cover simulation, do you also happen to maintain a webpage or a GitHub repo containing algorithms and description for such topics?

Best,
S.D.

#### MontBlanc2012

##### Member
Hey, MontBlanc.

I recently stumbled upon your youtube videos on similar topics, such as stable manifold trajectory & orbit insertion using stable manifold etc.
Since those videos only cover simulation, do you also happen to maintain a webpage or a GitHub repo containing algorithms and description for such topics?

Best,
S.D.
Hi, ShoorDrag

Sorry, I don't maintain a webpage or github repository.

If I recall, I generated the stable manifold trajectory example using a time reversal trick: let's suppose you that you know the centre-manifold trajectory (e.g.,, by using the L-P expansion scheme as discussed above.) Take a point on a given centre-manifold trajectory and then perform a time-reverse integration from that initial point using a numerical integration scheme of your choice. Because the centre-manifold is unstable, the integration will eventually depart the centre-manifold along the (time-reversed) unstable manifold. But that is just a stable manifold trajectory in positive time. This method doesn't allow one to explore the full stable manifold associated with a given centre-manifold trajectory. However, it's good enough to provide an illustration of a stable manifold solution.

A more fulsome treatment requires expanding the L-P framework to include stable manifold solutions as well as the centre-manifold solutions. If you've made headway with the L-P process, this should be doable. This approach is far more controlled and allows a full exploration of the stable manifold solutions. Note also that there is no need to perform a high-order L-P expansion for the stable manifold.

Regards

M.B.

#### pantanplan

##### New member
Hello MontBlanc2012, I have been reading your detailed posts about halo orbits around libration points and finding them very useful, thanks for all the work you put in them. I would like to use the algorithm you describe in your "Numerical calculation of the Earth-Moon L2 halo orbit family" post, but the dropbox link to your Raw halo orbit interpolation data spreadsheet has been deleted. Do you by any chance still have that file and can post a link or send it to me?

#### MontBlanc2012

##### Member
Pantanplan: sorry for the slow response - but, of late, I only check in here at O-F once a month or so. I will have a search through my files and see if I can resurrect the spreadsheet.

Last edited:
• pantanplan

#### MontBlanc2012

##### Member
Pantanplan: OK, I’ve located the deleted dropbox file. I’ll send to you via PM.

P.S. If in the future any one else needs this file, let me know and I will send to you.

• pantanplan

#### Rock_Doc

##### New member
Hi, I came across your numerical calculation of the earth moon L2 halo orbit etc and would love to introduce my planetary class to this methodology, but your dropbox file is no longer there. Might it be possible to get that 7MB file from you? You did some nice work. thanks in advance

#### MontBlanc2012

##### Member
Hi, Rock_Doc

I think the attached is the missing file. Let me know if you have any issues.

Regards

Mont B

#### Attachments

• em l2 - halo - southern - 70.zip
6.5 MB · Views: 5

#### Rock_Doc

##### New member
Thanks MB. I am a bit overwhelmed by the intellect on this forum. I am a simple geologist, who has always focused on planetary geology, geomorphology, surface processes, etc. I will dive in to this. I appreciate this all.

#### Ajaja

##### Active member
That may be useful too: