# 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 elongation
   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);
}