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