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

Thread Tools 
04122018, 06:56 AM  #1 
Orbinaut

Lambert Solver for Lagrange Points
For a while now, I've been thinking about the kinds of tools that one needs to navigate efficiently and effectivel in the vicinity of the Lagrange points of, say, the EarthMoon system. To be sure, Orbiter has a number of lesserknown (but quite sophisticated) tools in the Orbiter MFD stable  notably Lagrange MFD v1.5 for Orbiter 2016  which permit some basic navigation relative to the Lagrange points of the EarthMoon system and the SunEarth system. But none of these tools allow one to design optimal transfers in 'cislunar space'.
Surely we can do better than this. A comparison with twobody motion As most users already know, Orbiter has been very successful in enabling tools that allow for the design of efficient interplanetary transfers using the now very familiar 'patched conics' approximation . In particular, IMFD's Course Program can be used to design efficient Hohmannlike transfers between planets. It is able to do this because it makes use of a Lambert Solver to construct a Keplerian arc between a trajectory startpoint and a trajectory endpoint with a specified timeofflight. A Lambert Solver solves the twopoint boundary value problem "Given a starting point A, an end point B and a transfer time dt, what initial velocity is required at the point A to make this transfer". For the Kepler problem efficient solutions to this boundary value problem can be computed. And ny judicious adjustment of the departure date and arrival date, one can construct an optimal, twoimpulse transfer manoeuvre that minimises the amount of dV required to make the transfer. But what about motion in the vicinity of a Lagrange Point? In the vicinity of the Lagrange points of the EarthMoon system say, IMFD's Keplerian Lambert Solver built into its Course Program no longer works (because the gravitational field in the vicinity here is no longer Keplerian). And, yet, we do not have an equivalent Lambert Solver for cislunar space that we can use to carry out the same sort of trajectory optimisation that we would for the patchedconics model. For example, let's suppose that we are at some arbitrary position reasonably close to the L1 Lagrange point and travelling reasonably slowly relative to it. Let's also suppose that we wish to make a transfer to the Lagrange point that takes a specific time, dt. What velocity should our satellite have now in order to effect this transfer to L1? And what velocity will we have relative to the Lagrange point when we arrive? When is the optimal time to leave; and what is the optimal duration of our transfer? In the absence of a Lambert Solver for the restricted threebody problem, we struggle to answer basic questions of this kind. A different kind of Lamber Solver Evidently, the only way of answering questions of this kind is to have a Lambertlike Solver that has been designed to work in the vicinity of Lagrange points. Needless to say, I wouldn't be posting this note if a tool of this kind didn't exist. Let's make things a little more concrete by considering a simple (but arbitrary) example. And because we are thinking about navigation in the vicinity of Lagrange points, we are going to set this example up in the rotating reference frame of the threebody problem such that the location of the Lagrange points in that reference frame is fixed. In this rotating reference frame, imagine that we have a satellite located 4,000 km from L1 towards the Moon (the 'x' coordinate in our rotating reference frame); 4,000 km vertically above the orbital plane of the Moon in its orbit around the Earth (the 'z' coordinate in our rotating reference frame); and 4,000 km 'ahead' of the Lagrange point (the 'y' coordinate in our rotating reference frame). In other words, the satellite starts off at some arbitrary point about 7,000 km from L1. Let's also suppose, for sake of example only, that we wish our transfer to L1 to take a rather leasurely 71.3 days. Of course, there is no particular reason for having a transfer of this kind, but we choose these number somewhat arbitrarily simply to work through an example in this note. In the context of this scenario, we can ask a number of questions: * what is our transfer trajectory? * what velocity do we need to have now to send us on this trajectory? and * what velocity will we have when we finally arrive at L1? * if we vary our timeofflight, how do these things change? These questions parallel those that we would ask of a standard Lambert Solver. And we use answers to those questions to drive the search for optimality. So, if we can answer these questions for the threebody problem, we should be able to carry out the same sort of optimsisation tasks that we would normally apply to a standard Keplerian (twobody) Lambert Solver. And since I do have a solver of this kind I can, of course, use it to answer these questions. What is our transfer trajectory? First of all, what is the transfer trajectory that we need to effect the transfer from our (arbitrary) starting position to L1 such that the timeofflight is 71 days. Well, if we apply our new Lambert Solver to this problem, we find (without proof) that the threedimensional transfer trajectory looks like this: Clearly, the trajectory has a complicated structure. But, then again, this is a transfer orbit that takes 71.3 days (~ 2.5 lunar orbits of the Earth), so one might have expected that the transfer wouldn't be straightforward. So what's going on here in this trajectory. Roughly speaking, from the start point, the trajectory follows the stable manifold of Lissajous orbit around L1. Once 'on' the Lissajous orbit, the satellite circumnavigates the L1 Lagrange point six times before finally exiting the Lissajous orbit on the unstable manifold of that Lisaajous orbit that happens to pass through the Lagrange point arriving at L1 on the correct date. Of course, there is nothing particularly significant about the starting point or the duration of the transfer. They have merely been used here to demonstrate the capability of the Lambert Solver. As for the Lambert Solver itself, I'll focus on the mathematical details elsewhere but its foundations rest on sold maths and physics and, in effect, reduced the problem to an exercise in linear algebra. Although the end result is a trajectory solution that makes use of the stable, unstable and centre manifolds in the vicinity of the Lagrange point, the Lambert Solver does not solve the problem in terms of those structure. What velocity do we need to have now to send us on this trajectory? Naturally, we can find the velocity vector (in rotating coordinates) that we need to start ourselves on this trajectory. For our sample solution the xyz components of the velocity vector (in the rotating coordinate system) are as follows: * xdirection  34.65 m/s * ydirection  +3.03 m/s * zdirection  +3.19 m/s Just as with a standard Keplerian Lambert Solver, we are able to calculate the initial velocity vector with ease. However, whereas Keplerian trajectories are metastable so that once on a Keplerian arc, and absent perturbations, one stays on the same Keplerian arc, for the Lambert Solver for Lagrange points, we have a situation where we will need to undertake periodic corrections in order to remain on the prescribed transfer trajectory. However, this kind of 'station keeping' doesn't present any particular problem since we can periodically rerun our solver during the flight to calculate the minor correction needed to remain 'on track' for the entire 71.3 days. What velocity will we have when we finally arrive at L1? It is just as straightforward to calculate the velocity components upon arrival at the L1 Lagrange point. For our example, for our test problem, we find that the speed of the vessel upon arrival is: * xdirection  +8.75 m/s * ydirection  +19.02 m/s * zdirection  +18.31 m/s Upon arrival at L1, and assuming one wishes to 'park' at L1, this is the velocity vector that would need to be nulled out. If we take the length of the velocity vector, we find that the dV required to park at L1 is 27.8 m/s. If we vary our timeofflight, how does this change? We can also ask the question: is this solution optimal? Well, let's suppose that we only care about minimising the dV requirements at L1, given our starting point, what is the timeofflight that minimises this dV? At this point, we can use the same kind of techniques that we would use when using IMFD's Course Program to optimise the transfer trajectory. And without going into the details, but simply using the timeofflight as a control parameter, we find that if we target a timeofflight of 69.5 days rather than 71.3 days, we can reduce the dV requirements to 'park' at L1 from 27.8 m/s down to just 15.7 m/s. For those interested, the optimal trajectory is shown below: Of course, this optimality search was not intended to be exhaustive  but the point remains valid: once one has a threebody Lambert Solver, one has one of the basic instrument for finding optimal transfers in the vicinity of Lagrange points. Where to from here? By virtue of the trajectory planning tools that have been created thus far for Orbiter, Orbiter has largely been a vehicle for designing, optimising and executing patched conic trajectories. As such, it has been remarkably successful. But even though Orbiter includes a faithful rendition of Newtonian gravity for the Solar System, very few tools have been created that allow one to design, optimise and execute trajectories in the vicinity of Lagrange points. The Lambert Solver briefly showcased in those note goes some way to addressing that issue. As for the Lambert Solver, itself. At the moment it exists, and it works  but there is a substantial amount of work to be done to make this into a tool that is useful by other in the Orbiter community. Nonetheless, this is the longterm goal of this endeavour and this note is but the first small step towards achieving that end. 
Thanked by: 
04122018, 07:33 AM  #2 
Orbinaut

Not sure how your solver looks like in practice, but feel free to use the numerics of LEO Targeting for your needs (or to contribute to the code)  I'm currently at the point where I toyed some with the Earth/Moon Lagrange point dynamics and plan to expand on it. Basically it does fits over bruteforce numerical trajectory computations with arbitrary perturbations (3rd body gravity, J2 and J3,...)
Last edited by Thorsten; 04122018 at 07:37 AM. 
04122018, 07:57 AM  #3 
Orbinaut

Thorsten
I haven't described the solver in detail in part because the maths is probably a little too much for most on this forum; and partly because setting out the details will take considerable effort. But, in brief: The solver is configured as a singlestep very highorder symplectic discrete variational integrator using GaussLobatto quadrature as per the procedure described by Marsden & West (Marsden & West 2001) reconfigured as the solution of a boundary value problem. In effect, this fits a highorder (near minimax polynomial) to the motion of the satellite and, as such is capable of calculating the trajectory of a body in a arbitrarily complex gravitational field with machine precision over multiple orbits in a single time step. These integrators are fearsome and powerful beasts. The physical problem that is being modelled is the elliptical threebody problem in rotating, pulsating coordinates. In the vicinity of the Lagrange points, the equations of motion can be linearised and, in this case, and using the Marsden & West variational integrators, one can solve the twopoint boundary problem in rotating coordinates using linear algebra after having specified the transfer time. For the purpose of this note, I have restricted attention to the linear approximation. But this linear solution can be used as startingpoint for finding highly accurate solutions to the full, nonlinear elliptical threebody boundary value problem. At the moment, the modelling and calculations are done in Mathematica. Although these can be ported to c/c++ without much effort, I'm not sure how this endeavour fits with your LEO targeting tool  although I'm sure that one could solve a multiperiod perturbed Keplerian variant of the Lambert Solver for that case if one wanted to. 
04122018, 09:16 AM  #4 
Orbinaut

Okay, then it really seems you don't need any of my code. Good job accessing an interesting problem!

Thanked by: 
04132018, 12:54 AM  #5 
Scientist

It would be interesting to see your solver mated with my LagrangeMFD. The code in LagrangeMFD is a 4th order symplectic integrator, well able to plot really eccentric orbits that you find around LP's, and able to do hypothetical burn predictions with really good accuracy over several thousand seconds.
What it lacks is the ability to find burns into really beautiful orbits. An interesting solution would be to code your solver into an MFD, and then pipe results via ModuleMessaging into my utility to visualize the orbit and to execute the burns. It's nontrivial work though. 
04142018, 05:14 AM  #6 
Orbinaut

ADSWNJ
I'm broadly familiar with your Lagrange MFD  I've even used it a couple of times to construct ad hoc, lowenergy transfers from Each to the Moon and back again. Although I'm quite pleased with the Lambert Solver I've introduced in this thread, I'm conscious that I have still a fair bit of work to do on the mathematical side of things before I would contemplate housing the thing in any form of MFD. Having said that, I have in mind using this kind of tool to do just the thing that you have highlighted  i.e., "to find burns into really beautiful orbits." Two things are needed to do this. The first is to build a Lambert Solver of the kind described above; the second is to build a widget that describes mathematically the various kinds of Lagrange point orbits that one can have. On the plus side, the mathematical foundations of the latter component has already been largely done by Lei et al (Highorder solutions of invariant manifolds associated with libration point orbits in the elliptic restricted threebody system. Still, coding of Lei's work into a useful tool is still a fairly formidable (and time consuming) challenge. All that I can see is that I'm moving in the direction that you have outlined  but it is going to take a while to get a point where it makes sense to have meaningful discussions about how best to house in the Orbiter MFD stables.  Post added 041418 at 05:14 AM  Previous post was 041318 at 10:24 AM  To further illustrate what one can do with a Lambert Solver as outlined above (and to expand on comments in reply to ADSWNJ) consider the following scenario in the EarthMoon system: Imagine that we have a vessel initial at the point labelled 'Start Point' in the above diagram and it wishes (for some reason) to transfer to the planar Lagrange point orbit highlighted in red and finally achieving a final rendezvous with the target Lagrange point orbit at the point labelled 'A'. First of all, we need a way of describing the red target orbit; and second we need a way of describing the transfer from the Start Point to the rendezvous point 'A'. For simplicity for this illustration, we can work with the linear approximation to motion in the vicinity of the Lagrange point, in which case the motion in the xdirection (i.e., in the direction pointing from the Earth to the Moon) is given by: km and in the ydirection: km where and where is measured in days. From this, description we can determine any point on the target Lagrange point orbit and, moreover, the velocity that we need to have in order to stay on the target orbit. Lambert Solver solution for transfer time = 28.075 days Given this description of the start and endpoints, we can now feed this description into our Lambert Solver. To get this solver to work, we need to choose a transfer time. Let's try, for sake of argument, a transfer time from the Start Point to the rendezvous point 'A' of 28.07 days. The Lambert Solver solution (the think black line) is shown below: This shows that the transfer trajectory dropping down to a lower amplitude Lagrange point orbit before exiting to meet the rendezvous point with an approach speed of 13.18 m/s. Lambert Solver solution for transfer time = 28.966 days Let's try a different transfer time of 28.966 days, say, and rerun our Lambert Solver to calculate a new transfer trajectory with the new transfer time. The Lambert Solver solution (again, the think black line) is shown below: This shows that the transfer trajectory dropping down to a higher amplitude Lagrange point orbit before exiting to meet the rendezvous point with an approach speed of 14.78 m/s. Lambert Solver solution for transfer time = 28.623 days Now, we can play around with our transfer time to try and minimise our approach speed at the rendezvous point 'A'. It doesn't take long using a simple process of trial and error to estimate that if our transfer time is 28.623 days, then we can reduce our approach speed to a trivial 0.06 m/s. And with further tinkering we could probably reduce it even further. In this particular case, the Lambert Solver transfer trajectory solution is: What is happening here, is that because we have sought the minimum approach speed to the target Lagrange point orbit, we have indirectly found the stable manifold trajectory that connects our Start Point to the target Lagrange point orbit. Once we are on it, we will glide smoothy towards the reference orbit. And, once on it, we should need only minor stationkeeping corrections to keep us n the target orbit. Extensions Although the above has focused on a very simple Lagrange point orbit, we can apply exactly the same kind of technique to approach other, more complicated, Lagrange point orbits e.g., halo orbits. And although the above analysis has focused on the simple linear approximation, incorporating the various nonlinear corrections to calculating target and transfer orbits is a just a matter or incorporating the appropriate detail and complexity of the full elliptical threebody model. 
Thanked by: 
04142018, 01:03 PM  #7 
Scientist

It's such a shame that Dr Keith Gelling is no longer active in this forum. He would love to have discussed this topic.
In terms of the mechanics of coding such a solution, you basically have two choices: code the Lambert Solver as a freestanding C++ executable, or code it as an MFD. In both cases, you need the solver to run in multithreaded code, to unhook the deep calculation engine from the user interface. This technique is used in LagrangeMFD, so you can dial in a target calc time for its integrator  e.g. 2 seconds  without blocking the responsiveness of the Orbiter simulator  e.g. refreshing each 0.1 sec. Interested to see how you get on with this project! 
Thanked by: 
04162018, 10:44 AM  #8 
Orbinaut

Quote:
After the thread: https://www.orbiterforum.com/showthread.php?t=38788 I found the NAIF files for the Lagrange points of the SunEarth system: https://naif.jpl.nasa.gov/pub/naif/g...agrange_point/ If I do all the calculations told me by perseus in that thread, I get very different results from NAIF. Is there any “standard” way to obtain an accurate position of L1, L2 and L3 for any system? 
04162018, 11:57 AM  #9 
Orbinaut

Quote:
for where: and . and where is the mass of the first gravitating body; and where is the mass of the first gravitating body. This expression is quite general and works for all restricted threebody systems. The rest of the note referred to by dgatsoulis uses this value to calculate the position of the L1 point in Orbiter's inertial reference frame. First, you have to calculate: The first of these is the eccentricity vector (and the length of which is the orbital eccentricity, ); the second is the setmajor axis; and the third is the true anomaly of the elliptical orbit of the second gravitating body around the first gravitating body. In turn these quantities depend on the state vectors of the primary and secondary gravitating bodies  but you can always get this information by using an appropriate ephemeris calculator. If you work through calculations in dgatsoulis' referenced thread, it seems to me that you should be able to calculate the position of L1 (or L2 and L3) in the context of the elliptical restricted threebody problem with some precision. 
Thanked by: 
04162018, 02:51 PM  #10 
Orbinaut

I tried to use the Keithth’s method, but it seemed overly complex to me and not good for the heavily perturbed systems (like Pluto and Jupiter).
But now I coded the Keithth’s Lua snippet for the SPICE library; the result is good for the SunEarth system (the ERTBP L1 point is no more than 11000 km away from the L1 point obtained from the SPICE file). But what about ERTBP method in other systems? I thought to use the gravitationalcentripetal accelerations method explained in the other thread because it seems an “always good” method, but it also seems that I don’t understand what happens in a Lagrange point, because the error in my calculations is too big. This is a very simple method because I just need to find a point where centrifugal acceleration = gravitational acceleration, is that correct? 
04172018, 01:22 AM  #11 
Orbinaut

Cristiapi
Strictly speaking, Lagrange points are only welldefined 'points' in the CR3BP and the ER3BP  i.e., only in specific classes of restricted threebody mathematical models. In a general nbody system, the Lagrange points do not (strictly speaking) exist since you will never get a stable point where gravitational and centripetal forces exactly cancel. However, if perturbations are weak enough we can still make use of notional Lagrange points. Take, for example, the EarthMoon system which is quite heavily perturbed by solar tidal forces. Tidal forces fall of as 1/r^3 rather than the usual 1/r^2  so, at the EarthMoon system, solar tidal forces are quite weak but still strong enough to distort noticeably the Moon's orbit around the Earth. Here, it seems to me that mission designers use threebody concepts from the ER3BP to construct 'reference orbits'. In principle, if the EarthMoon system were a true restricted threebody system then tiny amounts of dV would be required to remain 'on station' on the reference orbit. In practice, mission designers have to allow for the fact that the perturbations mean that the satellite has to do more work to stay on the reference orbit to counter the effect of the perturbations and so more dV is used. As the magnitude of the perturbations increases, the satellite has to do more and more additional work to stay on the reference orbit. And, at some point if perturbations are high enough, the concept of a Lagrange point and a reference point breaks down since the dV cost of staying on it is simply too high. At this point mission designers need to go back to the drawing board and work, I presume, with more general nbody calculations. wrt the position of L1 for the SunEarth system, are you using the state vectors of the Earth for your Lagrange point calculations or the state vectors of the centre of mass of the Earth and the Moon? Because the SunEarth L1 point is 1.5 million km from Earth, it is better to model the combined gravitational contribution of the Earth and the Moon as if it were emanating from a point source located at the centre of mass of the EarthMoon system. Including the mass of the Moon in the calculations increases the mass of the secondary by about 1% and so will change the location of the L1 point by around 1%  i.e., by around 15,000 km. My guess is that you haven't included the Moon in your calculations and that doing so may well eliminate the bulk of the gap between your calculations and the spice file locations. You mentioned 'the heavily perturbed systems of Jupiter and Pluto'. In what respect do you regard these as heavily perturbed systems? 
Thanked by: 
04172018, 09:23 AM  #12 
Orbinaut

Quote:
Quote:
I did many simulations for the stability of L1 and L2, but it was 1 year ago... Anyway, the result was as you say: Quote:

04172018, 09:55 AM  #13 
Orbinaut

For PlutoCharon system, as you say, L2 is indeed heavily perturbed. Basically, trying to calculate the location of L2 is fairly meaningless here. L1 for PlutoCharon ought to behave a little better, though. But still, because the mass of the other moons aren't negligible, a satellite parked in the vicinity of the notional ER3BP L1 point is still going to get kicked around a lot.
wrt the Lagrange points of the moons of Jupiter  Io, Europa, Ganymede, and Callisto. Yes, perturbations (due in large measure to Jupiter's J2 field) are likely to force higher stationkeeping costs than one might ideally like. Callisto, in particular, ought to be reasonably well behaved. It's the furthest out so the J2 perturbation should be minimal. And if one replaces the state vectors of Jupiter with the state vectors of the centre of mass of Jupiter + Io + Europa +Ganymede, that also ought to improve things a little more. wrt "if I use EMB instead of just Earth I get exactly the same JPL result (the difference is on the order of 104 m)", that's probably good enough Last edited by MontBlanc2012; 04172018 at 10:14 AM. 
Thanked by: 

Thread Tools  


Quick Links  Need Help? 