package Astro::MoonPhaseModified;
use strict;
use vars qw($VERSION @ISA @EXPORT @EXPORT_OK);
require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw(phase phasehunt);
$VERSION = '0.10';
use Time::Local qw(timegm);
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);
}
# jdate - convert internal GMT date and time to Julian day and fraction
### jdate is no longer needed by jtime.
sub jdate {
use integer;
my ($sec, $min, $hour, $mday, $mon, $year, $wday, $yday, $isdst) = @_;
my ($c, $m, $y);
$y = $year + 1900;
$m = $mon + 1;
if ($m > 2) {
$m = $m - 3;
}
else {
$m = $m + 9;
$y--;
}
$c = $y / 100; # compute century
$y -= 100 * $c;
return ($mday + ($c * 146097) / 4 + ($y * 1461) / 4 + ($m * 153 + 2) / 5 + 1721119);
}
# jtime - convert internal date and time to astronomical Julian
# time (i.e. Julian date plus day fraction)
####### new jtime subroutine ###############
sub jtime {
my $t = shift;
my ($julian);
$julian = ($t / 86400) + 2440587.5;
return ($julian);
}
####### previous jtime subroutine ###############p
# sub jtime {
#
# my (@dt, $julian);
# @dt = localtime($t); ##### <---- should be: @dt = gmtime($t);
# return (( jdate(@dt) - 0.5 ) + ( $dt[0] + 60 * ( $dt[1] + 60 * $dt[2] ) ) / 86400.0 );
#
# }
# jyear - convert Julian date to year, month, day, which are
# returned via integer pointers to integers
sub jyear {
my $td = shift;
my ($yy, $mm, $dd) = @_;
my ($j, $d, $y, $m);
$td += 0.5; # astronomical to civil
$j = floor($td);
$j = $j - 1721119.0;
$y = floor(((4 * $j) - 1) / 146097.0);
$j = ($j * 4.0) - (1.0 + (146097.0 * $y));
$d = floor($j / 4.0);
$j = floor(((4.0 * $d) + 3.0) / 1461.0);
$d = ((4.0 * $d) + 3.0) - (1461.0 * $j);
$d = floor(($d + 4.0) / 4.0);
$m = floor(((5.0 * $d) - 3) / 153.0);
$d = (5.0 * $d) - (3.0 + (153.0 * $m));
$d = floor(($d + 5.0) / 5.0);
$y = (100.0 * $y) + $j;
if ($m < 10.0) {
$m = $m + 3;
}
else {
$m = $m - 9;
$y = $y + 1;
}
$$yy = int $y;
$$mm = int $m;
$$dd = int $d;
}
# jhms - convert Julian time to hour, minutes, and seconds
sub jhms {
my $j = shift;
my ($ij, $h, $m, $s);
$j += 0.5; # astronomical to civil
$ij = int (($j - floor($j)) * 86400.0);
$h = int ($ij / 3600);
$m = int (($ij / 60) % 60);
$s = int ($ij % 60);
return (($h, $m, $s));
}
# jdaytosecs - convert Julian time to a time returned by time() function
sub jdaytosecs {
my $jday = shift;
my @hms = jhms($jday);
my ($y, $m, $d);
jyear($jday, \$y, \$m, \$d);
return ( timegm($hms[2], $hms[1], $hms[0], $d, --$m, $y) );
}
# meanphase - calculates mean phase of the Moon for a given base date
# and desired phase:
# 0.0 New Moon
# 0.25 First quarter
# 0.5 Full moon
# 0.75 Last quarter
# Beware!!! This routine returns meaningless
# results for any other phase arguments. Don't
# attempt to generalise it without understanding
# that the motion of the moon is far more complicated
# that this calculation reveals.
sub meanphase {
my ($sdate, $phase, $usek) = @_;
my ($yy, $mm, $dd);
my ($k, $t, $t2, $t3, $nt1);
jyear($sdate, \$yy, \$mm, \$dd);
$k = ($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685;
# 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
$$usek = $k = floor($k) + $phase;
$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);
$adate = $sdate - 45;
$nt1 = meanphase($adate, 0.0, \$k1);
while (1) {
$adate += $Synmonth;
$nt2 = meanphase($adate, 0.0, \$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))
);
}
# 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
MoonPhase - Information about the phase of the moon
=head1 SYNOPSIS
use MoonPhase;
( $MoonPhase,
$MoonIllum,
$MoonAge,
$MoonDist,
$MoonAng,
$SunDist,
$SunAng ) = phase($seconds_since_1970);
@phases = phasehunt($seconds_since_1970);
=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