Also, you can get breathtakingly high accuracy (about 1 mm) on the position of the moon by using JPL's high-accuracy ephemeris. The catch is that the fit polynomials are only good for a few days at a time, so you need a big database to cover any significant length of time.
{
# Time in Julian days since Jan 1 2000 12:00 GMT
my $T = (time() - 946728000) / (24 * 60 * 60);
my ($xm, $ym, $zm) = calc_position($T, 6439.5, 6443.5,
[-0.2576477260382208e6, -0.1244143318652203e6, 0.157149496430333
+8e5,
0.1128910971393731e4, -0.1026042752851206e3, -0.3344243625567091
+e1,
0.5108570281514203, 0.7273850253892471e-2, -0.2976495487727882e-
+2,
-0.1508820799631562e-4, 0.1988387208061850e-4, -0.99448295755443
+34e-7,
-0.1420192252858349e-6],
[0.2161347296153234e6, -0.1235954423872278e6, -0.133633064276375
+8e5,
0.1347907704333395e4, 0.6144497680180469e2, -0.6876023550099470e
+1,
-0.1365228336729124, 0.3626374756466881e-1, 0.2847434884414048e-
+3,
-0.2277001670510105e-3, 0.2675604082417520e-6, 0.158232148064080
+5e-5,
-0.1939924674109133e-7],
[0.8988317387937439e5, -0.3579858449186091e5, -0.554716617427666
+7e4,
0.4032955946240293e3, 0.2725379907222771e2, -0.2195153465571546e
+1,
-0.7682906296761302e-1, 0.1217607039333434e-1, 0.270238930850790
+9e-3,
-0.7822887028838594e-4, -0.1051459145053017e-5, 0.55540570861312
+47e-6,
0.1433643809289891e-8]);
}
sub calc_position {
# Evaluate Chebyshev polynomials to get x,y,z position
my ($T, $Tstart, $Tend, $xpoly, $ypoly, $zpoly) = @_;
die 'time out of range' if $T < $Tstart || $T > $Tend;
my $t = ($T - $Tstart) / ($Tend - $Tstart) * 2 - 1;
my $x = $xpoly->[0] + $xpoly->[1] * $t;
my $y = $ypoly->[0] + $ypoly->[1] * $t;
my $z = $zpoly->[0] + $zpoly->[1] * $t;
my $T0 = 1;
my $T1 = $t;
$t *= 2;
for my $i (2 .. $#$xpoly) {
($T1, $T0) = ($t*$T1 - $T0, $T1);
$x += $T1 * $xpoly->[$i];
$y += $T1 * $ypoly->[$i];
$z += $T1 * $zpoly->[$i];
}
return ($x, $y, $z);
}
