| 1 |
<?php |
| 2 |
/** |
| 3 |
* Moon phase calculation class |
| 4 |
* Adapted for PHP from Moontool for Windows (http://www.fourmilab.ch/moontoolw/) |
| 5 |
* by Samir Shah (http://rayofsolaris.net) |
| 6 |
* License: MIT |
| 7 |
**/ |
| 8 |
|
| 9 |
class Solaris_MoonPhase { |
| 10 |
private $timestamp; |
| 11 |
private $phase; |
| 12 |
private $illum; |
| 13 |
private $age; |
| 14 |
private $dist; |
| 15 |
private $angdia; |
| 16 |
private $sundist; |
| 17 |
private $sunangdia; |
| 18 |
|
| 19 |
private $synmonth; |
| 20 |
|
| 21 |
private $quarters = null; |
| 22 |
|
| 23 |
function __construct( $pdate = null ) { |
| 24 |
if( is_null( $pdate ) ) |
| 25 |
$pdate = time(); |
| 26 |
|
| 27 |
/* Astronomical constants */ |
| 28 |
$epoch = 2444238.5; // 1980 January 0.0 |
| 29 |
|
| 30 |
/* Constants defining the Sun's apparent orbit */ |
| 31 |
$elonge = 278.833540; // Ecliptic longitude of the Sun at epoch 1980.0 |
| 32 |
$elongp = 282.596403; // Ecliptic longitude of the Sun at perigee |
| 33 |
$eccent = 0.016718; // Eccentricity of Earth's orbit |
| 34 |
$sunsmax = 1.495985e8; // Semi-major axis of Earth's orbit, km |
| 35 |
$sunangsiz = 0.533128; // Sun's angular size, degrees, at semi-major axis distance |
| 36 |
|
| 37 |
/* Elements of the Moon's orbit, epoch 1980.0 */ |
| 38 |
$mmlong = 64.975464; // Moon's mean longitude at the epoch |
| 39 |
$mmlongp = 349.383063; // Mean longitude of the perigee at the epoch |
| 40 |
$mlnode = 151.950429; // Mean longitude of the node at the epoch |
| 41 |
$minc = 5.145396; // Inclination of the Moon's orbit |
| 42 |
$mecc = 0.054900; // Eccentricity of the Moon's orbit |
| 43 |
$mangsiz = 0.5181; // Moon's angular size at distance a from Earth |
| 44 |
$msmax = 384401; // Semi-major axis of Moon's orbit in km |
| 45 |
$mparallax = 0.9507; // Parallax at distance a from Earth |
| 46 |
$synmonth = 29.53058868; // Synodic month (new Moon to new Moon) |
| 47 |
$this->synmonth = $synmonth; |
| 48 |
$lunatbase = 2423436.0; // Base date for E. W. Brown's numbered series of lunations (1923 January 16) |
| 49 |
|
| 50 |
/* Properties of the Earth */ |
| 51 |
// $earthrad = 6378.16; // Radius of Earth in kilometres |
| 52 |
// $PI = 3.14159265358979323846; // Assume not near black hole |
| 53 |
|
| 54 |
$this->timestamp = $pdate; |
| 55 |
|
| 56 |
// pdate is coming in as a UNIX timstamp, so convert it to Julian |
| 57 |
$pdate = $pdate / 86400 + 2440587.5; |
| 58 |
|
| 59 |
/* Calculation of the Sun's position */ |
| 60 |
|
| 61 |
$Day = $pdate - $epoch; // Date within epoch |
| 62 |
$N = $this->fixangle((360 / 365.2422) * $Day); // Mean anomaly of the Sun |
| 63 |
$M = $this->fixangle($N + $elonge - $elongp); // Convert from perigee co-ordinates to epoch 1980.0 |
| 64 |
$Ec = $this->kepler($M, $eccent); // Solve equation of Kepler |
| 65 |
$Ec = sqrt((1 + $eccent) / (1 - $eccent)) * tan($Ec / 2); |
| 66 |
$Ec = 2 * rad2deg(atan($Ec)); // True anomaly |
| 67 |
$Lambdasun = $this->fixangle($Ec + $elongp); // Sun's geocentric ecliptic longitude |
| 68 |
|
| 69 |
$F = ((1 + $eccent * cos(deg2rad($Ec))) / (1 - $eccent * $eccent)); // Orbital distance factor |
| 70 |
$SunDist = $sunsmax / $F; // Distance to Sun in km |
| 71 |
$SunAng = $F * $sunangsiz; // Sun's angular size in degrees |
| 72 |
|
| 73 |
/* Calculation of the Moon's position */ |
| 74 |
$ml = $this->fixangle(13.1763966 * $Day + $mmlong); // Moon's mean longitude |
| 75 |
$MM = $this->fixangle($ml - 0.1114041 * $Day - $mmlongp); // Moon's mean anomaly |
| 76 |
$MN = $this->fixangle($mlnode - 0.0529539 * $Day); // Moon's ascending node mean longitude |
| 77 |
$Ev = 1.2739 * sin(deg2rad(2 * ($ml - $Lambdasun) - $MM)); // Evection |
| 78 |
$Ae = 0.1858 * sin(deg2rad($M)); // Annual equation |
| 79 |
$A3 = 0.37 * sin(deg2rad($M)); // Correction term |
| 80 |
$MmP = $MM + $Ev - $Ae - $A3; // Corrected anomaly |
| 81 |
$mEc = 6.2886 * sin(deg2rad($MmP)); // Correction for the equation of the centre |
| 82 |
$A4 = 0.214 * sin(deg2rad(2 * $MmP)); // Another correction term |
| 83 |
$lP = $ml + $Ev + $mEc - $Ae + $A4; // Corrected longitude |
| 84 |
$V = 0.6583 * sin(deg2rad(2 * ($lP - $Lambdasun))); // Variation |
| 85 |
$lPP = $lP + $V; // True longitude |
| 86 |
$NP = $MN - 0.16 * sin(deg2rad($M)); // Corrected longitude of the node |
| 87 |
$y = sin(deg2rad($lPP - $NP)) * cos(deg2rad($minc)); // Y inclination coordinate |
| 88 |
$x = cos(deg2rad($lPP - $NP)); // X inclination coordinate |
| 89 |
|
| 90 |
$Lambdamoon = rad2deg(atan2($y, $x)) + $NP; // Ecliptic longitude |
| 91 |
$BetaM = rad2deg(asin(sin(deg2rad($lPP - $NP)) * sin(deg2rad($minc)))); // Ecliptic latitude |
| 92 |
|
| 93 |
/* Calculation of the phase of the Moon */ |
| 94 |
$MoonAge = $lPP - $Lambdasun; // Age of the Moon in degrees |
| 95 |
$MoonPhase = (1 - cos(deg2rad($MoonAge))) / 2; // Phase of the Moon |
| 96 |
|
| 97 |
// Distance of moon from the centre of the Earth |
| 98 |
$MoonDist = ($msmax * (1 - $mecc * $mecc)) / (1 + $mecc * cos(deg2rad($MmP + $mEc))); |
| 99 |
|
| 100 |
$MoonDFrac = $MoonDist / $msmax; |
| 101 |
$MoonAng = $mangsiz / $MoonDFrac; // Moon's angular diameter |
| 102 |
// $MoonPar = $mparallax / $MoonDFrac; // Moon's parallax |
| 103 |
|
| 104 |
// store results |
| 105 |
$this->phase = $this->fixangle($MoonAge) / 360; // Phase (0 to 1) |
| 106 |
$this->illum = $MoonPhase; // Illuminated fraction (0 to 1) |
| 107 |
$this->age = $synmonth * $this->phase; // Age of moon (days) |
| 108 |
$this->dist = $MoonDist; // Distance (kilometres) |
| 109 |
$this->angdia = $MoonAng; // Angular diameter (degrees) |
| 110 |
$this->sundist = $SunDist; // Distance to Sun (kilometres) |
| 111 |
$this->sunangdia = $SunAng; // Sun's angular diameter (degrees) |
| 112 |
} |
| 113 |
|
| 114 |
private function fixangle($a) { |
| 115 |
return ( $a - 360 * floor($a / 360) ); |
| 116 |
} |
| 117 |
|
| 118 |
// KEPLER -- Solve the equation of Kepler. |
| 119 |
private function kepler($m, $ecc) { |
| 120 |
$epsilon = 0.000001; // 1E-6 |
| 121 |
$e = $m = deg2rad($m); |
| 122 |
do { |
| 123 |
$delta = $e - $ecc * sin($e) - $m; |
| 124 |
$e -= $delta / ( 1 - $ecc * cos($e) ); |
| 125 |
} |
| 126 |
while ( abs($delta) > $epsilon ); |
| 127 |
return $e; |
| 128 |
} |
| 129 |
|
| 130 |
/* Calculates time of the mean new Moon for a given |
| 131 |
base date. This argument K to this function is the |
| 132 |
precomputed synodic month index, given by: |
| 133 |
K = (year - 1900) * 12.3685 |
| 134 |
where year is expressed as a year and fractional year. |
| 135 |
*/ |
| 136 |
private function meanphase($sdate, $k){ |
| 137 |
// Time in Julian centuries from 1900 January 0.5 |
| 138 |
$t = ( $sdate - 2415020.0 ) / 36525; |
| 139 |
$t2 = $t * $t; |
| 140 |
$t3 = $t2 * $t; |
| 141 |
|
| 142 |
$nt1 = 2415020.75933 + $this->synmonth * $k |
| 143 |
+ 0.0001178 * $t2 |
| 144 |
- 0.000000155 * $t3 |
| 145 |
+ 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) ); |
| 146 |
|
| 147 |
return $nt1; |
| 148 |
} |
| 149 |
|
| 150 |
/* Given a K value used to determine the mean phase of |
| 151 |
the new moon, and a phase selector (0.0, 0.25, 0.5, |
| 152 |
0.75), obtain the true, corrected phase time. |
| 153 |
*/ |
| 154 |
private function truephase($k, $phase){ |
| 155 |
$apcor = false; |
| 156 |
|
| 157 |
$k += $phase; // Add phase to new moon time |
| 158 |
$t = $k / 1236.85; // Time in Julian centuries from 1900 January 0.5 |
| 159 |
$t2 = $t * $t; // Square for frequent use |
| 160 |
$t3 = $t2 * $t; // Cube for frequent use |
| 161 |
$pt = 2415020.75933 // Mean time of phase |
| 162 |
+ $this->synmonth * $k |
| 163 |
+ 0.0001178 * $t2 |
| 164 |
- 0.000000155 * $t3 |
| 165 |
+ 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) ); |
| 166 |
|
| 167 |
$m = 359.2242 + 29.10535608 * $k - 0.0000333 * $t2 - 0.00000347 * $t3; // Sun's mean anomaly |
| 168 |
$mprime = 306.0253 + 385.81691806 * $k + 0.0107306 * $t2 + 0.00001236 * $t3; // Moon's mean anomaly |
| 169 |
$f = 21.2964 + 390.67050646 * $k - 0.0016528 * $t2 - 0.00000239 * $t3; // Moon's argument of latitude |
| 170 |
if ( $phase < 0.01 || abs( $phase - 0.5 ) < 0.01 ) { |
| 171 |
// Corrections for New and Full Moon |
| 172 |
$pt += (0.1734 - 0.000393 * $t) * sin( deg2rad( $m ) ) |
| 173 |
+ 0.0021 * sin( deg2rad( 2 * $m ) ) |
| 174 |
- 0.4068 * sin( deg2rad( $mprime ) ) |
| 175 |
+ 0.0161 * sin( deg2rad( 2 * $mprime) ) |
| 176 |
- 0.0004 * sin( deg2rad( 3 * $mprime ) ) |
| 177 |
+ 0.0104 * sin( deg2rad( 2 * $f ) ) |
| 178 |
- 0.0051 * sin( deg2rad( $m + $mprime ) ) |
| 179 |
- 0.0074 * sin( deg2rad( $m - $mprime ) ) |
| 180 |
+ 0.0004 * sin( deg2rad( 2 * $f + $m ) ) |
| 181 |
- 0.0004 * sin( deg2rad( 2 * $f - $m ) ) |
| 182 |
- 0.0006 * sin( deg2rad( 2 * $f + $mprime ) ) |
| 183 |
+ 0.0010 * sin( deg2rad( 2 * $f - $mprime ) ) |
| 184 |
+ 0.0005 * sin( deg2rad( $m + 2 * $mprime ) ); |
| 185 |
$apcor = true; |
| 186 |
} else if ( abs( $phase - 0.25 ) < 0.01 || abs( $phase - 0.75 ) < 0.01 ) { |
| 187 |
$pt += (0.1721 - 0.0004 * $t) * sin( deg2rad( $m ) ) |
| 188 |
+ 0.0021 * sin( deg2rad( 2 * $m ) ) |
| 189 |
- 0.6280 * sin( deg2rad( $mprime ) ) |
| 190 |
+ 0.0089 * sin( deg2rad( 2 * $mprime) ) |
| 191 |
- 0.0004 * sin( deg2rad( 3 * $mprime ) ) |
| 192 |
+ 0.0079 * sin( deg2rad( 2 * $f ) ) |
| 193 |
- 0.0119 * sin( deg2rad( $m + $mprime ) ) |
| 194 |
- 0.0047 * sin( deg2rad ( $m - $mprime ) ) |
| 195 |
+ 0.0003 * sin( deg2rad( 2 * $f + $m ) ) |
| 196 |
- 0.0004 * sin( deg2rad( 2 * $f - $m ) ) |
| 197 |
- 0.0006 * sin( deg2rad( 2 * $f + $mprime ) ) |
| 198 |
+ 0.0021 * sin( deg2rad( 2 * $f - $mprime ) ) |
| 199 |
+ 0.0003 * sin( deg2rad( $m + 2 * $mprime ) ) |
| 200 |
+ 0.0004 * sin( deg2rad( $m - 2 * $mprime ) ) |
| 201 |
- 0.0003 * sin( deg2rad( 2 * $m + $mprime ) ); |
| 202 |
if ( $phase < 0.5 ) // First quarter correction |
| 203 |
$pt += 0.0028 - 0.0004 * cos( deg2rad( $m ) ) + 0.0003 * cos( deg2rad( $mprime ) ); |
| 204 |
else // Last quarter correction |
| 205 |
$pt += -0.0028 + 0.0004 * cos( deg2rad( $m ) ) - 0.0003 * cos( deg2rad( $mprime ) ); |
| 206 |
$apcor = true; |
| 207 |
} |
| 208 |
if (!$apcor) // function was called with an invalid phase selector |
| 209 |
return false; |
| 210 |
|
| 211 |
return $pt; |
| 212 |
} |
| 213 |
|
| 214 |
/* Find time of phases of the moon which surround the current date. |
| 215 |
Five phases are found, starting and |
| 216 |
ending with the new moons which bound the current lunation. |
| 217 |
*/ |
| 218 |
private function phasehunt() { |
| 219 |
$sdate = $this->utctojulian( $this->timestamp ); |
| 220 |
$adate = $sdate - 45; |
| 221 |
$ats = $this->timestamp - 86400 * 45; |
| 222 |
$yy = (int) gmdate( 'Y', $ats ); |
| 223 |
$mm = (int) gmdate( 'n', $ats ); |
| 224 |
|
| 225 |
$k1 = floor( ( $yy + ( ( $mm - 1 ) * ( 1 / 12 ) ) - 1900 ) * 12.3685 ); |
| 226 |
$adate = $nt1 = $this->meanphase( $adate, $k1 ); |
| 227 |
|
| 228 |
while (true) { |
| 229 |
$adate += $this->synmonth; |
| 230 |
$k2 = $k1 + 1; |
| 231 |
$nt2 = $this->meanphase( $adate, $k2 ); |
| 232 |
// if nt2 is close to sdate, then mean phase isn't good enough, we have to be more accurate |
| 233 |
if( abs( $nt2 - $sdate ) < 0.5 ) |
| 234 |
$nt2 = $this->truephase( $k2, 0.0 ); |
| 235 |
if ( $nt1 <= $sdate && $nt2 > $sdate ) |
| 236 |
break; |
| 237 |
$nt1 = $nt2; |
| 238 |
$k1 = $k2; |
| 239 |
} |
| 240 |
|
| 241 |
// results in Julian dates |
| 242 |
$data = array( |
| 243 |
$this->truephase( $k1, 0.0 ), |
| 244 |
$this->truephase( $k1, 0.25 ), |
| 245 |
$this->truephase( $k1, 0.5 ), |
| 246 |
$this->truephase( $k1, 0.75 ), |
| 247 |
$this->truephase( $k2, 0.0 ) |
| 248 |
); |
| 249 |
|
| 250 |
$this->quarters = array(); |
| 251 |
foreach( $data as $v ) |
| 252 |
$this->quarters[] = ( $v - 2440587.5 ) * 86400; // convert to UNIX time |
| 253 |
} |
| 254 |
|
| 255 |
/* Convert UNIX timestamp to astronomical Julian time (i.e. Julian date plus day fraction). */ |
| 256 |
private function utctojulian( $ts ) { |
| 257 |
return $ts / 86400 + 2440587.5; |
| 258 |
} |
| 259 |
|
| 260 |
private function get_phase( $n ) { |
| 261 |
if( is_null( $this->quarters ) ) |
| 262 |
$this->phasehunt(); |
| 263 |
|
| 264 |
return $this->quarters[$n]; |
| 265 |
} |
| 266 |
|
| 267 |
/* Public functions for accessing results */ |
| 268 |
|
| 269 |
function phase(){ |
| 270 |
return $this->phase; |
| 271 |
} |
| 272 |
|
| 273 |
function illumination(){ |
| 274 |
return $this->illum; |
| 275 |
} |
| 276 |
|
| 277 |
function age(){ |
| 278 |
return $this->age; |
| 279 |
} |
| 280 |
|
| 281 |
function distance(){ |
| 282 |
return $this->dist; |
| 283 |
} |
| 284 |
|
| 285 |
function diameter(){ |
| 286 |
return $this->angdia; |
| 287 |
} |
| 288 |
|
| 289 |
function sundistance(){ |
| 290 |
return $this->sundist; |
| 291 |
} |
| 292 |
|
| 293 |
function sundiameter(){ |
| 294 |
return $this->sunangdia; |
| 295 |
} |
| 296 |
|
| 297 |
function new_moon(){ |
| 298 |
return $this->get_phase( 0 ); |
| 299 |
} |
| 300 |
|
| 301 |
function first_quarter(){ |
| 302 |
return $this->get_phase( 1 ); |
| 303 |
} |
| 304 |
|
| 305 |
function full_moon(){ |
| 306 |
return $this->get_phase( 2 ); |
| 307 |
} |
| 308 |
|
| 309 |
function last_quarter(){ |
| 310 |
return $this->get_phase( 3 ); |
| 311 |
} |
| 312 |
|
| 313 |
function next_new_moon(){ |
| 314 |
return $this->get_phase( 4 ); |
| 315 |
} |
| 316 |
|
| 317 |
function phase_name() { |
| 318 |
$names = array( 'New Moon', 'Waxing Crescent', 'First Quarter', 'Waxing Gibbous', 'Full Moon', 'Waning Gibbous', 'Third Quarter', 'Waning Crescent', 'New Moon' ); |
| 319 |
// There are eight phases, evenly split. A "New Moon" occupies the 1/16th phases either side of phase = 0, and the rest follow from that. |
| 320 |
return $names[ floor( ( $this->phase + 0.0625 ) * 8 ) ]; |
| 321 |
} |
| 322 |
} |
