/inc/moonphase.inc.php
PHP | 472 lines | 259 code | 89 blank | 124 comment | 25 complexity | f462c1d0c3d2b1e5e969bcdfb7bd5381 MD5 | raw file
- <?php
- //
- // 2008-10-27
- // File: moonphase.inc.php (http://www.sentry.net/~obsid/moonphase)
- // Calculate information about the phase of the moon at a given time.
- //
- // Based on the Perl module Astro::MoonPhase, version 0.60.
- // http://search.cpan.org/~brett/Astro-MoonPhase-0.60/MoonPhase.pm
- //
- //
- // License:
- //
- // Astro::MoonPhase module is distributed under the public domain,
- // and so is this PHP translation.
- //
- //
- // Credits:
- //
- // 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
- //
- // Translated to PHP by Stephen A. Zarkos <obsid@sentry.net>, 2007
- // Fixed broken phasehunt function, 2008-10-27
- // Added phaselist function, 2008-10-27
- //
- //
- // Documentation:
- //
- // http://search.cpan.org/~brett/Astro-MoonPhase-0.60/MoonPhase.pm
- // http://www.obsid.org/2008/05/calculate-moon-phase-data-with-php.html
- //
- // Error definitions.
- define( 'ERR_UNDEF', -1 );
- // Astronomical constants.
- define( 'EPOCH', 2444238.5 ); // 1980 January 0.0
- // Constants defining the Sun's apparent orbit.
- define( 'ELONGE', 278.833540 ); // ecliptic longitude of the Sun at epoch 1980.0
- define( 'ELONGP', 282.596403 ); // ecliptic longitude of the Sun at perigee
- define( 'ECCENT', 0.016718 ); // eccentricity of Earth's orbit
- define( 'SUNSMAX', 1.495985e8 ); // semi-major axis of Earth's orbit, km
- define( 'SUNANGSIZ', 0.533128 ); // sun's angular size, degrees, at semi-major axis distance
- // Elements of the Moon's orbit, epoch 1980.0.
- define( 'MMLONG', 64.975464 ); // moon's mean longitude at the epoch
- define( 'MMLONGP', 349.383063 ); // mean longitude of the perigee at the epoch
- define( 'MLNODE', 151.950429 ); // mean longitude of the node at the epoch
- define( 'MINC', 5.145396 ); // inclination of the Moon's orbit
- define( 'MECC', 0.054900 ); // eccentricity of the Moon's orbit
- define( 'MANGSIZ', 0.5181 ); // moon's angular size at distance a from Earth
- define( 'MSMAX', 384401.0 ); // semi-major axis of Moon's orbit in km
- define( 'MPARALLAX', 0.9507 ); // parallax at distance a from Earth
- define( 'SYNMONTH', 29.53058868 ); // synodic month (new Moon to new Moon)
- // Handy mathematical functions.
- function sgn ( $arg ) { return (($arg < 0) ? -1 : ($arg > 0 ? 1 : 0)); } // extract sign
- function fixangle ( $arg ) { return ($arg - 360.0 * (floor($arg / 360.0))); } // fix angle
- function torad ( $arg ) { return ($arg * (pi() / 180.0)); } // deg->rad
- function todeg ( $arg ) { return ($arg * (180.0 / pi())); } // rad->deg
- function dsin ( $arg ) { return (sin(torad($arg))); } // sin from deg
- function dcos ( $arg ) { return (cos(torad($arg))); } // cos from deg
- // jtime - convert internal date and time to astronomical Julian
- // time (i.e. Julian date plus day fraction)
- function jtime ( $timestamp ) {
- $julian = ( $timestamp / 86400 ) + 2440587.5; // (seconds / (seconds per day)) + julian date of epoch
- return $julian;
- }
- // jdaytosecs - convert Julian date to a UNIX epoch
- function jdaytosecs ( $jday=0 ) {
- $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
- function jyear ( $td, &$yy, &$mm, &$dd ) {
- $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.
- function meanphase ( $sdate, $k ) {
-
- // 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.
- function truephase ( $k, $phase ) {
- $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;
- }
- elseif ( (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 ) {
- print "truephase() called with invalid phase selector ($phase).\n";
- exit( ERR_UNDEF );
- }
- 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
- function phasehunt ( $time=-1 ) {
-
-
- $k1k2 = time2k1k2($time);
- $k1 = $k1k2[0];
- $k2 = $k1k2[0];
- return array ( 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) )
- );
- }
- function time2k1k2($time = -1) {
-
- if ( empty($time) || $time == -1 ) {
- $time = time();
- }
- $sdate = jtime( $time );
- $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)) {
- break;
- }
- $nt1 = $nt2;
- $k1 = $k2;
- }
- return array($k1,$k2);
- }
- function phasemoon($time = -1,$phase = 0.5,$second = false) {
- $k1k2 = time2k1k2($time);
- $k1 = $k1k2[0];
- if ($second) {
- $k1 = $k1k2[1];
- }
- return jdaytosecs( truephase($k1, $phase) );
- }
- // phaselist() - Find time of phases of the moon between two dates.
- // Times (in & out) are seconds_since_1970
- function phaselist ( $sdate, $edate ) {
- if ( empty($sdate) || empty($edate) ) {
- return array();
- }
-
- $sdate = jtime( $sdate );
- $edate = jtime( $edate );
-
- $phases = array();
- $d = $k = $yy = $mm = 0;
-
- jyear( $sdate, $yy, $mm, $d );
- $k = floor(($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685) - 2;
-
- while (1) {
- ++$k;
- foreach ( array(0.0, 0.25, 0.5, 0.75) as $phase ) {
- $d = truephase( $k, $phase );
- if ( $d >= $edate ) {
- return $phases;
- }
- if ( $d >= $sdate ) {
- if ( empty($phases) ) {
- array_push( $phases, floor(4 * $phase) );
- }
- array_push( $phases, jdaytosecs($d) );
- }
- }
- } // End while(1)
- }
- // kepler() - solve the equation of Kepler
- function kepler ( $m, $ecc ) {
- $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.
- function phase ( $time=0 ) {
- if ( empty($time) || $time == 0 ) {
- $time = time();
- }
- $pdate = jtime( $time );
-
-
- $pphase = 0; // illuminated fraction
- $mage = 0; // age of moon in days
- $dist = 0; // distance in kilometres
- $angdia = 0; // angular diameter in degrees
- $sudist = 0; // distance to Sun
- $suangdia = 0; // 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 array ( $mpfrac, $pphase, $mage, $dist, $angdia, $sudist, $suangdia );
- }