http://www.perlmonks.org?node_id=1197712

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 JPL for more information.
# 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); }

Replies are listed 'Best First'.
Re: moon illumination and eclipses
by Your Mother (Chancellor) on Aug 20, 2017 at 20:40 UTC

    Super cool. Related: Astro::MoonPhase. I have been using an alias/command with it for ages

    moo@cow[6]~>which moon aliased to perl -MAstro::MoonPhase -e "printf qq{The Moon is %.3f%% fu +ll\n},100*(phase)[1]" moo@cow[9]~>moon The Moon is 1.109% full

    Counting down to zero for tomorrow. :P

      I should maybe point out that Astro::MoonPhase computes an idealized illumination that's exactly 1 at the full moon. That's not exactly accurate, which is where the connection with eclipses comes in.
Re: moon illumination and eclipses
by no_slogan (Deacon) on Aug 20, 2017 at 22:58 UTC
    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.