MontBlanc2012
Active member
- Joined
- Aug 13, 2017
- Messages
- 138
- Reaction score
- 26
- Points
- 28
This post presents a tool for 'dialing-up' Orbiter-friendly coordinates of arbitrary centre-manifold orbits (Halo, Lissajous, planar and vertical Lyapunov) around L1 and L2 of the Earth-Moon system in the Elliptic Restricted Three-Body Problem (ER3BP).
Simply provide the MJD and the coordinates of the Earth and Moon as well as four parameters that describe which centre-manifold orbit that you want, and the tools return the Orbiter compatible state vectors that a vessel needs to have in order to be in the required centre-manifold orbit in the ER3BP theory.
This tool is based on a Lindstedt-Poincaré procedure for constructing a perturbation series solution of the equations of motion in the ER3BP model. In this case, the attached files are based on a 9th order series perturbation expansion yielding (for centre-manifold orbits less than circa 60,000 km in extent) high-fidelity solutions to those equations of motion. For
For those familiar with VSOP87 - the foundation of Orbiter's planetary ephemeris system - the tool works in much the same way in that it constructs the required state vectors from a set of data tables that encode a large number of cosine and sine terms.
In this post, I'm simply going to present the tool. In subsequent posts, I will give specific examples of how to use the tool to construct, for example, Halo orbits around EM L1 and EM L2. In subsequent work, I'll extend the datasets for this tool to the Sun-Earth system and even perhaps more exotic three-body systems - e.g., Sun-Mecury; Sun-Venus; Sun-Mars; Sun-Jupiter; and Jupiter-Galliean Moon.
Although this my seem a particularly 'dry' academic study, my guess that the datasets underpinning this tool lays the foundation for extending the range of Orbiter trajectory planing tools from largely Keplerian (two-body) basis to a more comprehensive three-body basis. As such, I'll probably also take the drastic step of posting the datasets together with supporting explanatory notes on Orbit Hangar.
The tool
The tool itself consists of an .xlsx based dataset used to house the sine and cosine tables of the Lindstedt-Poincaré solution; and a set of simple Python scripts to decode those tables and calculate the Orbiter compatible state vectors of the required centre-manifold orbit. The relevant datasets for EM L1 and EM L2 are in the attached .zip file which contains two files:
'EM L1 - Lissajous - ER3BP.xlsx'
'EM L2 - Lissajous - ER3BP.xlsx'
Not surprisingly, the first is the dataset that describes centre-manifold orbits about EM L1; and the second is the dataset that describes central-manifold orbits about EM L2. The reference to Lissajous in the title is a reminder that the encoded solutions are general Lissajous centre-manifold orbits. Specific cases of the general Lissajous orbits are: the planar and vertical Lyapunov orbits and Halo orbits.
To decode these datasets, one needs to load the datasets into a suitable computational engine; perform some simple operations on the data; and have the computational engine return the Orbiter-friendly state vectors. For (my) convenience, I've elected to use Python as the computational engine of choice: it's free, it's widely available; and it is supported by a wide range of robust mathematical packages. It would be possible to use a C/C++ based computational engine, but I have considerably less familiarity with that language so for the time being, I'll leave it to others to transcribe the Python tools into similar C/C++ ones. The relevant Python scripts are presented below - but these scripts assume that you have downloaded the datasets contained in the attached .zip file to some known location on your computer.
The first step in using the Python tool is to extract the dataset from the .xlsx spreadsheet and load it into Python. This step only needs to be done once. For the EM L2 dataset, we have:
naturally, you have to set the path to file for this script to work but it will simply be the path to wherever you stored 'EM L2 - Lissajous - ER3BP.xlsx' on your machine.
This script uses the 'openpyxl', 'numpy' and 'math' packages. The 'math' package is a default package available in most Python installations, but you will need to separately install the 'openpyxl' and 'numpy' packages. The former contains utilities for extracting data from an .xlsx file; and the latter contains some general purpose utilities for working with data lists. Basically, what the script does is extract the six worksheets from the .xlsx spreadsheet and store them as Python variables - so all this utility does is move data from the Excel dataset into Python.
Having loaded the dataset into Python, we now need some utilities for decoding the data so that we can construct relevant state vectors. There are three utility functions that we need to do this:
This Python code snippet defines four functions:
1. 'getOmega' and 'getLambda' return the two fundamental frequencies of the ER3BP - one (omega) in the x-y plane; and the other (lambda) in the vertical z direction. These frequencies are functions of three parameters - 'e', 'da1' and 'da2'. The first 'e' is the orbital eccentricity of the Moon's elliptical orbit around the Earth. As such, for any given epoch, it is a 'given' quantity. The other two parameters, 'da1' and 'da2' are user-defied scale parameters. The first controls the size of the Lissajous orbit in the x-y plane; and the second controls the size of the Lissajous orbit in the z-direction.
2. 'pos' calculates the state vectors of the required Lissajous orbit in the non-inertial coordinate system of the ER3BP theory. This coordinate system includes a scale transformation; a rotation; and a change of time variable. Consequently, the state vector that is returned isn't immediately useful to the Orbiter users. This function has six parameters - 'e', 'da1', 'da2', 'nu', 'phi1' and 'phi2'. These six parameters fully describe the Lissajous orbit and play the same role as the six orbital parameters of a standard Keplerian orbit. The first three are the same as described above but we now have three additional parameters - 'nu', 'phi1' and 'phi2'. 'nu' (in radians) is the true anomaly of the Moon in its elliptical orbit around the Earth; and 'phi1' and 'phi2' are user-defined phase angles for motion in the x-y plane and in the z-direction respectively. In total, then, there are just four user-defined quantities to describe a halo orbit: 'da1' and 'phi1' for motion in the x-y plane; and 'da2' and 'phi2' for motion in the z-direction.
3. 'getECIcoords' calculates the Orbiter-friendly state-vectors of the vessel in the relevant ECI reference frame calculated using the 'pos' function. The state-vectors can be used directly in Orbiter using the Scenario Editor tool. This function has 4 inputs - 'a', 'e', 'nu' and 'pos'. 'a' is the semi-major axis of the Moon's orbit around the Earth; 'e' is the orbital eccentricity of the Moon; and 'nu' is the true anomaly of the Moon.
Having now loaded the dataset and defined some utility functions, we can now do some computations. To get the ball rolling, we need to provide some basic information about the Moon's orbit and where it is in that orbit - i.e., 'a', 'e' and 'nu'. This information can either be calculated either by interrogating Orbiter's ephemeris system or from some other source (if available). We then need to provide the four user-defined parameters - 'da1', 'da2', 'phi1' and 'phi2'. Then we call the 'pos' utility function followed by the getECIcoords utility function. Having done this, we can then immediately read off the Orbiter-Friendly state vectors of the vessel in the ECI reference frame.
An example of this sequence is given below:
(For those interested, this example calculates the state vectors of a vessel in a near halo orbit about EM L2.)
As a follow on to this post, I'm going to provide some illustrations on how to use this tool - first to calculate planar Lyapunov orbits, then vertical Lyapunov orbits; general Lissajous orbits and, finally, halo orbits. I'll also work on a more formal posting on Orbit Hangar. Having done all of that, I'll then move onto using these tools to provide reference orbits so that (with a little bit of station-keeping) it should be possible to keep a vessel on a prescribed centre-manifold orbit indefinitely.
Simply provide the MJD and the coordinates of the Earth and Moon as well as four parameters that describe which centre-manifold orbit that you want, and the tools return the Orbiter compatible state vectors that a vessel needs to have in order to be in the required centre-manifold orbit in the ER3BP theory.
This tool is based on a Lindstedt-Poincaré procedure for constructing a perturbation series solution of the equations of motion in the ER3BP model. In this case, the attached files are based on a 9th order series perturbation expansion yielding (for centre-manifold orbits less than circa 60,000 km in extent) high-fidelity solutions to those equations of motion. For
For those familiar with VSOP87 - the foundation of Orbiter's planetary ephemeris system - the tool works in much the same way in that it constructs the required state vectors from a set of data tables that encode a large number of cosine and sine terms.
In this post, I'm simply going to present the tool. In subsequent posts, I will give specific examples of how to use the tool to construct, for example, Halo orbits around EM L1 and EM L2. In subsequent work, I'll extend the datasets for this tool to the Sun-Earth system and even perhaps more exotic three-body systems - e.g., Sun-Mecury; Sun-Venus; Sun-Mars; Sun-Jupiter; and Jupiter-Galliean Moon.
Although this my seem a particularly 'dry' academic study, my guess that the datasets underpinning this tool lays the foundation for extending the range of Orbiter trajectory planing tools from largely Keplerian (two-body) basis to a more comprehensive three-body basis. As such, I'll probably also take the drastic step of posting the datasets together with supporting explanatory notes on Orbit Hangar.
The tool
The tool itself consists of an .xlsx based dataset used to house the sine and cosine tables of the Lindstedt-Poincaré solution; and a set of simple Python scripts to decode those tables and calculate the Orbiter compatible state vectors of the required centre-manifold orbit. The relevant datasets for EM L1 and EM L2 are in the attached .zip file which contains two files:
'EM L1 - Lissajous - ER3BP.xlsx'
'EM L2 - Lissajous - ER3BP.xlsx'
Not surprisingly, the first is the dataset that describes centre-manifold orbits about EM L1; and the second is the dataset that describes central-manifold orbits about EM L2. The reference to Lissajous in the title is a reminder that the encoded solutions are general Lissajous centre-manifold orbits. Specific cases of the general Lissajous orbits are: the planar and vertical Lyapunov orbits and Halo orbits.
To decode these datasets, one needs to load the datasets into a suitable computational engine; perform some simple operations on the data; and have the computational engine return the Orbiter-friendly state vectors. For (my) convenience, I've elected to use Python as the computational engine of choice: it's free, it's widely available; and it is supported by a wide range of robust mathematical packages. It would be possible to use a C/C++ based computational engine, but I have considerably less familiarity with that language so for the time being, I'll leave it to others to transcribe the Python tools into similar C/C++ ones. The relevant Python scripts are presented below - but these scripts assume that you have downloaded the datasets contained in the attached .zip file to some known location on your computer.
The first step in using the Python tool is to extract the dataset from the .xlsx spreadsheet and load it into Python. This step only needs to be done once. For the EM L2 dataset, we have:
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].valuenaturally, you have to set the path to file for this script to work but it will simply be the path to wherever you stored 'EM L2 - Lissajous - ER3BP.xlsx' on your machine.
This script uses the 'openpyxl', 'numpy' and 'math' packages. The 'math' package is a default package available in most Python installations, but you will need to separately install the 'openpyxl' and 'numpy' packages. The former contains utilities for extracting data from an .xlsx file; and the latter contains some general purpose utilities for working with data lists. Basically, what the script does is extract the six worksheets from the .xlsx spreadsheet and store them as Python variables - so all this utility does is move data from the Excel dataset into Python.
Having loaded the dataset into Python, we now need some utilities for decoding the data so that we can construct relevant state vectors. There are three utility functions that we need to do this:
Code:
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]]
def getECIcoords(a, e, nu, pos):
[[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(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]]This Python code snippet defines four functions:
1. 'getOmega' and 'getLambda' return the two fundamental frequencies of the ER3BP - one (omega) in the x-y plane; and the other (lambda) in the vertical z direction. These frequencies are functions of three parameters - 'e', 'da1' and 'da2'. The first 'e' is the orbital eccentricity of the Moon's elliptical orbit around the Earth. As such, for any given epoch, it is a 'given' quantity. The other two parameters, 'da1' and 'da2' are user-defied scale parameters. The first controls the size of the Lissajous orbit in the x-y plane; and the second controls the size of the Lissajous orbit in the z-direction.
2. 'pos' calculates the state vectors of the required Lissajous orbit in the non-inertial coordinate system of the ER3BP theory. This coordinate system includes a scale transformation; a rotation; and a change of time variable. Consequently, the state vector that is returned isn't immediately useful to the Orbiter users. This function has six parameters - 'e', 'da1', 'da2', 'nu', 'phi1' and 'phi2'. These six parameters fully describe the Lissajous orbit and play the same role as the six orbital parameters of a standard Keplerian orbit. The first three are the same as described above but we now have three additional parameters - 'nu', 'phi1' and 'phi2'. 'nu' (in radians) is the true anomaly of the Moon in its elliptical orbit around the Earth; and 'phi1' and 'phi2' are user-defined phase angles for motion in the x-y plane and in the z-direction respectively. In total, then, there are just four user-defined quantities to describe a halo orbit: 'da1' and 'phi1' for motion in the x-y plane; and 'da2' and 'phi2' for motion in the z-direction.
3. 'getECIcoords' calculates the Orbiter-friendly state-vectors of the vessel in the relevant ECI reference frame calculated using the 'pos' function. The state-vectors can be used directly in Orbiter using the Scenario Editor tool. This function has 4 inputs - 'a', 'e', 'nu' and 'pos'. 'a' is the semi-major axis of the Moon's orbit around the Earth; 'e' is the orbital eccentricity of the Moon; and 'nu' is the true anomaly of the Moon.
Having now loaded the dataset and defined some utility functions, we can now do some computations. To get the ball rolling, we need to provide some basic information about the Moon's orbit and where it is in that orbit - i.e., 'a', 'e' and 'nu'. This information can either be calculated either by interrogating Orbiter's ephemeris system or from some other source (if available). We then need to provide the four user-defined parameters - 'da1', 'da2', 'phi1' and 'phi2'. Then we call the 'pos' utility function followed by the getECIcoords utility function. Having done this, we can then immediately read off the Orbiter-Friendly state vectors of the vessel in the ECI reference frame.
An example of this sequence is given below:
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])
# calcuale a, e, and nu
GE = 3.986004418e14
GM = 4.9048695e12
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
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 - right-handed coordinate system
p = pos(e, 0.036, 0.02144854785901588, nu, 0.0, 0.0)
[[Qx, Qy, Qz], [Px, Py, Pz]] = getECIcoords(a, e, nu, p)
Q = Qx * xhat + Qy * yhat + Qz * zhat + com - Qe
P = Px * xhat + Py * yhat + Pz * zhat + cov - Pe
# print the position vector - left-handed coordinate system
print(Q[0])
print(Q[2])
print(Q[1])
# print the velocity vector - left-handed coordinate system
print(P[0])
print(P[2])
print(P[1])(For those interested, this example calculates the state vectors of a vessel in a near halo orbit about EM L2.)
As a follow on to this post, I'm going to provide some illustrations on how to use this tool - first to calculate planar Lyapunov orbits, then vertical Lyapunov orbits; general Lissajous orbits and, finally, halo orbits. I'll also work on a more formal posting on Orbit Hangar. Having done all of that, I'll then move onto using these tools to provide reference orbits so that (with a little bit of station-keeping) it should be possible to keep a vessel on a prescribed centre-manifold orbit indefinitely.
Last edited: