Go Back   Orbiter-Forum > Far Side of the Moon > Math & Physics
Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

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

Thread Tools
Old 03-05-2019, 08:25 AM   #1
Default Numerical calculation of the Earth-Moon L2 halo orbit family

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 is offline   Reply With Quote
Thanked by:
Old 03-07-2019, 02:11 AM   #2

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:

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:

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))[0]
itp  = np.transpose(np.array((xl.parse('Interp. pts')).values))[0]

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

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]
        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:

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.

# 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:

# 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:

# 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:

# 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('MJD:   ', mjd)
print('Position x-coordinate (m):    ', Q[0])
print('Position y-coordinate (m):    ', Q[2])
print('Position z-coordinate (m):    ', Q[1])

# print the velocity vector - left-handed coordinate system
print('Velocity x-coordinate (m/s):  ', P[0])
print('Velocity y-coordinate (m/s):  ', P[2])
print('Velocity z-coordinate (m/s):  ', P[1])
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:

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.
MontBlanc2012 is offline   Reply With Quote
Thanked by:

  Orbiter-Forum > Far Side of the Moon > Math & Physics

Thread Tools

Posting Rules
BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
Forum Jump

All times are GMT. The time now is 12:53 PM.

Quick Links Need Help?

About Us | Rules & Guidelines | TOS Policy | Privacy Policy

Orbiter-Forum is hosted at Orbithangar.com
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2019, vBulletin Solutions Inc.
Copyright ©2007 - 2017, Orbiter-Forum.com. All rights reserved.