package Astro::MoonPhase;

use strict;
use vars qw($VERSION @ISA @EXPORT @EXPORT_OK);

require Exporter;

@ISA = qw(Exporter);
@EXPORT = qw(phase phasehunt phaselist);
$VERSION = '0.60';

use vars qw (
			$Epoch
			$Elonge $Elongp $Eccent	$Sunsmax $Sunangsiz
			$Mmlong $Mmlongp $Mlnode $Minc $Mecc $Mangsiz $Msmax $Mparallax $Synmonth
			$Pi
			);

# Astronomical constants.

$Epoch					= 2444238.5;		# 1980 January 0.0

# Constants defining the Sun's apparent orbit.

$Elonge					= 278.833540;		# ecliptic longitude of the Sun at epoch 1980.0
$Elongp					= 282.596403;		# ecliptic longitude of the Sun at perigee
$Eccent					= 0.016718;		# eccentricity of Earth's orbit
$Sunsmax				= 1.495985e8;		# semi-major axis of Earth's orbit, km
$Sunangsiz				= 0.533128;		# sun's angular size, degrees, at semi-major axis distance

# Elements of the Moon's orbit, epoch 1980.0.

$Mmlong					= 64.975464;		# moon's mean longitude at the epoch
$Mmlongp				= 349.383063;		# mean longitude of the perigee at the epoch
$Mlnode					= 151.950429;		# mean longitude of the node at the epoch
$Minc					= 5.145396;		# inclination of the Moon's orbit
$Mecc					= 0.054900;		# eccentricity of the Moon's orbit
$Mangsiz				= 0.5181;		# moon's angular size at distance a from Earth
$Msmax					= 384401.0;		# semi-major axis of Moon's orbit in km
$Mparallax				= 0.9507;		# parallax at distance a from Earth
$Synmonth				= 29.53058868;		# synodic month (new Moon to new Moon)

# Properties of the Earth.

$Pi					= 3.14159265358979323846;	# assume not near black hole nor in Tennessee

# Handy mathematical functions.

sub sgn		{ return (($_[0] < 0) ? -1 : ($_[0] > 0 ? 1 : 0)); } 	# extract sign
sub fixangle	{ return ($_[0] - 360.0 * (floor($_[0] / 360.0))); }	# fix angle
sub torad	{ return ($_[0] * ($Pi / 180.0)); }						# deg->rad
sub todeg	{ return ($_[0] * (180.0 / $Pi)); }						# rad->deg
sub dsin	{ return (sin(torad($_[0]))); }						# sin from deg
sub dcos	{ return (cos(torad($_[0]))); }						# cos from deg

sub tan		{ return sin($_[0])/cos($_[0]); }
sub asin	{ return ($_[0]<-1 or $_[0]>1) ? undef : atan2($_[0],sqrt(1-$_[0]*$_[0])); }
sub atan {
    if		($_[0]==0)	{ return 0; }
	elsif	($_[0]>0)	{ return atan2(sqrt(1+$_[0]*$_[0]),sqrt(1+1/($_[0]*$_[0]))); }
	else 				{ return -atan2(sqrt(1+$_[0]*$_[0]),sqrt(1+1/($_[0]*$_[0]))); }
}

sub floor {
  my $val   = shift;
  my $neg   = $val < 0;
  my $asint = int($val);
  my $exact = $val == $asint;

  return ($exact ? $asint : $neg ? $asint - 1 : $asint);
}



# jtime - convert internal date and time to astronomical Julian
# time (i.e. Julian date plus day fraction)


sub jtime {

        my $t = shift;
        my ($julian);
        $julian = ($t / 86400) + 2440587.5;	# (seconds /(seconds per day)) + julian date of epoch
        return ($julian);

}


# jdaytosecs - convert Julian date to a UNIX epoch

sub jdaytosecs {
  my $jday = shift;
  my $stamp;

  $stamp = ($jday - 2440587.5)*86400;   # (juliandate - jdate of unix epoch)*(seconds per julian day)
  return($stamp);

}



# jyear - convert Julian date to year, month, day, which are
# returned via integer pointers to integers
sub jyear {

	my ($td, $yy, $mm, $dd) = @_;
	my ($z, $f, $a, $alpha, $b, $c, $d, $e);

	$td += 0.5;				# astronomical to civil
	$z = floor($td);
	$f = $td - $z;

	if ($z < 2299161.0) {
		$a = $z;
	} else {
		$alpha = floor(($z - 1867216.25) / 36524.25);
		$a = $z + 1 + $alpha - floor($alpha / 4);
	}


	$b = $a + 1524;
	$c = floor(($b - 122.1) / 365.25);
	$d = floor(365.25 * $c);
	$e = floor(($b - $d) / 30.6001);

	$$dd = $b - $d - floor(30.6001 * $e) + $f;
	$$mm = $e < 14 ? $e - 1 : $e - 13;
	$$yy = $$mm > 2 ? $c - 4716 : $c - 4715;

}

##  meanphase  --  Calculates  time  of  the mean new Moon for a given
##                 base date.  This argument K to this function is the
##                 precomputed synodic month index, given by:
##
##                        K = (year - 1900) * 12.3685
##
##                 where year is expressed as a year and fractional year.


sub meanphase {
  my ($sdate, $k) = @_;
  my ($t, $t2, $t3, $nt1);

  ## Time in Julian centuries from 1900 January 0.5
  $t = ($sdate - 2415020.0) / 36525;
  $t2 = $t * $t;                       ## Square for frequent use
  $t3 = $t2 * $t;                      ## Cube for frequent use

  $nt1 = 2415020.75933 + $Synmonth * $k
         + 0.0001178 * $t2
         - 0.000000155 * $t3
         + 0.00033 * dsin(166.56 + 132.87 * $t - 0.009173 * $t2);

  return ($nt1);
}


# truephase - given a K value used to determine the mean phase of the
# new moon, and a phase selector (0.0, 0.25, 0.5, 0.75),
# obtain the true, corrected phase time

sub truephase {
	my ($k, $phase) = @_;
	my ($t, $t2, $t3, $pt, $m, $mprime, $f);
	my $apcor = 0;

	$k += $phase;				# add phase to new moon time
	$t = $k / 1236.85;			# time in Julian centuries from
								# 1900 January 0.5
	$t2 = $t * $t;				# square for frequent use
	$t3 = $t2 * $t;				# cube for frequent use

	# mean time of phase */
	$pt = 2415020.75933
	 + $Synmonth * $k
	 + 0.0001178 * $t2
	 - 0.000000155 * $t3
	 + 0.00033 * dsin(166.56 + 132.87 * $t - 0.009173 * $t2);

	# Sun's mean anomaly
	$m = 359.2242
	+ 29.10535608 * $k
	- 0.0000333 * $t2
	- 0.00000347 * $t3;

	# Moon's mean anomaly
	$mprime = 306.0253
	+ 385.81691806 * $k
	+ 0.0107306 * $t2
	+ 0.00001236 * $t3;

	# Moon's argument of latitude
	$f = 21.2964
	+ 390.67050646 * $k
	- 0.0016528 * $t2
	- 0.00000239 * $t3;

	if (($phase < 0.01) || (abs($phase - 0.5) < 0.01)) {
		# Corrections for New and Full Moon.

		$pt += (0.1734 - 0.000393 * $t) * dsin($m)
		 + 0.0021 * dsin(2 * $m)
		 - 0.4068 * dsin($mprime)
		 + 0.0161 * dsin(2 * $mprime)
		 - 0.0004 * dsin(3 * $mprime)
		 + 0.0104 * dsin(2 * $f)
		 - 0.0051 * dsin($m + $mprime)
		 - 0.0074 * dsin($m - $mprime)
		 + 0.0004 * dsin(2 * $f + $m)
		 - 0.0004 * dsin(2 * $f - $m)
		 - 0.0006 * dsin(2 * $f + $mprime)
		 + 0.0010 * dsin(2 * $f - $mprime)
		 + 0.0005 * dsin($m + 2 * $mprime);
	$apcor = 1;
	}
	elsif ((abs($phase - 0.25) < 0.01 || (abs($phase - 0.75) < 0.01))) {
		$pt += (0.1721 - 0.0004 * $t) * dsin($m)
		 + 0.0021 * dsin(2 * $m)
		 - 0.6280 * dsin($mprime)
		 + 0.0089 * dsin(2 * $mprime)
		 - 0.0004 * dsin(3 * $mprime)
		 + 0.0079 * dsin(2 * $f)
		 - 0.0119 * dsin($m + $mprime)
		 - 0.0047 * dsin($m - $mprime)
		 + 0.0003 * dsin(2 * $f + $m)
		 - 0.0004 * dsin(2 * $f - $m)
		 - 0.0006 * dsin(2 * $f + $mprime)
		 + 0.0021 * dsin(2 * $f - $mprime)
		 + 0.0003 * dsin($m + 2 * $mprime)
		 + 0.0004 * dsin($m - 2 * $mprime)
		 - 0.0003 * dsin(2 * $m + $mprime);
		if ($phase < 0.5) {
			# First quarter correction.
			$pt += 0.0028 - 0.0004 * dcos($m) + 0.0003 * dcos($mprime);
		}
		else {
			# Last quarter correction.
			$pt += -0.0028 + 0.0004 * dcos($m) - 0.0003 * dcos($mprime);
		}
		$apcor = 1;
	}
	if (!$apcor) {
		die "truephase() called with invalid phase selector ($phase).\n";
	}
	return ($pt);
}

# phasehunt - find time of phases of the moon which surround the current
# date.  Five phases are found, starting and ending with the
# new moons which bound the current lunation

sub phasehunt {
        my $sdate = jtime(shift || time());
        my ($adate, $k1, $k2, $nt1, $nt2);
        my ($yy, $mm, $dd);

        $adate = $sdate - 45;

	jyear($adate, \$yy, \$mm, \$dd);
	$k1 = floor(($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685);

	$adate = $nt1 = meanphase($adate,  $k1);

        while (1) {
                $adate += $Synmonth;
		$k2 = $k1 + 1;
                $nt2 = meanphase($adate, $k2);
                if (($nt1 <= $sdate) && ($nt2 > $sdate)) {
                        last;
                }
                $nt1 = $nt2;
                $k1 = $k2;

        }



        return  (
                        jdaytosecs(truephase($k1, 0.0)),
                        jdaytosecs(truephase($k1, 0.25)),
                        jdaytosecs(truephase($k1, 0.5)),
                        jdaytosecs(truephase($k1, 0.75)),
                        jdaytosecs(truephase($k2, 0.0))
                        );
}



# phaselist - find time of phases of the moon between two dates
# times (in & out) are seconds_since_1970

sub phaselist
{
  my ($sdate, $edate) = map { jtime($_) } @_;

  my (@phases, $d, $k, $yy, $mm);

  jyear($sdate, \$yy, \$mm, \$d);
  $k = floor(($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685) - 2;

  while (1) {
    ++$k;
    for my $phase (0.0, 0.25, 0.5, 0.75) {
      $d = truephase($k, $phase);

      return @phases if $d >= $edate;

      if ($d >= $sdate) {
        push @phases, int(4 * $phase) unless @phases;
        push @phases, jdaytosecs($d);
      } # end if date should be listed
    } # end for each $phase
  } # end while 1
} # end phaselist



# kepler - solve the equation of Kepler

sub kepler {
	my ($m, $ecc) = @_;
	my ($e, $delta);
	my $EPSILON = 1e-6;

	$m = torad($m);
	$e = $m;
	do {
		$delta = $e - $ecc * sin($e) - $m;
		$e -= $delta / (1 - $ecc * cos($e));
	} while (abs($delta) > $EPSILON);
	return ($e);
}



# phase - calculate phase of moon as a fraction:
#
# The argument is the time for which the phase is requested,
# expressed as a Julian date and fraction.  Returns the terminator
# phase angle as a percentage of a full circle (i.e., 0 to 1),
# and stores into pointer arguments the illuminated fraction of
# the Moon's disc, the Moon's age in days and fraction, the
# distance of the Moon from the centre of the Earth, and the
# angular diameter subtended by the Moon as seen by an observer
# at the centre of the Earth.

sub phase {
	my $pdate = jtime(shift || time());

	my $pphase;				# illuminated fraction
	my $mage;				# age of moon in days
	my $dist;				# distance in kilometres
	my $angdia;				# angular diameter in degrees
	my $sudist;				# distance to Sun
	my $suangdia;				# sun's angular diameter

	my ($Day, $N, $M, $Ec, $Lambdasun, $ml, $MM, $MN, $Ev, $Ae, $A3, $MmP,
	   $mEc, $A4, $lP, $V, $lPP, $NP, $y, $x, $Lambdamoon, $BetaM,
	   $MoonAge, $MoonPhase,
	   $MoonDist, $MoonDFrac, $MoonAng, $MoonPar,
	   $F, $SunDist, $SunAng,
	   $mpfrac);

	# Calculation of the Sun's position.

	$Day = $pdate - $Epoch;						# date within epoch
	$N = fixangle((360 / 365.2422) * $Day);				# mean anomaly of the Sun
	$M = fixangle($N + $Elonge - $Elongp);				# convert from perigee
									# co-ordinates to epoch 1980.0
	$Ec = kepler($M, $Eccent);					# solve equation of Kepler
	$Ec = sqrt((1 + $Eccent) / (1 - $Eccent)) * tan($Ec / 2);
	$Ec = 2 * todeg(atan($Ec));					# true anomaly
	$Lambdasun = fixangle($Ec + $Elongp);				# Sun's geocentric ecliptic
									# longitude
	# Orbital distance factor.
	$F = ((1 + $Eccent * cos(torad($Ec))) / (1 - $Eccent * $Eccent));
	$SunDist = $Sunsmax / $F;					# distance to Sun in km
	$SunAng = $F * $Sunangsiz;					# Sun's angular size in degrees


	# Calculation of the Moon's position.

	# Moon's mean longitude.
	$ml = fixangle(13.1763966 * $Day + $Mmlong);

	# Moon's mean anomaly.
	$MM = fixangle($ml - 0.1114041 * $Day - $Mmlongp);

	# Moon's ascending node mean longitude.
	$MN = fixangle($Mlnode - 0.0529539 * $Day);

	# Evection.
	$Ev = 1.2739 * sin(torad(2 * ($ml - $Lambdasun) - $MM));

	# Annual equation.
	$Ae = 0.1858 * sin(torad($M));

	# Correction term.
	$A3 = 0.37 * sin(torad($M));

	# Corrected anomaly.
	$MmP = $MM + $Ev - $Ae - $A3;

	# Correction for the equation of the centre.
	$mEc = 6.2886 * sin(torad($MmP));

	# Another correction term.
	$A4 = 0.214 * sin(torad(2 * $MmP));

	# Corrected longitude.
	$lP = $ml + $Ev + $mEc - $Ae + $A4;

	# Variation.
	$V = 0.6583 * sin(torad(2 * ($lP - $Lambdasun)));

	# True longitude.
	$lPP = $lP + $V;

	# Corrected longitude of the node.
	$NP = $MN - 0.16 * sin(torad($M));

	# Y inclination coordinate.
	$y = sin(torad($lPP - $NP)) * cos(torad($Minc));

	# X inclination coordinate.
	$x = cos(torad($lPP - $NP));

	# Ecliptic longitude.
	$Lambdamoon = todeg(atan2($y, $x));
	$Lambdamoon += $NP;

	# Ecliptic latitude.
	$BetaM = todeg(asin(sin(torad($lPP - $NP)) * sin(torad($Minc))));

	# Calculation of the phase of the Moon.

	# Age of the Moon in degrees.
	$MoonAge = $lPP - $Lambdasun;

	# Phase of the Moon.
	$MoonPhase = (1 - cos(torad($MoonAge))) / 2;

	# Calculate distance of moon from the centre of the Earth.

	$MoonDist = ($Msmax * (1 - $Mecc * $Mecc)) /
		(1 + $Mecc * cos(torad($MmP + $mEc)));

	# Calculate Moon's angular diameter.

	$MoonDFrac = $MoonDist / $Msmax;
	$MoonAng = $Mangsiz / $MoonDFrac;

	# Calculate Moon's parallax.

	$MoonPar = $Mparallax / $MoonDFrac;

	$pphase = $MoonPhase;
	$mage = $Synmonth * (fixangle($MoonAge) / 360.0);
	$dist = $MoonDist;
	$angdia = $MoonAng;
	$sudist = $SunDist;
	$suangdia = $SunAng;
	$mpfrac = fixangle($MoonAge) / 360.0;
	return wantarray ? ( $mpfrac, $pphase, $mage, $dist, $angdia, $sudist,$suangdia ) : $mpfrac;
}

1;
__END__

=head1 NAME

Astro::MoonPhase - Information about the phase of the Moon

=head1 SYNOPSIS

use Astro::MoonPhase;

	( $MoonPhase,
	  $MoonIllum,
	  $MoonAge,
	  $MoonDist,
	  $MoonAng,
	  $SunDist,
	  $SunAng ) = phase($seconds_since_1970);

	@phases  = phasehunt($seconds_since_1970);

	($phase, @times) = phaselist($start, $stop);

=head1 DESCRIPTION

MoonPhase calculates information about the phase of the moon
at a given time.

=head1 FUNCTIONS

=head2 phase()

	( $MoonPhase,
	  $MoonIllum,
	  $MoonAge,
	  $MoonDist,
	  $MoonAng,
	  $SunDist,
	  $SunAng )  = phase($seconds_since_1970);

	  $MoonPhase = phase($seconds_since_1970);

The argument is the time for which the phase is requested,
expressed as a time returned by the C<time> function. If C<$seconds_since_1970>
is omitted, it does C<phase(time)>.

Return value in scalar context is $MoonPhase,
the terminator phase angle as a percentage of a full circle (i.e., 0 to 1).

=over 4

=item B<Return values in array context:>

=item $MoonPhase:

the terminator phase angle as a percentage of a full circle (i.e., 0 to 1)

=item $MoonIllum:

the illuminated fraction of the Moon's disc

=item $MoonAge:

the Moon's age in days and fraction

=item $MoonDist:

the distance of the Moon from the centre of the Earth

=item $MoonAng:

the angular diameter subtended by the Moon as seen by
an observer at the centre of the Earth.

=item $SunDist:

the distance from the Sun in km

=item $SunAng:

the angular size of Sun in degrees

=back

Example:

   ( $MoonPhase,
     $MoonIllum,
     $MoonAge,
     $MoonDist,
     $MoonAng,
     $SunDist,
     $SunAng ) = phase();

     print "MoonPhase  = $MoonPhase\n";
     print "MoonIllum  = $MoonIllum\n";
     print "MoonAge    = $MoonAge\n";
     print "MoonDist   = $MoonDist\n";
     print "MoonAng    = $MoonAng\n";
     print "SunDist    = $SunDist\n";
     print "SunAng     = $SunAng\n";

could print something like this:

     MoonPhase  = 0.598939375319023
     MoonIllum  = 0.906458030827876
     MoonAge    = 17.6870323368022
     MoonDist   = 372479.357420033
     MoonAng    = 0.534682403555093
     SunDist    = 152078368.820205
     SunAng     = 0.524434538105092

=head2 phasehunt()

     @phases = phasehunt($seconds_since_1970);

Finds time of phases of the moon which surround the given
date.  Five phases are found, starting and ending with the
new moons which bound the current lunation.

The argument is the time, expressed as a time returned
by the C<time> function. If C<$seconds_since_1970>
is omitted, it does C<phasehunt(time)>.

Example:

    @phases = phasehunt();
    print "New Moon      = ", scalar(localtime($phases[0])), "\n";
    print "First quarter = ", scalar(localtime($phases[1])), "\n";
    print "Full moon     = ", scalar(localtime($phases[2])), "\n";
    print "Last quarter  = ", scalar(localtime($phases[3])), "\n";
    print "New Moon      = ", scalar(localtime($phases[4])), "\n";

could print something like this:

    New Moon      = Wed Jun 24 06:51:47 1998
    First quarter = Wed Jul  1 21:42:19 1998
    Full moon     = Thu Jul  9 19:02:47 1998
    Last quarter  = Thu Jul 16 18:15:18 1998
    New Moon      = Thu Jul 23 16:45:01 1998

=head2 phaselist()

    ($phase, @times) = phaselist($start, $stop);

Finds times of all phases of the moon which occur on or after
C<$start> but before C<$stop>.  Both the arguments and the return
values are expressed as seconds since 1970 (like the C<time> function
returns).

C<$phase> is an integer indicating the phase of the moon at
C<$times[0]>, as shown in this table:

    0  New Moon
    1  First quarter
    2  Full Moon
    3  Last quarter

The remaining values in C<@times> indicate subsequent phases of the
moon (in ascending order by time).  If there are no phases of the moon
between C<$start> and C<$stop>, C<phaselist> returns the empty list.

Example:

    @name = ("New Moon", "First quarter", "Full moon", "Last quarter");
    ($phase, @times) = phaselist($start, $stop);

    while (@times) {
      printf "%-14s= %s\n", $name[$phase], scalar localtime shift @times;
      $phase = ($phase + 1) % 4;
    }

could produce the same output as the C<phasehunt> example above (given
the appropriate start & stop times).

=head1 ABOUT THE ALGORITHMS

The algorithms used in this program to calculate the positions of Sun and
Moon as seen from the Earth are given in the book I<Practical Astronomy
With  Your  Calculator>  by  B<Peter  Duffett-Smith,   Second   Edition,
Cambridge University Press, 1981>.  Ignore the word "Calculator" in the
title;  this  is  an  essential  reference  if  you're  interested  in
developing  software  which  calculates  planetary  positions, orbits,
eclipses, and  the  like.   If  you're  interested  in  pursuing  such
programming, you should also obtain:

I<Astronomical  Formulae for Calculators> by B<Jean Meeus, Third Edition,
Willmann-Bell, 1985>.  A must-have.

I<Planetary  Programs  and  Tables  from  -4000  to  +2800>  by  B<Pierre
Bretagnon  and Jean-Louis Simon, Willmann-Bell, 1986>.  If you want the
utmost  (outside  of  JPL)  accuracy  for  the  planets,  it's   here.

I<Celestial BASIC> by B<Eric Burgess, Revised Edition, Sybex, 1985>.  Very
cookbook oriented, and many of the algorithms are hard to dig  out  of
the turgid BASIC code, but you'll probably want it anyway.

Many of these references can be obtained from Willmann-Bell, P.O.  Box
35025,  Richmond,  VA 23235, USA.  Phone: (804) 320-7016.  In addition
to their own publications, they stock most of the standard  references
for mathematical and positional astronomy.

=head1 LICENCE

This  program is in the public domain: "Do what thou wilt shall be the
whole of the law".

=head1 AUTHORS

The moontool.c Release 2.0:

    A Moon for the Sun
    Designed and implemented by John Walker in December 1987,
    revised and updated in February of 1988.

Initial Perl transcription:

    Raino Pikkarainen, 1998
    raino.pikkarainen@saunalahti.fi

The moontool.c Release 2.4:

    Major enhancements by Ron Hitchens, 1989

Revisions:

    Brett Hamilton  http://simple.be/
    Bug fix, 2003
    Second transcription and bugfixes, 2004

    Christopher J. Madsen  http://www.cjmweb.net/
    Added phaselist function, March 2007