Saturday, August 10, 2013

"Aim it at... The Moon!"

So I thought, "keplerian elements are pretty straightforward once you get past the jargon", but there's different jargon for comets as asteroids, so that's still a work in progress. And I assumed that plotting the position of the moon would be a similar exercise.

Oh boy. Utterly not.

Let's just say that the best and easiest way is probably to ping NASA and just ask their big ephemeris computer, called HORIZON.

I did it the masochistic way, which is to port the ELP/MPP002 code out of GAL (the 'General Astrodynamics Library') and into Javascript. The amazing thing is that it seems to have worked.

See how all the lines are off from each other? Like it's a couple of rounds of a strange attractor? Well, that's actually correct. The "orbit of the moon" is chaotic. Not a neat hoop. Practically everything in the solar system perturbs it. The size of the ocean (and hence tide) it's currently over affects it. It's orbit shrinks when the Earth is farther from the sun, and expands again when we're closest.

And since Earth is the Moon's dancing partner, anything it does, we do to, but considerably less so. Plotting Earth's true position (compared to the stable "earth-moon barycenter" which is what most programs approximate down to) means adding a suitably scaled opposing vector.

The equation which computes this orbit is literally the largest I have ever dealt with. It is a polynomial with upwards of 20,000 terms. The code which runs is barely 50 lines long, but the 'coefficients' file it accesses is 10Mb of data.

This 'Lunar Equation' is not a physical computation of any kind. The coefficients represent a 'best fit' curve to the JPL ephemeris computations. (which are physically based) The JPL supercomputer crunched the numbers to find the position over time (ignore that part) and then the ELP computer "compressed" that path down to a best-fit curve. The code then "uncompresses" the coordinates for a given point in time.

Even that 10Mb of data is not the true 'best fit' curve: but, they discovered that adding any more terms actually degrades the computation through rounding errors more than the accuracy gained. That's how big it is.

Once again, I'm surprised at the numerical stability I'm getting out of Javascript's native floats. It's agreeing with the test cases down to the meter, over hundreds of years. Again, it would be good to re-implement the native math with a BigDecimal class, but the performance hit will be intense. Applying those 20,000 trig coefficients means I can only compute about 250 points-per-second on a beefy modern machine. That year's worth of orbit took about four seconds to calculate.

So, I can compute moon positions client-side in the browser now. At the cost of downloading 10Mb of javascript code and intense processing. That's a tradeoff that only makes sense in some circumstances. A pre-computed table of a year's moon positions, three times a day, would only take 3,000 numbers, and would be instant-lookup.

Then again, the GAL table contains _two_ lunar models. It might be possible to split those out. 5Mb (on demand) starts being a very manageable quantity. And the code can surely be optimized.

No comments:

Post a Comment