Building an ephemeris for Hyperion

[Keithth G]

In an earlier post, I asked 'martins' a question:

"On another issue: in an earlier comment relating to making public integrators you mentioned that "this would come in handy in various places, e.g. for mission planning tools, or for creating our own power series ephemeris solutions for celestial bodies (I'm still looking for a solution for Hyperion)." What by way of ephemeris solutions are you looking for?"


To this, 'martins' replied:

"Essentially I was thinking about a machinery that can generate a series expansion of a perturbed orbit from a numerical solution over a reasonably long period centered at the present - essentially what VSOP does for the major planets. In the case of Hyperion, this would be referenced to the true position of Saturn, or barycenter position of Saturn and its moons (and ideally in the same frame of reference as VSOP, to avoid the frame transformations you mentioned earlier).

As far as I know, Hyperion's may be a particularly difficult orbit to model, because of resonances with Titan and high eccentricity, so a series expansion diverges either very quickly or requires a large number of terms. But since you have already created a precise numerical gravity simulator, it might be interesting to use it to see just how well it can be approximated with a series solution."


This exchange got me thinking about how one would go about building a VSOP-like ephemeris of, say, an 'irregular' body such as Hyperion. After some reflection, and a little bit of a literature search, I decided that the mathematical machinery to do this exists and so I thought I would 'give it a go'. So, this post is a 'place-holder' for a number of (yet to be produced) results arising from this exercise. Results will be posted here as I work through this exercise.

For those that might not be familiar with the vocabulary, an ephemeris is essentially a mathematical 'look-up' table that encodes the position of the Sun, the planets and moons. It serves two purposes, the first of which is entirely practical: although a lot of effort may go into putting the ephemeris together in the first place, ephemeris users need only employ a simple algorithm to extract the encoded information. An ephemeris, VSOP87, lies at the heart of Orbiter. It is used by the numerical integrator to work out the location of the principle gravitation bodies in the Solar System. The second is more theoretical. Because of the theoretical underpinnings of the ephemeris construction, the ephemeris solution may yield useful information about the nature of perturbations - e.g., the rate of perihelion advance of an orbit - that would otherwise not be self-evident. For Orbiter, this is less relevant - but in many respects the more interesting aspect.

There are two kinds of ephemerides in common usage. The first is the VSOP series, principally VSOP87 and ELP series, produced by the Bureau des Longitudes. These encode the positions of planets (and the Moon) as a series of cosine terms which, when summed, yield the planetary positions with high accuracy. The full VSOP87 solution contains many thousands of such terms for the major planets and although a little daunting when first encountered are actually quite simple (and fast) to use in practice.

The other kind of ephemeris is the 'DE' series produced by Nasa's Jet Propulsion laboratory - e.g., DE405/406. These ephemeris series are more recent and are used for mission planning purposes. However, the encoded solution of planetary positions consists of an unwieldy set of Chebyshev interpolation schemes. Data files for the full ephemeris solution are extensive to such an extent that they make the extensive VSOP cosine tables look positively succinct. Although, it is older, the VSOP87 is still widely used since for, most purposes, it is sufficiently accurate - and far more compact.


In this thread, then, I'm going to focus on the construction of a VSOP-like ephemeris solution - that is, a cosine table for the saturn-centric position of Hyperion over, say, one or two centuries from the present.

And why is this such a hard task? Surely, fitting a few cosine terms to an integrated solution of Hyperion can't be that hard. Basically, the issue is largely one of accuracy: to produce an ephemeris capable of accurately determining the position of Hyperion one hundred years from now requires a great deal of precision in determining the values of the ephemeris cosine terms: being off by 1 part in 10^10 in the frequency of a cosine term can lead to inaccuracies of thousands of kilometres. And then there is the sheer number of terms that may need to be included. Because the forces perturbing Hyperion's motion are considerable, and because it is a high eccentricity orbit, and because it has a resonance with Titan, the number of cosine is terms needed to accurately determine Hyperion's position may be very large. Or, indeed, Hyperion's orbit may be chaotic and predicting its position one hundred years from now is impossible. At any roads, having thought about it, the construction of a realistic and accurate ephemeris for Hyperion is not a trivial thing to do. (And if you don't believe me, give it a go yourself.)

In fact, the mathematical and computational challenges that this task imposes are such that it doesn't make much sense to attack the problem with a direct frontal assault. Better, methinks, to learn how to walk first by practicing on simpler, less demanding problems.

So, this is my general plan of attack: to work up to the task of an accurate ephemeris for Hyperion, by working progressively through a series of progressively more challenging problems. The specific sequence that I propose to work through is as follows:

1. Determine a VSOP-like ephemeris for a high-eccentricity two-body elliptical orbit around a single stationary gravitational source, e.g., the Sun. Although there is no need for an ephemeris solution here, it contains many of the same challenges of the Hyperion problem without making the problem overly complicated.

2. Determine a VSOP-like ephemeris for the Sun/Earth/Moon system. This is now a three-body problem and contains a number of significant perturbations that lead to phenomena such as 'orbital precession'. Although quite simplistic by modern standards, until well after the Lunar landings, ephemeris based on this three body problem were central to lunar mission planning.

3. Determine a VSOP-like ephemeris for the Sun/Saturn/Titan/Hyperion system. This is now a four body-problem. It captures the essence of resonance between Titan and Hyperion as well as Hyperion's high eccentricity orbit. It ignores significant perturbations due to the gravitational influence of the other planets - principally, Jupiter.

4. Determine a VSOP-like ephemeris for Sun/Jupiter/Saturn/Titan+other moons/Hyperion. This will be a semi-realistic ephemeris for Hyperion. It will include perturbations from Jupiter as well as those from the rest of Saturn's moons. It will not, however, be consistent with the VSOP87 ephemeris solution used by Orbiter

5. Determine a VSOP87-like ephemeris for Hyperion based on the full VSOP87 solution. The end result.


I've already completed Step 1 - and I'll report on this next. Step 2 is underway.

 

[Keithth G]

A toy model - ephemeris of Keplerian "2-body" motion

In this post, I'm going to go through the formal exercise of building a VSOP-like ephemeris for Keplerian 2-body motion. This is a dry-run for building an ephemeris for more complicated situations, but it serves to highlight the general procedures that one has to go through in order to build the ephemeris. The ultimate goal is to build an ephemeris for Hyperion, one of Saturn's moons, but to get there I have to develop a suite of techniques and tools that will allow me to complete that task. Examining 2-body motion is the first (baby) step towards that goal.

The specific Keplerian motion for which I want to build an ephemeris is for an object moving on an elliptical trajectory about the Sun. The perihelion of the orbit is close to the orbit of Earth; and the aphelion is close to the orbit of Mars. To be more specific, at time some arbitrary time, t = 0, we are given the information that the heliocentric Cartesian coordinates (referenced to the ecliptic plane) of the object are as follows:

x: 1.0000 AU
y: 0.0000 AU
z: 0.0000 AU

v_x 0.0000 AU/day
v_y -0.0190 AU/day
v_z -0.0010 AU/DAY

These numbers are arbitrary. The first three numbers give the position of the object at time t = 0; and the last three numbers give the velocity of the object at time t = 0. The z direction is assumed to be perpendicular to the plane of the ecliptic. Calculations will be done in Gaussian units - i.e., distance will be measured in Astronomical Units (AU), time will be measured in days, and so velocities will be measured in AU/day.

In this particular, example I'm going to pretend that the planets don't exist so that our 'Solar System' for this problem is very simple: it consists of the stationary Sun at the centre of our x-y-z coordinate, and our object moving in an elliptical orbit with these known co-ordinates at time t = 0. There are no perturbing forces. The task, then, is to build an ephemeris of the orbit of our object for one full century (36500 days) from time t = 0. And we want the accuracy of the ephemeris to be accurate to within about 1 km throughout all of that century. Moreover, we are going to pretend that we know nothing about elliptical 2-body motion and build our ephemeris in Cartesian coordinates ( x-y-z coordinates, rather than in terms of the more natural osculating orbital elements). And the ephemeris is going to consist of a (hopefully) small number of cosine terms consisting of amplitude, frequency and phase which, when summed, will yield the Cartesian coordinates of the object in its orbit with the desired degree of accuracy. This is essentially what a VSOP-like ephemeris does, and this is what we will do here.

Before building the ephemeris, it may be useful to get a feel for what this orbit looks like. Essentially, it is an elliptical orbit. The semi-major axis of the orbit is around 1.29 AU, it is inclined with respect to the plane of the ecliptic by about 3 degrees, and its eccentricity is around 0.223. The orbital period is around 533 days. In one century, our object orbits the Sun about 70 times. The orientation of the orbit with respect to the ecliptic doesn't matter - but what is important is that the eccentricity of the orbit is reasonably high - certainly nowhere near parabolic, but large enough to provide a reasonable test of the ephemeris engine being tested in this toy problem.

So, this is the ephemeris that we want to build. How does one go about building it?

Basically, there are three steps in the process.

1. Work out the equations of motions of the object for which we want to build an ephemeris. In this case, the equations of motion are simply the equations of motion for a body in an r^-1 potential.

2. Use a numerical integrator to integrate the equations of motion of the object over the period for which one wants to build an ephemeris. By changing times steps in that integration ensure that the numerical integration is carried out over the period of the ephemeris with the required degree of accuracy. In this case, the test of accuracy is that by starting from the same initial conditions (given above), does the integrator yield the same final position of the object one full century later to within, say, 10 meters if the time-step is doubled. In my case, I already have a fourth order symplectic integrator (left over from an accuracy test of Orbiter's integration engine). Although quite a simple integrator, it is a lattice integrator based on integer arithmetic and it can be pushed hard so as to eliminate the rounding errors that would normally accumulate if it were to use floating point arithmetic. [Aside: At some point, I will post details of this integrator - together with details of a similar 6th and 8th order versions - but not today]. With time steps of around 1 hour, so 1 million time steps over the course of one century, this simple fourth integrator yields a 'converged' trajectory accurate to < 1 m. Rather than present the next step of the ephemeris engine a time series with one million data points, this was reduced to just 3650 data points - i.e., one data point every ten days.

3. The third and crucial step in the construction of the ephemeris is to take the converged time series for the x, y and z coordinates and, in effect, perform an Fourier decomposition of those time series. For those with an engineering or scientific background, this doesn't sound particularly hard. And the principle of what we are trying to do here isn't difficult to fathom. The real issue is one of numerical accuracy: we need to determine the coefficients of the Fourier decomposition to very high accuracy if we are going to have an ephemeris accurate to within 1 km throughout the whole of the century. Just to get a feel for the magnitude of the problem. Very roughly, our orbiting object will travel about 6,000 million kilometres - and we want our ephemeris to be accurate to just 1 km.

So, how do we do this Fourier decomposition?

The first thing we note is that standard techniques such as the Fast Fourier Transform (FFT) are not up to the task. Basically, the accuracy with which frequencies can be determined falls off as 1/T, where T is the length of the time series. Since we want to determine frequencies to about 1 part in 10^10 (or thereabouts), we would have to use an extremely long time series to get the desired results. So, we can rule these 'basic' techniques out.

Then there is the 'frequency map analysis' of Jacques Laskar. This was a technique that he pioneered in the 1990s to demonstrate the chaotic nature of the long-run evolution of the Solar System. (See, for example http://arxiv.org/pdf/math/0305364.pdf). With this frequency map analysis, he claims, he can determine frequencies with an accuracy that falls off as 1/T^4 - i.e., very much faster than the FFT. Basically, the procedure that he uses is to use a modified Fourier transform of the time series to extract information about the dominant frequency. Then remove the frequency from the data and orthoganalises the remainder. And then repeat the process until one has extracted all of the relevant frequency information. Undoubtedly, this technique works, but I am not very familiar with it.

The third technique that I know of, and the one that I have used here, is a Bayesian spectral analyser. This uses Bayesian probability theory to estimate frequencies in batches. Whereas both the FFT and 'frequency map analysis' effectively single frequency analysers, the Bayesian approach allows one to optimise over frequency estimation for a batch of frequencies - ten, one hundred or even one thousand - at a time. And it is this 'batch' approach to frequency estimation that allows one to get very accurate estimates of frequencies. Largely because I'm more practised with Bayesain methods and would have to invest a lot of time to learn Laskar's frequency map analysis, this is the approach I've used here. Now, I could write a long note on how this method works - and if anyone is sufficiently interested, I will provide - but in this note I'm just going to focus on the basic physics and the results.

As for the physics, let's begin by doing the "poor man's" spectral analysis and simply take the FFT of the x-coordinate of the time series of our object's orbit (the 3650 values spanning one full century) and look at the spectrum:



What we immediately notice is that our orbit 'rings like a bell'. There is a fundamental frequency at around 0.0118 (corresponding to the orbital period of 533 days - i.e., 533 = 2 * pi / 0.0118); and there are series of harmonics at integer multiples of the fundamental frequency. These harmonics are visible to fourth order but they can easily be tracked down in the dynamics to 15th order (and beyond) if one wishes to push the analysis hard enough. This kind of 'fundamental frequency' and 'harmonics' is a general feature of physical systems of this kind. This knowledge can be used to good effect to refine one's estimate of the fundamental frequency (a.k.a. inverse orbital period ) and the harmonic frequencies. And this is exactly what my Bayesian spectral estimator does. It estimates all of these frequencies in one-fell swoop (at least down to the 12th order harmonic which is all that was necessary to achieve the 1 km accuracy target).

So what were the results? Well, first of all the Bayesian spectral analyser produced a very precise estimate of the fundamental frequency. Based on just the x-coordinate time series (with just 3650 data points), the entirely empirical estimate of the fundamental frequency was:

[MATH] \omega'_0 = 0.0117741988579948 [/MATH]
Based on well-established Keplerian physics, the inverse of the theoretical orbital period (multiplied by 2*pi) can be calculated to be:

[MATH] \omega_0 = 0.0117741988579962 [/MATH]
which is remarkably close to the empirical result. Personally, I regard this as a little computational triumph and was quite pleased when this result popped out of the analysis.

Now, with knowledge of the fundamental frequency, it becomes straightforward to obtain maximum likelihood estimates of the amplitudes of the sine and cosine terms of the fundamental frequency and harmonics. And when converted to cosine amplitude and phase, our ephemeris construction for the x-coordinate is complete. The first five terms of this ephemeris for the x-coordinate are given below.

[MATH] \begin{array}{crrr} \text{Order} &\text{Amplitude (AU)} & \text{Freq. (rads/day)}& \text{Phase (rads)} \\ 0 & -0.4313366223 & 0.0000000000 & 0.0000000000 \\ 1 & 1.2635578190 & 0.0117741989 & -0.0000000000 \\ 2 & 0.1390423702 & 0.0235483977 & -0.0000000000 \\ 3 & 0.0229746144 & 0.0353225966 & -0.0000000000 \\ 4 & 0.0045011225 & 0.0470967954 & -0.0000000014 \\ \end{array} [/MATH]
The first column denotes the order number of the harmonic. 'Zero' indicates a constant term; and '1' indicates the fundamental. '2', '3' and higher order denotes the order of the harmonic. The second column gives the amplitude (in AU) of the cosine term corresponding to that order; and the next two columns give the frequency and phase of the cosine term. The 'full' ephemeris for this orbit contains contributions down to 12 order but limitations on posting maths in this blog limits the number of lines that I can include here.

Having completed the analysis of the x-coordinate, one can repeat the analysis for the y and z components. Doing so completes the ephemeris construction.

So that interested readers can see what the full ephemeris solution looks like, I've attached an .xlsx file to this post (somewhere). Feel free to download and examine.

Although this was but a simple ephemeris, it works remarkably well for the task at hand. The next step is to look at a more complicated physical system where perturbations are most decidedly significant - the Sun/Earth/Moon system.

 

[Keithth G]

Building a lunar ephemeris (two dimensions)

This is the next update on an on-going effort to build an ephemeris for Hyperion, one of Saturn's moons and locked in a highly elliptical orbit with a 3:4 resonance with Titan.

So, far I've outlined a general methodology, built appropriate tools, and applied them to the simple case of Keplerian 2-body motion. All has gone well.

The next step, though, is more challenging: namely the construction of an ephemeris for the Sun/Earth/Moon system. Those with a feel for history will immediately recognise this as a tough problem. It was first studied (unsuccessfully) by Newton and has been studied by some of the greatest physicists and mathematicians in the intervening three and a half century. To a degree, it is still an area of active research. To say that the construction of an accurate ephemeris for, say, the Moon is what may be reasonably be called 'a hard problem'. And it is a hard problem because the Moon's motion through the night sky is very complicated filled with a host of surprising nuances.

Although the basic motions of the Moon have been known since antiquity, developing an accurate ephemeris for the Moon has taken science a few hundred years. Fortunately three centuries of hard slog culminated in an ephemeris of sufficient accuracy to take Man to the Moon in the late 1960s. It would be sheer hubris on my part,however, to claim that I'm going to construct an ephemeris of comparable accuracy in this note. But
it is possible to produce an ephemeris of reasonable accuracy in relatively short order. Here, my goal is to produce an ephemeris of the Moon's position (with respect to the Earth's geocentre) to an root-mean-square accuracy of around 5 km over the course of one Julian century. As it turns out, once can do this with under 50 cosine terms for each of the x, y and z coordinates. But it is also clear that to develop an ephemeris to resolve the Moon's motion to, say, 1 km will require a few hundred more terms. Now, it is quite possible to do this with the tools at hand but, frankly, its a tedious challenge with very little benefit - particularly given that this is just a warm-up exercise for calculating an ephemeris for Hyperion; and given that very accurate ephemerides for the Moon already exist. So, here, I've (thankfully) limited myself to an accuracy of just 5 km or so.

My analysis of the motion of the Moon will be undertaken in two steps. The first will be to take the Sun/Earth/Moon system and squash it so that it motions of the three bodies all lie in the same plane. This removes a degree of freedom from the problem and reduces the problem's complexity. The second step will be to repeat the exercise, but this time in three dimensions. In both cases, and to further simplify the problem, I've assumed that the Sun is sufficiently massive that its motion is not affected by the motion of the Earth/Moon system

The Moon's motion in 2 dimensions
Before going though the exercise of building an ephemeris for the Moon in the case were all motion is squashed down into a plane, it's worthwhile going though a little counting exercise to see if we can divine the kind of motions that we expect the Moon to exhibit. Let's begin by counting degrees of freedom.

OK, we know that we have three bodies in our mini Solar System: the Sun, the Earth and the Moon. By assumption, we're going to say that the position of the Sun is fixed, so there are no degrees of freedom associated with its motion. This leaves the Earth and the Moon. For the Earth, moving in the plane of the Ecliptic, we need two variables to describe its position (an x-coordinate and a y-coordinate, say) and we need two variables to describe its velocity - so four in total. The same is true for the Moon - so another four. This means that overall, we need eight variables. However, we know that for this system, the total energy (kinetic and potential) of the combined system is constant, so we lose one degree of freedom. We also lose another degree of freedom because angular momentum is also conserved. So, overall, we expect to be able to describe the motion with just six variables.

Now, because the system is close to being integrable, we can think of these six variables as naturally coming in three action-angle pairs. For any specific trajectory, the 'action' component will be more or less constant - and so, we expect that (very nearly) the whole system can be described by the motion of just three 'angle' variables. Each of these 'angle' variables will be controlled by a separate frequency - so three frequencies in total. Overall, then, we expect the system to be driven by just three frequencies (and harmonics of those frequencies). So, even though the motion of the Moon in this system is likely to be complicated, we expect its motion to be reduced to a series of harmonics of just three numbers. And, indeed, this is what we find. Quite remarkable in a way - but also absolutely true.

So, what are these three frequencies? Well, we don't know exactly what they are (and in a way it doesn't really matter what they represent) but we can guess what they might correspond to. Intuitively, we might expect that one frequency correspond to the orbital period of the Moon around the Earth. Another might obviously be related to the orbital period of the Earth/Moon system around the Sun. But what about the third frequency? What is that? The gravitational attraction of the Sun on the Moon exerts a torque on the Earth/Moon system and like a spinning top on a table, the orientation of the orbit of the Moon around the Earth precesses. The rate of precession of the Moon's orbit about the Earth's orbit is the third frequency. In summary, then, we expect to be able to reduce the motion of the Moon about Earth's centre to harmonics of these three fundamental frequency. And with a sufficiently accurate frequency analyser, one can.

Building the lunar ephemeris
So, how do we go about building the an ephemeris then? Really, there are just two steps. The first is to take a numerical integrator and then accurately integrate the trajectory of the Earth and the Moon over the period for which one wants to build an ephemeris (in this case one Julian century - 36500 days). The second is to take this integrated path and decompose it into a series of cosine terms that are just harmonics of the (in this case) just three fundamental frequencies.

The first step requires some assumption about initial conditions for the Earth and Moon (i.e., positions and velocities). Because I've squashed the system into two dimensions, choosing realistic initial conditions is a bit tricky. Mainly due to laziness on my part, I've chosen initial coordinates that lead to the Moon orbiting the Earth slightly faster than the real Moon (and with a substantially higher eccentricity). Having said that, my artificial 2D Earth/Moon system exhibits many of the same qualitative features of motion as the real thing and so the resulting ephemeris is semi-realistic.

The first step also requires an integrator. The integrator that I've used is just the fourth order symplectic integrator described else where in this forum. A time step was used to ensure that the Moon's position over the course of the Julian century was accurate to within about 1 km.

The second step is the analysis of frequencies. Here, I'm not going to give a blow-by-blow account suffice to say that the three fundamental frequencies were identified. The periods corresponding to these frequencies are as follows:

a. 23.8475 days. Given my slightly perverse initial conditions for the Moon, this is essentially the period of revolution of the Moon around the Earth, slightly shorter than one (real) lunar month.

b. 365.718 days. Again, given my slightly perverse initial conditions, this is the period of revolution of the Earth/Moon system around the Sun - very nearly one (real) calendar year.

c. 4069.06 days. In my artificial Earth/Moon system, this is the rate of precession of the Moon's orbit around the Earth due to the Sun's torque. It corresponds to about 11.14 years. In other word, the orientation of the Moon's orbit does one complete pirouette around the Earth every 11.14 years. How does this correspond to the real Moon? Well, because the real Moon moves in three dimensions, there are actually two periods to describe this motion. One is the period of precession of the lunar perigee. This is about 8.85 years. The other is period of precession of the orbital plane. And that is about 18.60 years. In two dimensions, these two collapse to being the same, so that one ends up with just one period which is somewhere between these two.

The final ephemeris, then consists of a series of cosine terms (amplitude, phase and frequency) for the x-coordinate of the Moon's position and the y-coordinate of the Moon's position. In turn, all of the frequencies of these cosine terms can be written as integer multiples of the three fundamental frequencies. So, for example, the first 8 lines of the ephemeris for the y-coordinate of the Moon's position are:

[math] \begin{array}{rrrrrrr} \text{Order} &n_1 &n_2 &n_3 &\text{Amplitude (AU)} & \text{Freq. (rads/day)}& \text{Phase (rads)} \\ 1& 1 & 0 & 0 &+0.0023354617 & +0.2634737904 & -1.570796448 \\ 2& 0 & 1 & 0 &+0.0003458840 & +0.0015441386 & -1.570796062 \\ 3& 0 & -1 & 2 &+0.0000566833 & +0.0328167165 & -1.570798382 \\ 4& 1 & 0 & -2 &+0.0000138571 & +0.2291129354 & -1.570795865 \end{array} [/math][math] \begin{array}{rrrrrrr} \text{Order} &n_1 &n_2 &n_3 &\text{Amplitude (AU)} & \text{Freq. (rads/day)}& \text{Phase (rads)} \\ 5& 2& -1& 0& +0.0001150510 & +0.5254034423 & +1.570795917 \\ 6& 2& 1& -2& +0.0000208703 & +0.4941308644 & +1.570799513 \\ 7& 3& -2& 0& +0.0000084408 & +0.7873330941 & -1.570797331 \\ 8& 3& 0& -2& +0.0000054780 & +0.7560605163 & -1.570792619 \\ \end{array} [/math]
In total, the full ephemeris consists of 41 such lines for the x-coordinate of position; and 41 for the y-coordinate of the Moon's position. These are contained in attached spreadsheet for those interested.

The first column contains an index similarly labelling the ephemeris term. The next three columns denote the relevant harmonic of the three fundamental frequency such that the term's frequency can be written as:

[MATH]\nu \to n_1 \omega _1+n_2 \omega _2+n_3 \omega _3[/MATH]
where

[MATH]\omega_1 = 0.263473790444[/MATH][MATH]\omega_2 = 0.001544138609[/MATH][MATH]\omega_3 = 0.017180427535[/MATH]
the final three terms contain the amplitude, the frequency and phase of the cosine term such that the overall coordinate value can be written as:

[MATH]x = \sum _{i=1}^{41} A_i \cos \left(t \nu _i+\phi _i\right)[/MATH]
And that about sums up that - now to repeat the exercise on a three dimensional version of the Sun/Earth/Moon system - if only to see how close we can get to the 'real' values - and then on to looking at the Saturn/Titan/Hyperion system.

P.S. I'll upload the ephemeris spreadsheet in a separate post, methinks.