Astro::Coord::ECI::Moon - Compute the position of the Moon.

 use Astro::Coord::ECI;
 use Astro::Coord::ECI::Moon;
 use Astro::Coord::ECI::Utils qw{deg2rad};

 # 1600 Pennsylvania Ave, Washington DC USA
 # latitude 38.899 N, longitude 77.038 W,
 # altitude 16.68 meters above sea level
 my $lat = deg2rad (38.899);    # Radians
 my $long = deg2rad (-77.038);  # Radians
 my $alt = 16.68 / 1000;        # Kilometers
 my $moon = Astro::Coord::ECI::Moon->new ();
 my $sta = Astro::Coord::ECI->
     universal (time ())->
     geodetic ($lat, $long, $alt);
 my ($time, $rise) = $sta->next_elevation ($moon);
 print "Moon @{[$rise ? 'rise' : 'set']} is ",
     scalar localtime $time, "\n";

package Astro::Coord::ECI::Moon;

use strict;
use warnings;

our $VERSION = '0.039';

use base qw{Astro::Coord::ECI};

use Astro::Coord::ECI::Sun;	# Need for phase of moon calc.
use Astro::Coord::ECI::Utils qw{:all};
use Carp;
use POSIX qw{floor strftime};


#	Load the periodic terms from the table.

my %terms;
{	# Begin local symbol block.
my $where;
while (<DATA>) {
    chomp;
    s/^\s+//;
    s/\s+$//;
    next unless $_;
    next if m/^#/;
    s/^-// and do {
	last if $_ eq 'end';
	$where = $terms{$_} ||= [];
	next;
	};
    s/_//g;
    push @$where, [split '\s+', $_];
    }
}	# End local symbol block.

my %static = (
    id => 'Moon',
    name => 'Moon',
    diameter => 3476,
    );

my $weaken = eval {
    require Scalar::Util;
    Scalar::Util->can('weaken');
    };
my $object;

our $Singleton = $weaken;

sub new {
    my ($class, @args) = @_;
    ref $class and $class = ref $class;
    if ( $Singleton && $weaken && _classisa( $class, __PACKAGE__ ) ) {
	if ($object) {
	    $object->set (@args) if @args;
	    return $object;
	} else {
	    my $self = $object = $class->SUPER::new (%static, @args);
	    $weaken->($object);
	    return $self;
	}
    } else {
	return $class->SUPER::new (%static, @args);
    }
}

my @quarters = ('New Moon', 'First quarter Moon', 'Full Moon',
    'Last quarter Moon');

sub almanac {
    my $self = shift;
    my $location = shift;
    embodies($location, 'Astro::Coord::ECI') or
	croak <<eod;
Error - The first argument of the almanac() method must be a member of
        the Astro::Coord::ECI class, or a subclass thereof.
eod

    my $start = shift || $self->universal;
    my $end = shift || $start + 86400;

    my @almanac;

    foreach (
	[$location, next_elevation => [$self, 0, 1], 'horizon',
		['Moon set', 'Moon rise']],
	[$location, next_meridian => [$self], 'transit',
		[undef, 'Moon transits meridian']],
	[$self, next_quarter => [], 'quarter',
		[@quarters]],
    ) {
	my ($obj, $method, $arg, $event, $descr) = @$_;
	$obj->universal ($start);
	while (1) {
	    my ($time, $which) = $obj->$method (@$arg);
	    last if $time >= $end;
	    push @almanac, [$time, $event, $which, $descr->[$which]]
		if $descr->[$which];
	}
    }
    return (sort {$a->[0] <=> $b->[0]} @almanac);
}

sub almanac_hash {
    return map {
	body => $_[0],
	station => $_[1],
	time => $_->[0],
	almanac => {
	    event => $_->[1],
	    detail => $_->[2],
	    description => $_->[3],
	}
    }, almanac(@_);
}

use constant QUARTER_INC => 86400 * 6;	# 6 days.

sub next_quarter {
    my ($self, $quarter) = @_;
    $quarter = (defined $quarter ? $quarter :
	floor ($self->phase () / PIOVER2) + 1) % 4;
    my $begin;
    while (floor ($self->phase () / PIOVER2) == $quarter) {
	$begin = $self->dynamical;
	$self->dynamical ($begin + QUARTER_INC);	# Bump 6 days.
    }
    while (floor ($self->phase () / PIOVER2) != $quarter) {
	$begin = $self->dynamical;
	$self->dynamical ($begin + QUARTER_INC);	# Bump 6 days.
    }
    my $end = $self->dynamical ();

    while ($end - $begin > 1) {
	my $mid = floor (($begin + $end) / 2);
	my $qq = floor ($self->dynamical ($mid)->phase () / PIOVER2);
	($begin, $end) = $qq == $quarter ?
	    ($begin, $mid) : ($mid, $end);
    }

    $self->dynamical ($end);

    return wantarray ?
	($self->universal, $quarter, $quarters[$quarter]) :
	$self->universal;
}

sub next_quarter_hash {
    my ($self, @args) = @_;
    my ($time, $quarter, $desc) = $self->next_quarter(@args);
    my %hash = (
	body => $self,
	almanac => {
	    event => 'quarter',
	    detail => $quarter,
	    description => $desc,
	},
	time => $time,
    );
    return wantarray ? %hash : \%hash;
}

sub period {return 2360591.5968}	# 27.321662 * 86400

sub phase {
    my ($self, $time) = @_;

    (ref $self || $time) or croak <<eod;
Error - You must specify a time if you call phase() as a class method.
eod

    $self = $self->new () unless ref $self;

    $self->universal ($time) if $time;

    my $sun = Astro::Coord::ECI::Sun->universal ($self->universal);

    my (undef, $longs) = $sun->ecliptic ();
    my (undef, $longm) = $self->ecliptic ();

    my $phase = mod2pi ($longm - $longs);
    return wantarray ? ($phase, (1 + cos ($self->PI - $phase)) / 2) : $phase;
}

sub time_set {
    my $self = shift;

    my $time = $self->dynamical;

    my $T = jcent2000 ($time);			# Meeus (22.1)

#	Moon's mean longitude.

    my $Lprime = mod2pi (deg2rad (			# Meeus (47.1)
	(((- ($T / 65_194_000) +
	1 / 538_841) * $T - 0.0015786) * $T +
	481267.88123421) * $T + 218.3164477));

#	Moon's mean elongation.

    my $D = mod2pi (deg2rad (((($T / 113_065_000 +	# Meeus (47.2)
	1 / 545_868) * $T - 0.0018819) * $T +
	445267.1114034) * $T + 297.8501921));

#	Sun's mean anomaly.

    my $M = mod2pi (deg2rad ((($T / 24_490_000 -		# Meeus (47.3)
	0.000_1536) * $T + 35999.050_2909) * $T +
	357.5291092));

#	Moon's mean anomaly.

    my $Mprime = mod2pi (deg2rad ((((- $T / 14_712_000 +	# Meeus (47.4)
	1 / 69_699) * $T + 0.008_7414) * $T +
	477198.867_5055) * $T + 134.963_3964));

#	Moon's argument of latitude (mean distance
#	from ascending node).

    my $F = mod2pi (deg2rad (((($T / 863_310_000 -	# Meeus (47.5)
	1 / 3_526_000) * $T - 0.003_6539) * $T +
	483202.017_5233) * $T + 93.272_0950));

#	Eccentricity correction factor.

    my $E = (- 0.000_0074 * $T - 0.002_516) * $T + 1;	# Meeus (47.6)
    my @efac = (1, $E, $E * $E);

#	Compute "further arguments".

    my $A1 = mod2pi (deg2rad (131.849 * $T + 119.75));		# Venus
    my $A2 = mod2pi (deg2rad (479264.290 * $T + 53.09));	# Jupiter
    my $A3 = mod2pi (deg2rad (481266.484 * $T + 313.45));	# undocumented

#	Compute periodic terms for longitude (sigma l) and
#	distance (sigma r).

    my ($sigmal, $sigmar) = (0, 0);
    foreach (@{$terms{lr}}) {
	my ($mulD, $mulM, $mulMprime, $mulF, $sincof, $coscof) = @$_;
	if ($mulM) {
	    my $corr = $efac[abs $mulM] || confess <<eod;
Programming error - M multiple greater than 2.
eod
	    $sincof *= $corr;
	    $coscof *= $corr;
	}
	my $arg = $D * $mulD + $M * $mulM + $Mprime * $mulMprime +
	    $F * $mulF;
	$sigmal += $sincof * sin ($arg);
	$sigmar += $coscof * cos ($arg);
    }

#	Compute periodic terms for latitude (sigma b).

    my $sigmab = 0;
    foreach (@{$terms{b}}) {
	my ($mulD, $mulM, $mulMprime, $mulF, $sincof) = @$_;
	if ($mulM) {
	    my $corr = $efac[abs $mulM] || confess <<eod;
Programming error - M multiple greater than 2.
eod
	    $sincof *= $corr;
	}
	my $arg = $D * $mulD + $M * $mulM + $Mprime * $mulMprime +
	    $F * $mulF;
	$sigmab += $sincof * sin ($arg);
    }

#	Add other terms

    $sigmal += 3958 * sin ($A1) + 1962 * sin ($Lprime - $F) +
	318 * sin ($A2);
    $sigmab += - 2235 * sin ($Lprime) + 382 * sin ($A3) +
	175 * sin ($A1 - $F) + 175 * sin ($A1 + $F) +
	127 * sin ($Lprime - $Mprime) - 115 * sin ($Lprime + $Mprime);

#	Coordinates of Moon (finally!)

    my $lambda = mod2pi( deg2rad ($sigmal / 1_000_000) + $Lprime );
    my $beta = deg2rad ($sigmab / 1_000_000);
    my $delta = $sigmar / 1000 + 385_000.56;

#	Correct longitude for nutation (from Chapter 22, pg 144).

    $lambda += nutation_in_longitude ($time);

    $self->ecliptic ($beta, $lambda, $delta);
    ## $self->set (equinox_dynamical => $time);
    $self->equinox_dynamical ($time);
    return $self;
}

# The moon is normally positioned in inertial coordinates.

sub _initial_inertial { return 1 }

1;

__DATA__
# Source: Jean Meeus' "Astronomical Algorithms", Table 47.A pg. 339-340.
-lr
0  0  1  0	6_288_774	-20_905_355
2  0 -1  0	1_274_027	 -3_699_111
2  0  0  0	  658_314	 -2_955_968
0  0  2  0	  213_618	   -569_925
0  1  0  0	 -185_116	     48_888
0  0  0  2	 -114_332	     -3_149
2  0 -2  0	   58_793	    246_158
2 -1 -1  0	   57_066	   -152_138
2  0  1  0	   53_322	   -170_733
2 -1  0  0	   45_758	   -204_586
0  1 -1  0	  -40_923	   -129_620
1  0  0  0	  -34_720	    108_743
0  1  1  0	  -30_383	    104_755
2  0  0 -2	   15_327	     10_321
0  0  1  2	  -12_528	          0
0  0  1 -2	   10_980	     79_661
4  0 -1  0	   10_675	    -34_782
0  0  3  0	   10_034	    -23_210
4  0 -2  0	    8_548	    -21_636
2  1 -1  0	   -7_888	     24_208
2  1  0  0	   -6_766	     30_824
1  0 -1  0	   -5_163	     -8_379
1  1  0  0	    4_987	    -16_675
2 -1  1  0	    4_036	    -12_831
2  0  2  0	    3_994	    -10_445
4  0  0  0	    3_861	    -11_650
2  0 -3  0	    3_665	     14_403
0  1 -2  0	   -2_689	     -7_003
2  0 -1  2	   -2_602	          0
2 -1 -2  0	    2_390	     10_056
1  0  1  0	   -2_348	      6_322
2 -2  0  0	    2_236	     -9_884
0  1  2  0	   -2_120	      5_751
0  2  0  0	   -2_069	          0
2 -2 -1  0	    2_048	     -4_950
2  0  1 -2	   -1_773	      4_130
2  0  0  2	   -1_595	          0
4 -1 -1  0	    1_215	     -3_958
0  0  2  2	   -1_110	          0
3  0 -1  0	     -892	      3_258
2  1  1  0	     -810	      2_616
4 -1 -2  0	      759	     -1_897
0  2 -1  0	     -713	     -2_117
2  2 -1  0	     -700	      2_354
2  1 -2  0	      691	          0
2 -1  0 -2	      596	          0
4  0  1  0	      549	     -1_423
0  0  4  0	      537	     -1_117
4 -1  0  0	      520	     -1_571
1  0 -2  0	     -487	     -1_739
2  1  0 -2	     -399	          0
0  0  2 -2	     -381	     -4_421
1  1  1  0	      351	          0
3  0 -2  0	     -340	          0
4  0 -3  0	      330	          0
2 -1  2  0	      327	          0
0  2  1  0	     -323	      1_165
1  1 -1  0	      299	          0
2  0  3  0	      294	          0
2  0 -1 -2	        0	      8_752
# Source: Jean Meeus' "Astronomical Algorithms", Table 47.B pg. 341.
-b
#>>> I have an error in the sigma b result versus the book
#>>> solution when working Meeus' example 47.a on page 342.
#>>> The book gives -3229126, and I get -3228669.6
0  0  0  1	5_128_122
0  0  1  1	  280_602
0  0  1 -1	  277_693
2  0  0 -1	  173_237
2  0 -1  1	   55_413
2  0 -1 -1	   46_271
2  0  0  1	   32_573
0  0  2  1	   17_198
2  0  1 -1	    9_266
0  0  2 -1	    8_822
2 -1  0 -1	    8_216
2  0 -2 -1	    4_324
2  0  1  1	    4_200
2  1  0 -1	   -3_359
2 -1 -1  1	    2_463
2 -1  0  1	    2_211
2 -1 -1 -1	    2_065
0  1 -1 -1	   -1_870
4  0 -1 -1	    1_828
0  1  0  1	   -1_794
0  0  0  3	   -1_749
0  1 -1  1	   -1_565
1  0  0  1	   -1_491
0  1  1  1	   -1_475
0  1  1 -1	   -1_410
0  1  0 -1	   -1_344
1  0  0 -1	   -1_335
0  0  3  1	    1_107
4  0  0 -1	    1_021
4  0 -1  1	      833
0  0  1 -3	      777
4  0 -2  1	      671
2  0  0 -3	      607
2  0  2 -1	      596
2 -1  1 -1	      491
2  0 -2  1	     -451
0  0  3 -1	      439
2  0  2  1	      422
2  0 -3 -1	      421
2  1 -1  1	     -366
2  1  0  1	     -351
4  0  0  1	      331
2 -1  1  1	      315
2 -2  0 -1	      302
0  0  1  3	     -283
2  1  1 -1	     -229
1  1  0 -1	      223
1  1  0  1	      223
0  1 -2 -1	     -220
2  1 -1 -1	     -220
1  0  1  1	     -185
2 -1 -2 -1	      181
0  1  2  1	     -177
4  0 -2 -1	      176
4 -1 -1 -1	      166
1  0  1 -1	     -164
4  0  1 -1	      132
1  0 -1 -1	     -119
4 -1  0 -1	      115
2 -2  0  1	      107