I've thought a bit more about getting to EM L1 from Low Earth Orbit (LEO) in Orbiter and this post sets out my take on what you need to do to get there.
I've tested this approach in Orbiter and, although it lacks pin-point precision, it is a perfectly serviceable way of getting to L1. Starting from a 200 x 200 km LEO, the L1 injection burn requires around 3,100 m/s. There are no mid-course corrections needed and, at apoapsis (and rendezvous with L1), a burn of around 700 m/s is needed to match the velocity of the spacecraft with L1. At this point, one is close to matching L1 speed and position, and after that, its just a matter of making small adjustments to velocity so that one eventually arrives 'on station' (to within a few metres) at L1. This manoeuvring doesn't take much Delta-V but it does take a while.
Things you need to be able to calculate
This is what you need:
1. You need to be able to calculate the coordinates of the Moon (relative to the Earth) for any time and date - cartesian, ecliptic coordinates. This calculation doesn't need to be high precision. but knowing the position of the Moon to within 100 km is about the right kind of accuracy needed here. There are a number of algorithms that you can download from the internet to do this - or you can look up the algorithms in books such as Meeus' "Astronomical Algorithms".
2. You need to be able to calculate the coordinates of L1 (relative to the Earth) for any time and date - again, cartesian ecliptic coordinates. Based on knowledge of the position of the Moon from '1' above; and the reference orbit calculations from
http://www.orbiter-forum.com/showthread.php?t=36110, this is doable.
3. You need to be able to calculate the coordinates of your spacecraft (relative to the Earth) for any date and time - again cartesian, ecliptic coordinates while in its parking orbit around the Earth. Here, I am going to assume that you have a spacecraft in circular 300 km x 300 km orbit, coplanar with the Moon.
4. You need to be able to calculate your speed and velocity with respect to L1 at any point in time. This information is used to identify when you need to execute the burn to match speed with L1; and to work out how much delta-V to apply. It can also be used to move progressively 'on station' once you have executed this velocity matching burn. I will attach the necessary lua script that calculates this below in a separate post.
This is the algorithm
Step 1: the transfer from Earth to L1 takes about 4.5 days. Call 't1' a date 4.5 days ahead of your current date. At 't1', calculate the position of the Moon and, hence, the position of L1 at time 't1'.
Step 2: calculate the distance of L1 from Earth on that date. Subtract 5500 km from it. (This is a crude adjustment to take into account 3-body affects close to apoapsis of the insertion trajectory.) Call this number 'ra'. Assuming that you are in a 300 x 300 km parking orbit, set 'rp' to 6671 km.
Step 3: Calculate the transfer time as:
[MATH]\Delta\tau = \pi\,\sqrt{\frac{(r_a+r_p)^3}{8\,\mu}}[/MATH]
where [MATH]\mu[/MATH] is the gravitational constant for the Moon.
Step 4: Subtract [MATH]\Delta\tau[/MATH] from 't1' and call this 't0'.
Step 5: Calculate the position of the spacecraft at 't0'.
Step 6. Calculate the angle between the position of the spacecraft at 't0' and the position of the L1 at 't1'. You want to execute the injection burn on the opposite side of Earth from where L1 will be). So, if this angle is not 180 degrees, make small adjustments to 't1' and repeat the above calculations until it is. 't0' then gives you your L1 injection date.
Step 7: Calculate the amount of dV required to execute the escape burn from:
[MATH]\Delta V = \sqrt{\frac{r_a}{r_p}\frac{2\,\mu}{r_a+r_p}}-\sqrt{\frac{\mu}{6671000}}[/MATH]
assuming, again, that you are in a 300 x 300 km parking orbit around the Earth.
Step 8: Set up a manoeuvre in TransX (or, if you prefer, IMFD) to execute the burn.
Step 9: execute the trans L1 injection burn and coast out to orbital apoapsis. This should get you to with a few hundred kms of L1.
Step 10: On approach to apoapsis, use information from the lua script to identify when you are close to apoapsis, and the use the velocity information to execute a velocity matching manoeuvre. To do this, you will need to orient the vessel in the lunar retrograde orientation and execute approx 650 - 700 m/s burn.
Step 11: By now you should be 'close' to L1 (to within a few hundred km) with a velocity that more or less matches that of L1. Now, use the distance/speed information from the lua script display to slowly manoeuvre eliminate the residual distance and speed difference. This may take a while, but you should be able to get distance differentials to with a few metres; and speed differentials down to around 0.1 cm/s without much difficult. Small manoeuvring thrusts are needed to make these adjustments.
Although this isn't an ideal way of getting to L1, it is a workable solution and avoids most of the complexities of three-body calculations relying instead primarily on 2-body Keplerian physics. In short, it works.
Because I'm sure that this description isn't entirely clear, I'll work on setting up a short video to show implementation of this approach.
---------- Post added at 06:05 AM ---------- Previous post was at 05:06 AM ----------
Here is the script for calculating the position and velocity of L1. To use, copy the script into a text file and call it, say. 'L1script.lua'. Save the script to Orbiter's 'Script' directory. Open Orbiter and then open a Lua Console window. Type in "run 'L1sript' ". You should see some blue text on the right-hand side updating in real-time. This text gives the position and velocity of the focus vessel relative to L1. The coordinate system is an IMFD style moon-entric coordinate system giving position and velocity vectors in terms of 'prograde', 'plane' and 'inward'.
Code:
note = oapi.create_annotation()
note:set_pos (0.80,0.09,1,1)
note:set_colour ({r=0,g=100/255,b=1})
note:set_size (0.7)
-- get the handle for the spacecraft
ves = vessel.get_focushandle()
-- get the handles for the Earth and the Moon
earth = oapi.get_objhandle("Earth")
moon = oapi.get_objhandle("Moon")
-- define a set of constants relevant to the Earth-Moon L1 libration point
GM1 = 398600440157821.0 -- gravitational constant for the Earth (SI units)
GM2 = 4902794935300.0 -- gravitational constant for the Moon (SI units)
GM = GM1 + GM2
MU1 = GM1 / GM
MU2 = GM2 / GM
alpha = 0.8369151948720568 -- the libration parameter for the Earth-Moon L2 point
goal = 0
while goal < 1 do
-- get the current location of the vessel
q = oapi.get_globalpos(ves)
p = oapi.get_globalvel(ves)
-- 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)
-- calculate the weighted average position and velocity of the Earth and Moon
com = vec.add( vec.mul( q_ear, MU1 ), vec.mul( q_mon, MU2 ) )
cov = vec.add( vec.mul( p_ear, MU1 ), vec.mul( p_mon, MU2 ) )
-- calculate the position of the Moon relative to the Earth
r = vec.sub( q_mon, q_ear )
v = vec.sub( p_mon, p_ear )
-- calculate some quantities that are used multiple times in the ensuing calculations
vsq = vec.dotp( v, v)
rln = vec.length( r )
rv = vec.dotp( r, v)
-- calculate:
-- 'e' - the eccentricity vector of the Moon relative to Earth
-- 'ecc' - the eccentricity
-- 'a' - the semi-major axis
-- 'nu' - the mean anomaly
e = vec.sub( vec.sub( vec.mul( r, vsq / GM ), vec.mul( v, rv / GM) ), vec.mul( r, 1.0 / rln ) )
ecc = vec.length(e)
a = GM / (2.0 * GM / rln - vsq)
nu = math.acos( vec.dotp(e, r) / ecc / rln)
if rv < 0 then
nu = 2.0 * math.pi - nu
end
-- calculate the unit vectors of a dextral reference frame aligned with the Moon's
-- orbital plane and orbital orientation:
-- 'xhat' - a unit vector pointing in the direction to the Moon's orbital periapsis
-- 'zhat' - a unit vector point normal to the Moon's orbital plane
-- 'yhat' - the third unit vector to complete the trio
xhat = vec.unit( e )
zhat = vec.unit( vec.crossp( r , v ) )
yhat = vec.unit( vec.crossp( zhat, xhat ) )
-- calculate some more intermediate values
k1 = a * (1.0 - ecc * ecc )
k2 = math.sqrt( GM / k1 )
cnu = math.cos(nu)
snu = math.sin(nu)
k3 = 1.0 + ecc * cnu
-- calculate the position of the Lagrange point in the dextral reference frame:
-- 'qx' - the position of the Lagrange point in the 'xhat' direction
-- 'qy' - the position of the Lagrange point in the 'yhat' direction
-- 'qz' = 0 by definition
-- 'px' - the speed of the Lagrange point in the 'xhat' direction
-- 'py' - the speed of the Lagrange point in the 'yhat' direction
-- 'pz' = 0 by definition
qx = alpha * k1 * cnu / k3
qy = alpha * k1 * snu / k3
px = -alpha * snu * k2
py = alpha * (ecc + cnu) * k2
-- caluclate the position of the Lagrange point in Orbiter's global reference frame
l2qa = vec.mul( xhat, qx )
l2qb = vec.mul( yhat, qy )
l2qc = vec.add( l2qa, l2qb )
l2q = vec.add( l2qc, com )
-- calculate the velocity of the Lagrange point in Orbiter's global reference frame
l2pa = vec.mul( xhat, px )
l2pb = vec.mul( yhat, py )
l2pc = vec.add( l2pa, l2pb )
l2p = vec.add( l2pc, cov )
-- calculate the distance from the reference L1 orbit
qrel = vec.sub( q, l2q)
prel = vec.sub( p, l2p)
-- calculate a new set of axes
rm = vec.sub( q, q_mon )
vm = vec.sub( p, p_mon )
xm_hat = vec.unit( vm )
ym_hat = vec.crossp( vm, rm )
ym_hat = vec.unit( ym_hat)
zm_hat = vec.crossp( xm_hat, ym_hat)
q1 = vec.dotp( qrel, xm_hat )
q2 = vec.dotp( qrel, ym_hat )
q3 = vec.dotp( qrel, zm_hat )
p1 = vec.dotp( prel, xm_hat )
p2 = vec.dotp( prel, ym_hat )
p3 = vec.dotp( prel, zm_hat )
-- print out the difference between the ship's psoition and the Lagrange point's position
msg = string.format("Position relative to L1\n-----------------------\nPrograde %.3f km\nPlane %.3f km\nInward %.3f km\n\nVelocity relative to L1\n-----------------------\nPrograde %.3f m/s\nPlane %.3f m/s\nInward %.3f m/s\n", q1/1000, q2/1000, q3/1000, p1, p2, p3 )
note:set_text(msg)
if goal > 0 then end
proc.skip()
end
---------- Post added 07-14-16 at 04:55 AM ---------- Previous post was 07-13-16 at 06:05 AM ----------
A couple of videos are now up on Youtube walking through implementation of the method described above.
LEO to Earth-Moon L1
