OK, time now to say a few words about how to convert the PyKEP / PyGMO solution to the Aldrin Cycler problem into a 'high fidelity' trajectory for Orbiter 2010.
Let's start with the assertion that the PyKEP / PyGMO serves as a template or cartoon sketch of the trajectory: it sets out the basic structural elements of the trajectory but the underlying assumptions in PyKEP are of a zero-radius Sphere of Influence, 2-body Keplerian motion - i.e., a 'linked conics' approximation. If one were to try and fly this trajectory in Orbiter, one would find that one would need to apply plenty of mid-course corrections to keep the spacecraft 'on track'. Converting to a 'high fidelity' trajectory basically means tweaking the trajectory slightly so as to (with a high degree of precision) remove almost entirely the need for those mid-course corrections. What this means is that the high-fidelity trajectory will accurately model the 'true' trajectory of the spacecraft so that staying on the trajectory does not introduce any non-Newtonian forces on the spacecraft.
In other words, suppose that one had a Cycler craft being forced to run along the rails of the high-fidelity trajectory track. Like the planets in Orbiter, the Cycler is not in 'free fall'. The ephemeris solution forces it to run along a particular trajectory regardless of the true gravitational forces acting upon it. Imagine now an undocked craft immediately alongside the Cycler craft allowed to free fall in the gravitational field. Then, if the high-fidelity trajectory has been calculated correctly, there should be no change in the position of the free-falling spacecraft relative to the Cycler as it runs along its Cycler trajectory. Of course, where Deep Space Manoeuvres are required, this notion breaks down and there is a 'jump' in the speed and direction of the Cycler which is most decidedly non-ballistic. But for vast majority of the trajectory, the Cycler trajectory should accurately simulate a ballistic trajectory.
Realistic ballistic trajectory are important
This notion of accurately simulating a ballistic trajectory is most important near planetary rendezvous. Imagine that you are in a craft attempting to dock with the Cycler and suppose that you have worked your way to with, say, 100 m of the Cycler and have 'nulled out' your relative velocity with respect to the Cycler. A bit like docking with the ISS, what you hope and expect to see happen is the Cycler remain more or less stationary with respect to you. Not only does this make docking easier but it is also physically realistic. What you do not want to see is the Cycler accelerate away from you for no good reason simply because the Cycler is being forced to run along a set of rails that do not truly reflect the gravitational forces acting upon it. For you, in the nearby spacecraft, it would look as if some mysterious force were acting upon the Cycler and to dock with it you would have to continually match this mysterious acceleration. It would be weird, frustrating and completely unrealistic. Consequently, for the Cycler concept to work, we need to calculate an accurate, realistic ballistic trajectory for it.
Spicing and dicing
So, how do we go about doing that? The first step is to slice and dice the trajectory up into a series of smaller pieces. We note that the Cycler trajectory is punctuated by a series of DSMs. These occur in deep space - a long, long way from any planetary bodies. In between these DSMs, there is a ballistic arc that connects one DSM to another and with one or more ballistic planetary encounters. Between those planetary encounters where there is no DSM, we can imagine that there is a DSM at the half-way point along this arc - albeit one of zero magnitude. In this way, we can break up the entire Cycler trajectory into a series of ballistic arcs with exactly one planetary encounter in between. Rather than examine the whole trajectory in one go, we can develop ballistic trajectories for each of these arcs in turn and then splice them together. The net result will be a trajectory that is almost entirely ballistic everywhere except at the splice points where an instantaneous change in velocity (in effect, a DSM) takes place. In the PyKEP /PyGMO Cycler trajectory solution identified above, we have eleven such ballistic arcs - six contain encounters with Earth; and the remaining five encounters with Mars.
OK, so let's look at one of these encounters - one that goes from DSM - Earth - DSM, say. The PyKEP /PyGMO trajectory solution defines the 'when' and 'where' of the DSMs, the end-points of the ballistic arc. Now, PyKEP uses a slightly different coordinate system from that used in Orbiter but because these manoeuvres take place in deep-space a difference of 1,000 km, 10,000 km or even 100,000 km isn't going to make a lot of difference. So, the first thing that we can do is calculate the time and position of the DSMs from the PyKEP solution and immediately transcribe them into Orbiter's global reference frame. Now, we have what is known as boundary-value problem - we know where/when we start and where/when we finish; and we want to know the path we need to take in between. In many ways, this is similar to the classic 2-body Lambert problem. Here, though, we have a more complicated arrangement in which there is a generally low-altitude ballistic flyby of Earth in between. Moreover, we need to take be careful to take into account the perturbations induced by other Solar System gravitating bodies and so we need to use a full n-body gravity model. Here, simple semi-analytical solutions such as a Lambert Solver aren't really of much use. Instead, we need to use a high-fidelity n-body integrator to help us calculate trajectories. Fortunately, I have one of those.
The right ephemeris colution
What does the n-body integrator need to work effectively? To calculate a ballistic trajectory that would be calculated by Orbiter's own integration engine, we need to use the same gravity model as Orbiter. Since Orbiter uses a Newtonian gravity model, we need to use the same coordinate system as Orbiter, and we need to know to use the same planetary positions that Orbiter uses. In other words, we have to use the same planetary ephemeris model that Orbiter uses. Now, Orbiter planetary ephemeris is based on the VSOP87 ephemeris solution. This ephemeris can be downloaded from this code generator:
http://www.neoprogrammics.com/vsop87/source_code_generator_tool/
and so, for the most part, building an ephemeris model for Orbiter amounts to a 'cut and paste' job from this code generator. However, for the Earth, life is a little more complicated. The VSOP87 series calculates the position of the Earth-Moon barycentre and not the position of the Earth,
per se. So, what we need to do here is to also calculate the position of the Moon relative to the Earth in the same way that Orbiter does - namely, by using the ELP82 ephemeris solution for the Moon. Once we have that, and the VSOP87 solution for the position of the EarthMoon barycentre, then we can back-solve to find the separate positions of the Earth and the Moon in Orbiter's global coordinate system.
With these ephemeris solutions and a symplectic integrator, we are good to go.
Completing the ballistic arc
So, how do we actually go about building the ballistic trajectory between the DSMs? Well, we start in the middle - at the planetary encounter. Let' suppose that we can guess our periapsis radius and the orientation of our encounter plane. Let's suppose also that we can guess our velocity at periapsis and the time when the encounter occurs. Using our n-body integrator, we can then propagate forward to the time of the second DSM at the end of the ballistic trajectory and we can see where we end up. In the same way, we can use our integrator to propagate backwards in time to the time of the first DSM at the start of the ballistic trajectory and see where we end up. Now, unless we are very good at guessing, there is going to be mismatch between where we the n-body integrator calculates that the DSMs take place and where the PyKEP / PyGMO trajectory design says that they should. So, we update our guess and tray again. And we keep updating our guess until we guess right. The formal procedure for doing this a multidimensional Newton-Raphson root-finding calculation but basically it consists of making a series of intelligent guesses and then updating those guesses.
Once one has solved this problem (and guessed right), one then has a fully ballistic arc that goes from one DSM to another DSM with a planetary encounter in between. If one repeats this exercise for each of the separate 11 legs of the Cycler trajectory, then one has defined the high-fidelity, almost entirely ballistic Cycler trajectory. And once this is done, the sole task to remain is to encode this trajectory as a C function call that will return the Cycler position and velocity for any time between the start of the trajectory and the end of the trajectory.
Do we have enough free variables
In solving problems of this kind, one of the things that we need to work out is do we have enough free variables to solve this problem. In this case, we have six constraints that we need to satisfy:
1. the n-body integrator has to integrate forwards and backwards from the planetary encounter so that the position of the two ends of the trajectory match the PyKEP / PyGMO DSM positions - three constraints apiece.
What about the variables? The forward/backwards integration of the ballistic trajectory is fully defined by planetary periapsis encounter. What variables define this encounter? We have:
1. three variables that define the planet-centric location of the point of periapsis,
2. two variables that define the velocity of the Cycler in the plane that is perpendicular to the position vector that connects the planet centre to the point of periapsis, and
3. the MJD of periapsis.
In total, then, we have six free variables to solve six equations. This should give us enough free variables with which to solve the boundary-value problem. To start the guessing game, we need an initial guess. Again, we can use the PyKEP /PyGMO solution to provide this
It's going to take a while
All of this amounts to some fairly heavy-duty computations of my laptop. And it will take a while. First step, though is to understand the ELP82 ephemeris solution and make sure that I have an Earth-Moon ephemeris tool that closely matches that of Orbiter's. I have done this in the past for the Galilean moons by implementing the Galsat routines used by Orbiter. Her, as there, I expect that the result will be an online Earth-Moon ephemeris tool capable of replicating Orbiter's ephemeris calculations to with a few metres.
