Here's a relatively short program that calculates the fraction of illumination of the moon at a given time. This is relevant right now because it's a function of the angle between the sun and the moon. When the angle is small enough, there's a solar eclipse, as we will get a chance to see tomorrow. Unfortunately, the moon wobbles around in the sky too much for a simple program like this to cope with, so it can't produce high-accuracy eclipse predictions, but it might be interesting to some people. Visit
# Calculate illumination of the moon, and approximate eclipses.
# Source: "Numerical expressions for precession formulae and mean
# elements for the Moon and the planets" by J.L. Simon et.al.,
# Astronomy and Astrophysics vol. 282, no. 2, p.663-683 (1994)
use strict;
use warnings;
use POSIX qw( floor );
use constant PI => 4*atan2(1,1);
{
# Time in Julian centuries since Jan 1 2000 12:00 GMT
my $T = (time() - 946728000) / (36525 * 24 * 60 * 60);
my ($xe, $ye, $ze) = earth($T);
my ($xm, $ym, $zm) = moon($T);
my $illum = 0.5 + 0.5 * ($xe*$xm + $ye*$ym + $ze*$zm)
/ sqrt($xe*$xe + $ye*$ye + $ze*$ze)
/ sqrt($xm*$xm + $ym*$ym + $zm*$zm);
print "illumination = $illum\n";
print "lunar eclipse\n" if $illum > 0.999924;
print "solar eclipse\n" if $illum < 0.000076;
# Moon position is only accurate to about half a lunar radius, not
# good enough to determine if there's a partial or total eclipse.
}
sub earth {
# Approximate heliocentric coordinates of the earth-moon barycenter
+ in km
my ($T) = @_;
my $a = 149598022.96
+ 2440*cos(5.11567+575.32983*$T)
+ 2376*cos(3.43546+786.05399*$T)
+ 1377*cos(0.83152+1150.65965*$T)
+ 817*cos(1.71787+393.00902*$T)
+ 523*cos(0.41241+522.37014*$T)
+ 711*cos(2.61200+588.48885*$T)
+ 506*cos(5.30859+550.73755*$T)
+ 368*cos(2.03444+1179.06301*$T);
my $e = 0.0167086342 - 0.00004203654*$T;
my $L = 1.7534704595 + 628.331965406500*$T
+ 0.0000342*cos(3.45408+0.35954*$T)
+ 0.0000350*cos(3.54386+575.32983*$T)
+ 0.0000270*cos(1.86793+786.05399*$T)
+ 0.0000235*cos(0.14973+393.00902*$T)
+ 0.0000127*cos(4.29097+52.95968*$T)
+ 0.0000131*cos(5.54640+1150.65965*$T)
+ 0.0000125*cos(5.16887+157.72853*$T)
+ 0.0000090*cos(4.23879+2.62461*$T);
my $p = 1.7965956473 + 0.03001023491*$T;
return orbit($a, $e, 0, $L, $p, 0);
}
sub moon {
# Approximate geocentric coordinates of the moon in km
my ($T) = @_;
my $D = 5.1984667410 + 7771.3771468120*$T; # moon mean solar elonga
+tion
my $F = 1.6279052334 + 8433.4661581308*$T; # moon mean argument of
+latitude
my $lm = 2.3555558983 + 8328.6914269554*$T; # moon mean anomaly, l
my $ls = 6.2400601269 + 628.3019551714*$T; # sun mean anomaly, l'
my $a = 383397.7916 # semimajor axis
+ 3400.4*cos(2*$D)
- 635.6*cos(2*$D-$lm)
- 235.6*cos($lm)
+ 218.1*cos(2*$D-$ls)
+ 181.0*cos(2*$D+$lm);
my $e = 0.055545526 # eccentricity
+ 0.014216*cos(2*$D-$lm)
+ 0.008551*cos(2*$D-2*$lm)
- 0.001383*cos($lm)
+ 0.001356*cos(2*$D+$lm)
- 0.001147*cos(4*$D-3*$lm)
- 0.000914*cos(4*$D-2*$lm)
+ 0.000869*cos(2*$D-$ls-$lm)
- 0.000627*cos(2*$D)
- 0.000394*cos(4*$D-4*$lm)
+ 0.000282*cos(2*$D-$ls-2*$lm)
- 0.000279*cos($D-$lm)
- 0.000236*cos(2*$lm)
+ 0.000231*cos(4*$D)
+ 0.000229*cos(6*$D-4*$lm)
- 0.000201*cos(2*$lm-2*$F);
my $I = 0.0900012160 # inclination
+ 0.00235746*cos(2*$D-2*$F)
+ 0.00012474*cos(2*$lm-2*$F)
+ 0.00009682*cos(2*$D-$ls-2*$F)
- 0.00001270*cos(2*$D);
my $L # mean longitude, lambda
= 3.8103444305 + 8399.70911352222*$T
- 0.0161584*sin(2*$D)
+ 0.0058052*sin(2*$D-$lm)
- 0.0032119*sin($ls)
+ 0.0019213*sin($lm)
- 0.0010569*sin(2*$D-$ls);
my $p # longitude of perigee, omega bar
= 1.4547885324 + 71.0176865667*$T
- 0.269600*sin(2*$D-$lm)
- 0.168284*sin(2*$D-2*$lm)
- 0.047473*sin($lm)
+ 0.045500*sin(4*$D-3*$lm)
+ 0.036385*sin(4*$D-2*$lm)
+ 0.025782*sin(2*$D+$lm)
+ 0.016891*sin(4*$D-4*$lm)
- 0.016566*sin(2*$D-$ls-$lm)
- 0.012266*sin(6*$D-4*$lm)
- 0.011519*sin(2*$D)
- 0.010060*sin(2*$D-3*$lm)
- 0.009129*sin(2*$lm)
- 0.008416*sin(6*$D-5*$lm)
+ 0.007883*sin($ls)
- 0.006642*sin(6*$D-3*$lm);
my $N # longitude of ascending node, Omega
= 2.1824391971 - 33.7570446083*$T
- 0.026141*sin(2*$D-2*$F)
- 0.002618*sin($ls)
- 0.002138*sin(2*$D)
+ 0.002051*sin(2*$F)
- 0.001396*sin(2*$lm-2*$F);
return orbit($a, $e, $I, $L, $p, $N);
}
sub orbit {
my ($a, $e, $I, $L, $p, $N) = @_;
my $w = $p - $N;
my $M = $L - $p;
$M -= (2*PI) * floor($M / (2*PI) + 0.5);
# Kepler's equation
my $E = $M + $e * sin($M);
while (1) {
my $dM = $M - $E + $e * sin($E);
my $dE = $dM / (1 - $e * cos($E));
$E += $dE;
last if abs($dE) < 1e-6;
}
my $x = $a * (cos($E) - $e);
my $y = $a * sqrt(1 - $e*$e) * sin($E);
my ($c, $s) = (cos($w), sin($w));
($x, $y) = ($c*$x - $s*$y, $s*$x + $c*$y);
my $z = $y * sin($I);
$y *= cos($I);
($c, $s) = (cos($N), sin($N));
($x, $y) = ($c*$x - $s*$y, $s*$x + $c*$y);
return ($x, $y, $z);
}