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

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!