Skip to content

Instantly share code, notes, and snippets.

@aaryadev
Created April 3, 2015 16:10
Show Gist options
  • Save aaryadev/862c8e7706d497a22716 to your computer and use it in GitHub Desktop.
Save aaryadev/862c8e7706d497a22716 to your computer and use it in GitHub Desktop.
MoonPhase
<?php
class MoonPhase {
private $timestamp;
private $phase;
private $illum;
private $age;
private $dist;
private $angdia;
private $sundist;
private $sunangdia;
private $synmonth;
private $quarters = null;
function __construct( $pdate = null ) {
if( is_null( $pdate ) )
$pdate = time();
/* 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; // 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)
$this->synmonth = $synmonth;
$lunatbase = 2423436.0; // Base date for E. W. Brown's numbered series of lunations (1923 January 16)
/* Properties of the Earth */
// $earthrad = 6378.16; // Radius of Earth in kilometres
// $PI = 3.14159265358979323846; // Assume not near black hole
$this->timestamp = $pdate;
// pdate is coming in as a UNIX timstamp, so convert it to Julian
$pdate = $pdate / 86400 + 2440587.5;
/* Calculation of the Sun's position */
$Day = $pdate - $epoch; // Date within epoch
$N = $this->fixangle((360 / 365.2422) * $Day); // Mean anomaly of the Sun
$M = $this->fixangle($N + $elonge - $elongp); // Convert from perigee co-ordinates to epoch 1980.0
$Ec = $this->kepler($M, $eccent); // Solve equation of Kepler
$Ec = sqrt((1 + $eccent) / (1 - $eccent)) * tan($Ec / 2);
$Ec = 2 * rad2deg(atan($Ec)); // True anomaly
$Lambdasun = $this->fixangle($Ec + $elongp); // Sun's geocentric ecliptic longitude
$F = ((1 + $eccent * cos(deg2rad($Ec))) / (1 - $eccent * $eccent)); // Orbital distance factor
$SunDist = $sunsmax / $F; // Distance to Sun in km
$SunAng = $F * $sunangsiz; // Sun's angular size in degrees
/* Calculation of the Moon's position */
$ml = $this->fixangle(13.1763966 * $Day + $mmlong); // Moon's mean longitude
$MM = $this->fixangle($ml - 0.1114041 * $Day - $mmlongp); // Moon's mean anomaly
$MN = $this->fixangle($mlnode - 0.0529539 * $Day); // Moon's ascending node mean longitude
$Ev = 1.2739 * sin(deg2rad(2 * ($ml - $Lambdasun) - $MM)); // Evection
$Ae = 0.1858 * sin(deg2rad($M)); // Annual equation
$A3 = 0.37 * sin(deg2rad($M)); // Correction term
$MmP = $MM + $Ev - $Ae - $A3; // Corrected anomaly
$mEc = 6.2886 * sin(deg2rad($MmP)); // Correction for the equation of the centre
$A4 = 0.214 * sin(deg2rad(2 * $MmP)); // Another correction term
$lP = $ml + $Ev + $mEc - $Ae + $A4; // Corrected longitude
$V = 0.6583 * sin(deg2rad(2 * ($lP - $Lambdasun))); // Variation
$lPP = $lP + $V; // True longitude
$NP = $MN - 0.16 * sin(deg2rad($M)); // Corrected longitude of the node
$y = sin(deg2rad($lPP - $NP)) * cos(deg2rad($minc)); // Y inclination coordinate
$x = cos(deg2rad($lPP - $NP)); // X inclination coordinate
$Lambdamoon = rad2deg(atan2($y, $x)) + $NP; // Ecliptic longitude
$BetaM = rad2deg(asin(sin(deg2rad($lPP - $NP)) * sin(deg2rad($minc)))); // Ecliptic latitude
/* Calculation of the phase of the Moon */
$MoonAge = $lPP - $Lambdasun; // Age of the Moon in degrees
$MoonPhase = (1 - cos(deg2rad($MoonAge))) / 2; // Phase of the Moon
// Distance of moon from the centre of the Earth
$MoonDist = ($msmax * (1 - $mecc * $mecc)) / (1 + $mecc * cos(deg2rad($MmP + $mEc)));
$MoonDFrac = $MoonDist / $msmax;
$MoonAng = $mangsiz / $MoonDFrac; // Moon's angular diameter
// $MoonPar = $mparallax / $MoonDFrac; // Moon's parallax
// store results
$this->phase = $this->fixangle($MoonAge) / 360; // Phase (0 to 1)
$this->illum = $MoonPhase; // Illuminated fraction (0 to 1)
$this->age = $synmonth * $this->phase; // Age of moon (days)
$this->dist = $MoonDist; // Distance (kilometres)
$this->angdia = $MoonAng; // Angular diameter (degrees)
$this->sundist = $SunDist; // Distance to Sun (kilometres)
$this->sunangdia = $SunAng; // Sun's angular diameter (degrees)
}
private function fixangle($a) {
return ( $a - 360 * floor($a / 360) );
}
// KEPLER -- Solve the equation of Kepler.
private function kepler($m, $ecc) {
$epsilon = 0.000001; // 1E-6
$e = $m = deg2rad($m);
do {
$delta = $e - $ecc * sin($e) - $m;
$e -= $delta / ( 1 - $ecc * cos($e) );
}
while ( abs($delta) > $epsilon );
return $e;
}
/* 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.
*/
private function meanphase($sdate, $k){
// Time in Julian centuries from 1900 January 0.5
$t = ( $sdate - 2415020.0 ) / 36525;
$t2 = $t * $t;
$t3 = $t2 * $t;
$nt1 = 2415020.75933 + $this->synmonth * $k
+ 0.0001178 * $t2
- 0.000000155 * $t3
+ 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) );
return $nt1;
}
/* 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.
*/
private function truephase($k, $phase){
$apcor = false;
$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
$pt = 2415020.75933 // Mean time of phase
+ $this->synmonth * $k
+ 0.0001178 * $t2
- 0.000000155 * $t3
+ 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) );
$m = 359.2242 + 29.10535608 * $k - 0.0000333 * $t2 - 0.00000347 * $t3; // Sun's mean anomaly
$mprime = 306.0253 + 385.81691806 * $k + 0.0107306 * $t2 + 0.00001236 * $t3; // Moon's mean anomaly
$f = 21.2964 + 390.67050646 * $k - 0.0016528 * $t2 - 0.00000239 * $t3; // Moon's argument of latitude
if ( $phase < 0.01 || abs( $phase - 0.5 ) < 0.01 ) {
// Corrections for New and Full Moon
$pt += (0.1734 - 0.000393 * $t) * sin( deg2rad( $m ) )
+ 0.0021 * sin( deg2rad( 2 * $m ) )
- 0.4068 * sin( deg2rad( $mprime ) )
+ 0.0161 * sin( deg2rad( 2 * $mprime) )
- 0.0004 * sin( deg2rad( 3 * $mprime ) )
+ 0.0104 * sin( deg2rad( 2 * $f ) )
- 0.0051 * sin( deg2rad( $m + $mprime ) )
- 0.0074 * sin( deg2rad( $m - $mprime ) )
+ 0.0004 * sin( deg2rad( 2 * $f + $m ) )
- 0.0004 * sin( deg2rad( 2 * $f - $m ) )
- 0.0006 * sin( deg2rad( 2 * $f + $mprime ) )
+ 0.0010 * sin( deg2rad( 2 * $f - $mprime ) )
+ 0.0005 * sin( deg2rad( $m + 2 * $mprime ) );
$apcor = true;
} else if ( abs( $phase - 0.25 ) < 0.01 || abs( $phase - 0.75 ) < 0.01 ) {
$pt += (0.1721 - 0.0004 * $t) * sin( deg2rad( $m ) )
+ 0.0021 * sin( deg2rad( 2 * $m ) )
- 0.6280 * sin( deg2rad( $mprime ) )
+ 0.0089 * sin( deg2rad( 2 * $mprime) )
- 0.0004 * sin( deg2rad( 3 * $mprime ) )
+ 0.0079 * sin( deg2rad( 2 * $f ) )
- 0.0119 * sin( deg2rad( $m + $mprime ) )
- 0.0047 * sin( deg2rad ( $m - $mprime ) )
+ 0.0003 * sin( deg2rad( 2 * $f + $m ) )
- 0.0004 * sin( deg2rad( 2 * $f - $m ) )
- 0.0006 * sin( deg2rad( 2 * $f + $mprime ) )
+ 0.0021 * sin( deg2rad( 2 * $f - $mprime ) )
+ 0.0003 * sin( deg2rad( $m + 2 * $mprime ) )
+ 0.0004 * sin( deg2rad( $m - 2 * $mprime ) )
- 0.0003 * sin( deg2rad( 2 * $m + $mprime ) );
if ( $phase < 0.5 ) // First quarter correction
$pt += 0.0028 - 0.0004 * cos( deg2rad( $m ) ) + 0.0003 * cos( deg2rad( $mprime ) );
else // Last quarter correction
$pt += -0.0028 + 0.0004 * cos( deg2rad( $m ) ) - 0.0003 * cos( deg2rad( $mprime ) );
$apcor = true;
}
if (!$apcor) // function was called with an invalid phase selector
return false;
return $pt;
}
/* 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.
*/
private function phasehunt() {
$sdate = $this->utctojulian( $this->timestamp );
$adate = $sdate - 45;
$ats = $this->timestamp - 86400 * 45;
$yy = (int) gmdate( 'Y', $ats );
$mm = (int) gmdate( 'n', $ats );
$k1 = floor( ( $yy + ( ( $mm - 1 ) * ( 1 / 12 ) ) - 1900 ) * 12.3685 );
$adate = $nt1 = $this->meanphase( $adate, $k1 );
while (true) {
$adate += $this->synmonth;
$k2 = $k1 + 1;
$nt2 = $this->meanphase( $adate, $k2 );
// if nt2 is close to sdate, then mean phase isn't good enough, we have to be more accurate
if( abs( $nt2 - $sdate ) < 0.5 )
$nt2 = $this->truephase( $k2, 0.0 );
if ( $nt1 <= $sdate && $nt2 > $sdate )
break;
$nt1 = $nt2;
$k1 = $k2;
}
// results in Julian dates
$data = array(
$this->truephase( $k1, 0.0 ),
$this->truephase( $k1, 0.25 ),
$this->truephase( $k1, 0.5 ),
$this->truephase( $k1, 0.75 ),
$this->truephase( $k2, 0.0 )
);
$this->quarters = array();
foreach( $data as $v )
$this->quarters[] = ( $v - 2440587.5 ) * 86400; // convert to UNIX time
}
/* Convert UNIX timestamp to astronomical Julian time (i.e. Julian date plus day fraction). */
private function utctojulian( $ts ) {
return $ts / 86400 + 2440587.5;
}
private function get_phase( $n ) {
if( is_null( $this->quarters ) )
$this->phasehunt();
return $this->quarters[$n];
}
/* Public functions for accessing results */
function phase(){
return $this->phase;
}
function illumination(){
return $this->illum;
}
function age(){
return $this->age;
}
function distance(){
return $this->dist;
}
function diameter(){
return $this->angdia;
}
function sundistance(){
return $this->sundist;
}
function sundiameter(){
return $this->sunangdia;
}
function new_moon(){
return $this->get_phase( 0 );
}
function first_quarter(){
return $this->get_phase( 1 );
}
function full_moon(){
return $this->get_phase( 2 );
}
function last_quarter(){
return $this->get_phase( 3 );
}
function next_new_moon(){
return $this->get_phase( 4 );
}
}
@aaryadev
Copy link
Author

aaryadev commented Apr 3, 2015

MoonPhase calculation

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment