MontBlanc2012
Active member
- Joined
- Aug 13, 2017
- Messages
- 138
- Reaction score
- 26
- Points
- 28
This is a follow-on post from Dialing-up arbitrary Lagrange Point orbits. In it, I'm going to demonstrate how to place a vessel in a planar Lyapunov orbit around EM L2 using the tools described in that earlier post. (If you want to see an illustration of a planar Lyapunov orbit, see here: More on Lagrange point orbits)
The first thing we need to decide is 'when' we want to place a vessel in a planar Lyapunov orbit - i.e., we need to choose a starting MJD for our Lyapunov orbit. For convenience, I'm going to choose the MJD of the example given in Dialing-up arbitrary Lagrange Point orbits, namely MJD = 52013.754909351. There is nothing particularly special about this date - it's just the date that I happened to extract from Orbiter 2016 when I was putting together the example for my earlier post.
The next thing that we need to do is to establish some of the orbital parameters of the Moon's motion around the Earth. In particular, we need to determine the semi-major axis; the orbital eccentricity and the true anomaly. We need this information because calculating a Lyapunov orbit is done in the context of a restricted three-body theory - the Earth, the Moon and the vessel - so we need to say something about the Moon's orbit around the Earth. The semi-major axis; the orbital eccentricity and the true anomaly captures the relevant information. So, where do we get this information from? One way is to open up Orbit MFD and just copy down the relevant numbers when referencing the Earth and targeting the Moon. However, this approach doesn't provide much precision. So, my preferred approach is to interrogate the Orbiter 2016 ephemeris system directly using a Lua script and calculate the state vectors (position and velocity) of the Earth and Moon at the epoch that we have chosen and then calculate the Moon's orbital parameters directly.
Interrogating Orbiter's ephemeris system is simple using the assorted api constructs provided with Orbiter. In terms of a Lua script, one would run something like:
If we use a script like this, then we should find that the script will write the current state vectors of the Earth and Moon to the fine 'statevectors.txt'. These state-vectors will be expressed in a left-handed coordinate system. To convert to a more conventional right-handed coordinate system, we just swap the 'y' and 'z' components of each vector. And when we have finished doing conventional calculations in a right-handed coordinate system and want to ship vectors back into Orbiter, we need to remember to swap the 'y' and 'z' components back again to an Orbiter-friendly format.
Anyway, if we run a Lua script like this for MJD = 52013.754909351, we should find that we get something like:
Putting all this together in Python, we would write:
Qe and Pe hold the position and velocity vectors of the Earth in a right-handed coordinate system in Orbiter's global (inertial) reference frame for the given MJD. Similarly, Qm and Pm hold the position and velocity vectors of the Moon in the same reference frame at the same epoch.
From this information, we calculate in Python the required orbital parameters as follows:
We should find that e = 0.06527555911340753; a = 380105578.6068151 metres; and nu = 2.575177554472375 radians.
When we eventually calculate the state vectors of the starting point of our Lyapunov orbit, we also need to specify the orientation of the Moon's orbit. This is given by calculating the following unit vectors:
Together these three vectors - xhat, yhat and zhat - encode the three Euler angles (inclination, Longitude of Ascending Node and Argument of Periapsis) conventionally used to encode the orientation of an orbital plane.
OK, so now we have all the relevant information about the (current) orbit of the Moon around the Earth. Next we use this information to calculate the state vectors of the Lyapunov orbit. From the earlier post, one needs to specify six parameters - 'e', 'da1', 'da2', 'nu', 'phi1' and 'phi2'. Now, we already know 'e' and 'nu' from the preceding analysis of the Moon's orbit. So, e = 0.06527555911340753; and nu = 2.575177554472375. What about the other parameters? Now, in this example, we want to focus on planar Lyapunov orbits. These are defined such that 'da2' = 0. Since we have 'da2' to zero, it doesn't matter what value we use for 'phi2' and so, for convenience, we will set it to 0.0.
This just leaves us with two free parameters, 'da1' and 'phi1'. The latter is a phase angle which we are able to choose freely between 0 and 2*pi. Here, we'll choose 'phi1' = 0.0 if for no other reason than it convenient. The remaining parameter 'da1'>0 is a scale parameter which sets the size of the planar Lyapunov orbit. Here we have to be a little careful in choosing the value of 'da1'. Because the underlying mathematical theory is based on a Lindstedt-Poincaré perturbation analysis, if we choose too large a value, the perturbative solution that is encoded in the Lagrange Point datasets will not provide a good approximation to the true solution of the equations of motion of the elliptical restricted tree-body problem (ER3BP). So, how big is too big? Well, if we choose a value of 'da1' in the range 0.00 to about 0.04 we should achieve satisfactory results - but for values of 'da1' less than 0.02, the underlying Lindstedt-Poincaré perturbation analysis is highly accurate. But how big is an orbit corresponding to 'da1'=0.02? Well, we can roughly estimate the size of the Lyapunov orbit by calculating: da1 * 3.0* 2 * 385000 km. This is the end-to-end length of the orbit. For da1 = 0.02, we estimate that the end-to-end-length of the orbit is around 46,000 km. This is pretty large - and with an upper bound of da1 = 0.04, the end-to-end length of the orbit is over 90,000 km. Frankly, I can't imagine circumstances where one would need to work with Lyapunov orbits larger than this.
Here, we'll work with da1 = 0.02. Using the Python script introduced in the previous post, we run the following script:
and this should yield the following Orbiter-friendly state-vectors:
The position components - x, y and z - are given in metres and are Earth-centric inertial (ECI); and the velocity components are given in m/s and also Earth-centric inertial (ECI). This values can be cut and paste directly into Orbiter using the Scenario Editor tool. First load Orbiter and start a scenario running. Pause Orbiter. Then open the Scenario Editor and set the date:
After setting the date, cut and paste the state vectors in the State Vectors dialog (remembering to set the Orbit reference to the Earth; the Frame to the Ecliptic; using Cartesian coordinates):
Congratulations! You have just placed a vessel in a prescribe planar Lyapunov orbit.
IMPORTANT NOTE: Although the Lindstedt-Poincare perturbation analysis gives high-precision to the elliptical restricted three-body problem, the model assumes that the Moon's orbit is not perturbed. In reality, it is quite heavily perturbed by the Sun, and so we do not expect the vessel to stay on the reference Lyapunov orbit for very long without station-keeping. Later on, I'm going to show how you can overcome these perturbative effects and keep a vessel on the reference Lyapunov orbit indefinitely.
(And in case anyone was wondering what the planar Lyapunov orbit looks like, (in the somewhat obscure rotating coordinates of the ER3BP theory) here it is (for three 3 months beyond the epoch MJD at least):
Note that the Lyapunov orbit is not a closed loop. If the Moon's orbit were circular, it would be - but the Moon's orbit is quite highly eccentric, and that eccentricity knocks the centre-manifold around a bit so that it no longer forms closed loops. The result is a quasi-periodic structure. And you should also note that if we chose a different value of 'phi1' our orbit would look different - basically the same size and shape, but still, a different quasi-periodic orbit.
The first thing we need to decide is 'when' we want to place a vessel in a planar Lyapunov orbit - i.e., we need to choose a starting MJD for our Lyapunov orbit. For convenience, I'm going to choose the MJD of the example given in Dialing-up arbitrary Lagrange Point orbits, namely MJD = 52013.754909351. There is nothing particularly special about this date - it's just the date that I happened to extract from Orbiter 2016 when I was putting together the example for my earlier post.
The next thing that we need to do is to establish some of the orbital parameters of the Moon's motion around the Earth. In particular, we need to determine the semi-major axis; the orbital eccentricity and the true anomaly. We need this information because calculating a Lyapunov orbit is done in the context of a restricted three-body theory - the Earth, the Moon and the vessel - so we need to say something about the Moon's orbit around the Earth. The semi-major axis; the orbital eccentricity and the true anomaly captures the relevant information. So, where do we get this information from? One way is to open up Orbit MFD and just copy down the relevant numbers when referencing the Earth and targeting the Moon. However, this approach doesn't provide much precision. So, my preferred approach is to interrogate the Orbiter 2016 ephemeris system directly using a Lua script and calculate the state vectors (position and velocity) of the Earth and Moon at the epoch that we have chosen and then calculate the Moon's orbital parameters directly.
Interrogating Orbiter's ephemeris system is simple using the assorted api constructs provided with Orbiter. In terms of a Lua script, one would run something like:
Code:
earth = oapi.get_objhandle("Earth")
moon = oapi.get_objhandle("Moon")
-- get the current location of Earth
q_ear = oapi.get_globalpos(earth)
p_ear = oapi.get_globalvel(earth)
-- get the current location of the Moon
q_mon = oapi.get_globalpos(moon)
p_mon = oapi.get_globalvel(moon)
file = io.open("statevectors.txt", "w")
io.output(file)
file:write("qe = \{", q_ear.x, ", ", q_ear.y, ", ", q_ear.z, "\}\n")
file:write("pe = \{", p_ear.x, ", ", p_ear.y, ", ", p_ear.z, "\}\n")
file:write("qm = \{", q_mon.x, ", ", q_mon.y, ", ", q_mon.z, "\}\n")
file:write("pm = \{", p_mon.x, ", ", p_mon.y, ", ", p_mon.z, "\}\n")
file:close()If we use a script like this, then we should find that the script will write the current state vectors of the Earth and Moon to the fine 'statevectors.txt'. These state-vectors will be expressed in a left-handed coordinate system. To convert to a more conventional right-handed coordinate system, we just swap the 'y' and 'z' components of each vector. And when we have finished doing conventional calculations in a right-handed coordinate system and want to ship vectors back into Orbiter, we need to remember to swap the 'y' and 'z' components back again to an Orbiter-friendly format.
Anyway, if we run a Lua script like this for MJD = 52013.754909351, we should find that we get something like:
Code:
qe = {-136757075937.97, 20730124.687839, -63826187339.785}
pe = {12032.791283836, 0.8388150790952, -27150.845638622}
qm = {-136653091584.32, 17225090.561136, -64212987439.536}
pm = {12980.539253522, -85.194965428945, -26932.733069022}Putting all this together in Python, we would write:
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])Qe and Pe hold the position and velocity vectors of the Earth in a right-handed coordinate system in Orbiter's global (inertial) reference frame for the given MJD. Similarly, Qm and Pm hold the position and velocity vectors of the Moon in the same reference frame at the same epoch.
From this information, we calculate in Python the required orbital parameters as follows:
Code:
# 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 - nuWe should find that e = 0.06527555911340753; a = 380105578.6068151 metres; and nu = 2.575177554472375 radians.
When we eventually calculate the state vectors of the starting point of our Lyapunov orbit, we also need to specify the orientation of the Moon's orbit. This is given by calculating the following unit vectors:
Code:
xhat= evec / e
zhat= np.cross(r,v)
zhat= zhat / np.linalg.norm(zhat)
yhat= np.cross(zhat,xhat)Together these three vectors - xhat, yhat and zhat - encode the three Euler angles (inclination, Longitude of Ascending Node and Argument of Periapsis) conventionally used to encode the orientation of an orbital plane.
OK, so now we have all the relevant information about the (current) orbit of the Moon around the Earth. Next we use this information to calculate the state vectors of the Lyapunov orbit. From the earlier post, one needs to specify six parameters - 'e', 'da1', 'da2', 'nu', 'phi1' and 'phi2'. Now, we already know 'e' and 'nu' from the preceding analysis of the Moon's orbit. So, e = 0.06527555911340753; and nu = 2.575177554472375. What about the other parameters? Now, in this example, we want to focus on planar Lyapunov orbits. These are defined such that 'da2' = 0. Since we have 'da2' to zero, it doesn't matter what value we use for 'phi2' and so, for convenience, we will set it to 0.0.
This just leaves us with two free parameters, 'da1' and 'phi1'. The latter is a phase angle which we are able to choose freely between 0 and 2*pi. Here, we'll choose 'phi1' = 0.0 if for no other reason than it convenient. The remaining parameter 'da1'>0 is a scale parameter which sets the size of the planar Lyapunov orbit. Here we have to be a little careful in choosing the value of 'da1'. Because the underlying mathematical theory is based on a Lindstedt-Poincaré perturbation analysis, if we choose too large a value, the perturbative solution that is encoded in the Lagrange Point datasets will not provide a good approximation to the true solution of the equations of motion of the elliptical restricted tree-body problem (ER3BP). So, how big is too big? Well, if we choose a value of 'da1' in the range 0.00 to about 0.04 we should achieve satisfactory results - but for values of 'da1' less than 0.02, the underlying Lindstedt-Poincaré perturbation analysis is highly accurate. But how big is an orbit corresponding to 'da1'=0.02? Well, we can roughly estimate the size of the Lyapunov orbit by calculating: da1 * 3.0* 2 * 385000 km. This is the end-to-end length of the orbit. For da1 = 0.02, we estimate that the end-to-end-length of the orbit is around 46,000 km. This is pretty large - and with an upper bound of da1 = 0.04, the end-to-end length of the orbit is over 90,000 km. Frankly, I can't imagine circumstances where one would need to work with Lyapunov orbits larger than this.
Here, we'll work with da1 = 0.02. Using the Python script introduced in the previous post, we run the following script:
Code:
# now calculate the ECI state vectors - right-handed coordinate system
p = pos(e, 0.02, 0.0, 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('x = ', Q[0])
print('y = ', Q[2])
print('z = ', Q[1])
print()
# print the velocity vector - left-handed coordinate system
print('dx/dt = ', P[0])
print('dy/dt = ', P[2])
print('dz/dt = ', P[1])and this should yield the following Orbiter-friendly state-vectors:
Code:
x = 142630732.39672852
y = -6060181.297925532
z = -443869440.7295227
dx/dt = 1091.2763792888745
dy/dt = -99.33037705787646
dz/dt = 269.6587591310308The position components - x, y and z - are given in metres and are Earth-centric inertial (ECI); and the velocity components are given in m/s and also Earth-centric inertial (ECI). This values can be cut and paste directly into Orbiter using the Scenario Editor tool. First load Orbiter and start a scenario running. Pause Orbiter. Then open the Scenario Editor and set the date:
After setting the date, cut and paste the state vectors in the State Vectors dialog (remembering to set the Orbit reference to the Earth; the Frame to the Ecliptic; using Cartesian coordinates):
Congratulations! You have just placed a vessel in a prescribe planar Lyapunov orbit.
IMPORTANT NOTE: Although the Lindstedt-Poincare perturbation analysis gives high-precision to the elliptical restricted three-body problem, the model assumes that the Moon's orbit is not perturbed. In reality, it is quite heavily perturbed by the Sun, and so we do not expect the vessel to stay on the reference Lyapunov orbit for very long without station-keeping. Later on, I'm going to show how you can overcome these perturbative effects and keep a vessel on the reference Lyapunov orbit indefinitely.
(And in case anyone was wondering what the planar Lyapunov orbit looks like, (in the somewhat obscure rotating coordinates of the ER3BP theory) here it is (for three 3 months beyond the epoch MJD at least):
Note that the Lyapunov orbit is not a closed loop. If the Moon's orbit were circular, it would be - but the Moon's orbit is quite highly eccentric, and that eccentricity knocks the centre-manifold around a bit so that it no longer forms closed loops. The result is a quasi-periodic structure. And you should also note that if we chose a different value of 'phi1' our orbit would look different - basically the same size and shape, but still, a different quasi-periodic orbit.
Last edited:
