There's more than one way to do things PerlMonks

moon illumination and eclipses

by no_slogan (Deacon)
 on Aug 20, 2017 at 20:00 UTC ( #1197712=CUFP: print w/replies, xml ) Need Help??

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.

Create A New User
Node Status?
node history
Node Type: CUFP [id://1197712]
Front-paged by stevieb
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others surveying the Monastery: (9)
As of 2017-10-19 07:53 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My fridge is mostly full of:

Results (252 votes). Check out past polls.

Notices?