rad function todeg($x) { return ($x * (180.0 / M_PI)); } # rad->deg function jtime($t) { $julian = ($t / 86400) + 2440587.5; # (seconds /(seconds per day)) + julian date of epoch 2440587.5 / 86400 = 28,24753472222 Days return ($julian); } function kepler($m, $ecc) { $EPSILON = 1e-6; $m = Lune::torad($m); $e = $m; while (abs($delta) > $EPSILON) { $delta = $e - $ecc * sin($e) - $m; $e -= $delta / (1 - $ecc * cos($e)); } return ($e); } /** * Calculates the moon rise/set for a given location and day of year */ public static function calculateMoonTimes($month, $day, $year, $lat, $lon) { $utrise = $utset = 0; $timezone = (int)($lon / 15); $date = self::modifiedJulianDate($month, $day, $year); $date -= $timezone / 24; $latRad = deg2rad($lat); $sinho = 0.0023271056; $sglat = sin($latRad); $cglat = cos($latRad); $rise = false; $set = false; $above = false; $hour = 1; $ym = self::sinAlt($date, $hour - 1, $lon, $cglat, $sglat) - $sinho; $above = $ym > 0; while ($hour < 25 && (false == $set || false == $rise)) { $yz = self::sinAlt($date, $hour, $lon, $cglat, $sglat) - $sinho; $yp = self::sinAlt($date, $hour + 1, $lon, $cglat, $sglat) - $sinho; $quadout = self::quad($ym, $yz, $yp); $nz = $quadout[0]; $z1 = $quadout[1]; $z2 = $quadout[2]; $xe = $quadout[3]; $ye = $quadout[4]; if ($nz == 1) { if ($ym < 0) { $utrise = $hour + $z1; $rise = true; } else { $utset = $hour + $z1; $set = true; } } if ($nz == 2) { if ($ye < 0) { $utrise = $hour + $z2; $utset = $hour + $z1; } else { $utrise = $hour + $z1; $utset = $hour + $z2; } } $ym = $yp; $hour += 2.0; } // Convert to unix timestamps and return as an object $retVal = new stdClass(); $utrise = self::convertTime($utrise); $utset = self::convertTime($utset); $retVal->moonrise = $rise ? mktime($utrise['hrs'], $utrise['min'], 0, $month, $day, $year) : mktime(0, 0, 0, $month, $day + 1, $year); $retVal->moonset = $set ? mktime($utset['hrs'], $utset['min'], 0, $month, $day, $year) : mktime(0, 0, 0, $month, $day + 1, $year); return $retVal; } /** * finds the parabola throuh the three points (-1,ym), (0,yz), (1, yp) * and returns the coordinates of the max/min (if any) xe, ye * the values of x where the parabola crosses zero (roots of the self::quadratic) * and the number of roots (0, 1 or 2) within the interval [-1, 1] * * well, this routine is producing sensible answers * * results passed as array [nz, z1, z2, xe, ye] */ private static function quad($ym, $yz, $yp) { $nz = $z1 = $z2 = 0; $a = 0.5 * ($ym + $yp) - $yz; $b = 0.5 * ($yp - $ym); $c = $yz; $xe = -$b / (2 * $a); $ye = ($a * $xe + $b) * $xe + $c; $dis = $b * $b - 4 * $a * $c; if ($dis > 0) { $dx = 0.5 * sqrt($dis) / abs($a); $z1 = $xe - $dx; $z2 = $xe + $dx; $nz = abs($z1) < 1 ? $nz + 1 : $nz; $nz = abs($z2) < 1 ? $nz + 1 : $nz; $z1 = $z1 < -1 ? $z2 : $z1; } return array($nz, $z1, $z2, $xe, $ye); } /** * this rather mickey mouse function takes a lot of * arguments and then returns the sine of the altitude of the moon */ private static function sinAlt($mjd, $hour, $glon, $cglat, $sglat) { $mjd += $hour / 24; $t = ($mjd - 51544.5) / 36525; $objpos = self::minimoon($t); $ra = $objpos[1]; $dec = $objpos[0]; $decRad = deg2rad($dec); $tau = 15 * (self::lmst($mjd, $glon) - $ra); return $sglat * sin($decRad) + $cglat * cos($decRad) * cos(deg2rad($tau)); } /** * returns an angle in degrees in the range 0 to 360 */ private static function degRange($x) { $b = $x / 360; $a = 360 * ($b - (int)$b); $retVal = $a < 0 ? $a + 360 : $a; return $retVal; } private static function lmst($mjd, $glon) { $d = $mjd - 51544.5; $t = $d / 36525; $lst = self::degRange(280.46061839 + 360.98564736629 * $d + 0.000387933 * $t * $t - $t * $t * $t / 38710000); return $lst / 15 + $glon / 15; } /** * takes t and returns the geocentric ra and dec in an array mooneq * claimed good to 5' (angle) in ra and 1' in dec * tallies with another approximate method and with ICE for a couple of dates */ private static function minimoon($t) { $p2 = 6.283185307; $arc = 206264.8062; $coseps = 0.91748; $sineps = 0.39778; $lo = self::frac(0.606433 + 1336.855225 * $t); $l = $p2 * self::frac(0.374897 + 1325.552410 * $t); $l2 = $l * 2; $ls = $p2 * self::frac(0.993133 + 99.997361 * $t); $d = $p2 * self::frac(0.827361 + 1236.853086 * $t); $d2 = $d * 2; $f = $p2 * self::frac(0.259086 + 1342.227825 * $t); $f2 = $f * 2; $sinls = sin($ls); $sinf2 = sin($f2); $dl = 22640 * sin($l); $dl += -4586 * sin($l - $d2); $dl += 2370 * sin($d2); $dl += 769 * sin($l2); $dl += -668 * $sinls; $dl += -412 * $sinf2; $dl += -212 * sin($l2 - $d2); $dl += -206 * sin ($l + $ls - $d2); $dl += 192 * sin($l + $d2); $dl += -165 * sin($ls - $d2); $dl += -125 * sin($d); $dl += -110 * sin($l + $ls); $dl += 148 * sin($l - $ls); $dl += -55 * sin($f2 - $d2); $s = $f + ($dl + 412 * $sinf2 + 541 * $sinls) / $arc; $h = $f - $d2; $n = -526 * sin($h); $n += 44 * sin($l + $h); $n += -31 * sin(-$l + $h); $n += -23 * sin($ls + $h); $n += 11 * sin(-$ls + $h); $n += -25 * sin(-$l2 + $f); $n += 21 * sin(-$l + $f); $L_moon = $p2 * self::frac($lo + $dl / 1296000); $B_moon = (18520.0 * sin($s) + $n) / $arc; $cb = cos($B_moon); $x = $cb * cos($L_moon); $v = $cb * sin($L_moon); $w = sin($B_moon); $y = $coseps * $v - $sineps * $w; $z = $sineps * $v + $coseps * $w; $rho = sqrt(1 - $z * $z); $dec = (360 / $p2) * atan($z / $rho); $ra = (48 / $p2) * atan($y / ($x + $rho)); $ra = $ra < 0 ? $ra + 24 : $ra; return array($dec, $ra); } /** * returns the self::fractional part of x as used in self::minimoon and minisun */ private static function frac($x) { $x -= (int)$x; return $x < 0 ? $x + 1 : $x; } /** * Takes the day, month, year and hours in the day and returns the * modified julian day number defined as mjd = jd - 2400000.5 * checked OK for Greg era dates - 26th Dec 02 */ private static function modifiedJulianDate($month, $day, $year) { if ($month <= 2) { $month += 12; $year--; } $a = 10000 * $year + 100 * $month + $day; $b = 0; if ($a <= 15821004.1) { $b = -2 * (int)(($year + 4716) / 4) - 1179; } else { $b = (int)($year / 400) - (int)($year / 100) + (int)($year / 4); } $a = 365 * $year - 679004; return $a + $b + (int)(30.6001 * ($month + 1)) + $day; } /** * Converts an hours decimal to hours and minutes */ private static function convertTime($hours) { $hrs = (int)($hours * 60 + 0.5) / 60.0; $h = (int)($hrs); $m = (int)(60 * ($hrs - $h) + 0.5); return array('hrs'=>$h, 'min'=>$m); } } ?>