Orbiter-Forum Numerical calculation of the Earth-Moon L2 halo orbit family
 Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

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

 03-05-2019, 08:25 AM #1 MontBlanc2012 Orbinaut 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)
 Thanked by:
 03-07-2019, 02:11 AM #2 MontBlanc2012 Orbinaut Well, in the preceding post in this thread I wrote: Quote: 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))[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'``` 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[0]) print('Position y-coordinate (m): ', Q[2]) print('Position z-coordinate (m): ', Q[1]) print() # 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]) 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.
 Thanked by:

 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 User Control Panel Private Messages Subscriptions Who's Online Search Forums Forums Home Orbiter-Forum.com     Announcements     Meets & Greets Orbiter Space Flight Simulator     Orbiter Web Forum         OFMM         Orbiter Forum Space Station         Simpit Forum     General Questions & Help     MFD Questions & Help     Hardware & Software Help     Tutorials & Challenges     Orbiter SDK     Orbiter Visualization Project     Orbiter Beta » Orbiter Project Orbiter Addons     OrbitHangar Addons & Comments     Addons     Addon Development     Addon Requests     Addon Support & Bugs         Addon Developer Forums             Project Apollo - NASSP     Orbiter Lua Scripting Far Side of the Moon     Spaceflight News     Math & Physics     Astronomy & the Night Sky     Backyard Rocketry     Brighton Lounge     International Forum

All times are GMT. The time now is 01:40 AM.