moonphase.inc.php | searchcode

/inc/moonphase.inc.php

https://github.com/chregu/isitfullmoon
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 );
}