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

Thread Tools 
03052019, 08:25 AM  #1 
Orbinaut

Numerical calculation of the EarthMoon 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 centremanifold orbits (planar and vertical Lyapunov, Lissajous and halo orbits) of the L1 and L2 Lagrange points  particularly those of the EarthMoon and the SunEarth systems.
So far I've posted only maps using a form of perturbation theory (known as LindstedtPoincaré 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 Excelbased 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 LOPG (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 LOPGlike 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) 
Thanked by: 
03072019, 02:11 AM  #2 
Orbinaut

Well, in the preceding post in this thread I wrote:
Quote:
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 GaussLobatto style variational integrator used to construct the 1300 or so individual halo orbits. On the first six sheet, the GaussLobatto 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 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))[0] itp = np.transpose(np.array((xl.parse('Interp. pts')).values))[0] 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' Two inputs into this function are required because the halo orbit manifold is a curved twodimensional structure in threedimensional 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 /(1e*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]] 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  righthanded 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  righthanded 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 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 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 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  righthanded 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  lefthanded coordinate system print('Mooncentric, ecliptic, Orbiter 2016 lefthanded frame') print('') print() print('MJD: ', mjd) print() print('Position xcoordinate (m): ', Q[0]) print('Position ycoordinate (m): ', Q[2]) print('Position zcoordinate (m): ', Q[1]) print() # print the velocity vector  lefthanded coordinate system print('Velocity xcoordinate (m/s): ', P[0]) print('Velocity ycoordinate (m/s): ', P[2]) print('Velocity zcoordinate (m/s): ', P[1]) print() 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:
Mooncentric, ecliptic, Orbiter 2016 lefthanded frame  MJD: 52013.754909351 Position xcoordinate (m): 16115050.988265991 Position ycoordinate (m): 3050794.458271753 Position zcoordinate (m): 5524352.302307129 Velocity xcoordinate (m/s): 384.45187900043493 Velocity ycoordinate (m/s): 534.7812752345075 Velocity zcoordinate (m/s): 146.61220797426722 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. 
Thanked by: 

Thread Tools  


Quick Links  Need Help? 