MontBlanc2012
Active member
- Joined
- Aug 13, 2017
- Messages
- 138
- Reaction score
- 26
- Points
- 28
This is a second follow-on post from Dialing-up arbitrary Lagrange Point orbits. In it, I'm going to demonstrate how to place a vessel in a vertical Lyapunov orbit around EM L2 using the tools described in the earlier post. The first follow-on post Dialing-up a planar Lyapunov orbit demonstrated how to place a vessel in a planar Lyapunov orbit. For a discussion on Lyapunov orbits see "More on Lagrange point orbits".
(If you want to see an illustration of a vertical Lyapunov orbit, see the following short vide clip which illustrates vertical Lyapunov orbits and contrasts them with planar Lyapunov orbits -
.)
As with planar Lyapunov orbits, the first thing we need to decide is 'when' we want to place a vessel in a vertical Lyapunov orbit - i.e., we need to choose a starting MJD for our orbit. Again for convenience, I'm going to choose the MJD of the example given in Dialing-up arbitrary Lagrange Point orbits, namely MJD = 52013.754909351. Again, there is nothing particularly special about this date - it just happens to be the date that I happened to extract from Orbiter 2016 when I was putting together the example for my earlier post. But any MJD will work just as well.
And, again, the next thing that we need to do is to establish some of the orbital parameters of the Moon's motion around the Earth on the afore mentioned MJD that we have just determined. In particular, we need to determine the semi-major axis; the orbital eccentricity and the true anomaly. As before, for the given MJD, we use the following Lua script to interrogate Orbiter's ephemeris system ad extract the relevant state vectors of the Earth and the Moon:
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 file 'statevectors.txt'. These state-vectors will be 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 into their Orbiter-friendly form. 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 finally 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 vertical Lyapunov orbits. These are defined such that 'da1' = 0.0. And with 'da1' = 0.0 we can set 'phi1' to zero for our convenience (although it doesn't really matter if we take any value between 0 and 2*pi).
This just leaves us with two free parameters, 'da2' and 'phi2', with which to define an vertical Lyapunov orbit. The latter is a phase angle which we are able to choose freely between 0 and 2*pi. In keeping with the approach adopted in Dialing-up a planar Lyapunov orbit, we'll choose phi2 = 0.0 if for no other reason than it is convenient to do so. The remaining parameter 'da2'>0 is a scale parameter which sets the size of the vertical Lyapunov orbit. As before, we have to be a little careful in our choice of value here. Because the underlying mathematical theory is based on a Lindstedt-Poincaré perturbation analysis, if we choose too large a value, the perturbation solution will not be a good approximation to the true solution of the equations of motion. So, how big is too big? Well, if we choose a value of da2 in the range 0.00 to about 0.08 we should achieve satisfactory results - but for values of 'da2' less than 0.04 or so, the underlying Lindstedt-Poincaré perturbation analysis is highly accurate. But how big is a vertical Lyapunov orbit with 'da2'=0.02? Well, we can roughly estimate the size of the Lyapunov orbit by calculating: da2 * 1.00 * 2 * 385000 km. This is the end-to-end length of the orbit. For da2 = 0.02, we estimate that the end-to-end-length of the orbit is around 15,400 km. With an upper bound of 'da2' around 0.10, the end-to-end length of the vertical Lyapunov orbit that we can reasonably approximate using our perturbative theory is around 75,000 km.
Here, we'll work with da2 = 0.02. Using the Python script introduced in the previous post, we run the following script:
and this should yield the following Orbiter-friendly ECI 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:
Then 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 prescribed vertical Lyapunov orbit.
IMPORTANT NOTE: Although the Lindstedt-Poincare perturbation analysis gives high-precision solution 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 some form of 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 you were wondering, here is what the prescribed Lyapunov orbit looks like:
(If you want to see an illustration of a vertical Lyapunov orbit, see the following short vide clip which illustrates vertical Lyapunov orbits and contrasts them with planar Lyapunov orbits -
As with planar Lyapunov orbits, the first thing we need to decide is 'when' we want to place a vessel in a vertical Lyapunov orbit - i.e., we need to choose a starting MJD for our orbit. Again for convenience, I'm going to choose the MJD of the example given in Dialing-up arbitrary Lagrange Point orbits, namely MJD = 52013.754909351. Again, there is nothing particularly special about this date - it just happens to be the date that I happened to extract from Orbiter 2016 when I was putting together the example for my earlier post. But any MJD will work just as well.
And, again, the next thing that we need to do is to establish some of the orbital parameters of the Moon's motion around the Earth on the afore mentioned MJD that we have just determined. In particular, we need to determine the semi-major axis; the orbital eccentricity and the true anomaly. As before, for the given MJD, we use the following Lua script to interrogate Orbiter's ephemeris system ad extract the relevant state vectors of the Earth and the Moon:
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 file 'statevectors.txt'. These state-vectors will be 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 into their Orbiter-friendly form. 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:
# calculate 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 finally 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 vertical Lyapunov orbits. These are defined such that 'da1' = 0.0. And with 'da1' = 0.0 we can set 'phi1' to zero for our convenience (although it doesn't really matter if we take any value between 0 and 2*pi).
This just leaves us with two free parameters, 'da2' and 'phi2', with which to define an vertical Lyapunov orbit. The latter is a phase angle which we are able to choose freely between 0 and 2*pi. In keeping with the approach adopted in Dialing-up a planar Lyapunov orbit, we'll choose phi2 = 0.0 if for no other reason than it is convenient to do so. The remaining parameter 'da2'>0 is a scale parameter which sets the size of the vertical Lyapunov orbit. As before, we have to be a little careful in our choice of value here. Because the underlying mathematical theory is based on a Lindstedt-Poincaré perturbation analysis, if we choose too large a value, the perturbation solution will not be a good approximation to the true solution of the equations of motion. So, how big is too big? Well, if we choose a value of da2 in the range 0.00 to about 0.08 we should achieve satisfactory results - but for values of 'da2' less than 0.04 or so, the underlying Lindstedt-Poincaré perturbation analysis is highly accurate. But how big is a vertical Lyapunov orbit with 'da2'=0.02? Well, we can roughly estimate the size of the Lyapunov orbit by calculating: da2 * 1.00 * 2 * 385000 km. This is the end-to-end length of the orbit. For da2 = 0.02, we estimate that the end-to-end-length of the orbit is around 15,400 km. With an upper bound of 'da2' around 0.10, the end-to-end length of the vertical Lyapunov orbit that we can reasonably approximate using our perturbative theory is around 75,000 km.
Here, we'll work with da2 = 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.00, 0.02, 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 ECI state-vectors:
Code:
x = 121321168.44146729
y = -5241818.015493814
z = -451601123.21481323
dx/dt = 1109.0182756245813
dy/dt = -65.54264423291947
dz/dt = 254.72878296539056The 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:
Then 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 prescribed vertical Lyapunov orbit.
IMPORTANT NOTE: Although the Lindstedt-Poincare perturbation analysis gives high-precision solution 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 some form of 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 you were wondering, here is what the prescribed Lyapunov orbit looks like:
