daily pastebin goal
85%
SHARE
TWEET

PHP Moon rise/set

dxprog Jun 28th, 2011 3,637 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. <?php
  2.  
  3. /******************************************************************************
  4. * The following is a PHP implementation of the JavaScript code found at:      *
  5. * http://bodmas.org/astronomy/riset.html                                      *
  6. *                                                                             *
  7. * Original maths and code written by Keith Burnett <bodmas.org>               *
  8. * PHP port written by Matt "dxprog" Hackmann <dxprog.com>                     *
  9. *                                                                             *
  10. * This program is free software: you can redistribute it and/or modify        *
  11. * it under the terms of the GNU General Public License as published by        *
  12. * the Free Software Foundation, either version 3 of the License, or           *
  13. * (at your option) any later version.                                         *
  14. *                                                                             *
  15. * This program is distributed in the hope that it will be useful,             *
  16. * but WITHOUT ANY WARRANTY; without even the implied warranty of              *
  17. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               *
  18. * GNU General Public License for more details.                                *
  19. *                                                                             *
  20. * You should have received a copy of the GNU General Public License           *
  21. * along with this program.  If not, see <http://www.gnu.org/licenses/>.       *
  22. ******************************************************************************/
  23.  
  24. class Moon {
  25.        
  26.         /**
  27.          * Calculates the moon rise/set for a given location and day of year
  28.          */
  29.         public static function calculateMoonTimes($month, $day, $year, $lat, $lon) {
  30.        
  31.                 $utrise = $utset = 0;
  32.                
  33.                 $timezone = (int)($lon / 15);
  34.                 $date = self::modifiedJulianDate($month, $day, $year);
  35.                 $date -= $timezone / 24;
  36.                 $latRad = deg2rad($lat);
  37.                 $sinho = 0.0023271056;
  38.                 $sglat = sin($latRad);
  39.                 $cglat = cos($latRad);
  40.                
  41.                 $rise = false;
  42.                 $set = false;
  43.                 $above = false;
  44.                 $hour = 1;
  45.                 $ym = self::sinAlt($date, $hour - 1, $lon, $cglat, $sglat) - $sinho;
  46.                
  47.                 $above = $ym > 0;
  48.                 while ($hour < 25 && (false == $set || false == $rise)) {
  49.                
  50.                         $yz = self::sinAlt($date, $hour, $lon, $cglat, $sglat) - $sinho;
  51.                         $yp = self::sinAlt($date, $hour + 1, $lon, $cglat, $sglat) - $sinho;
  52.                        
  53.                         $quadout = self::quad($ym, $yz, $yp);
  54.                         $nz = $quadout[0];
  55.                         $z1 = $quadout[1];
  56.                         $z2 = $quadout[2];
  57.                         $xe = $quadout[3];
  58.                         $ye = $quadout[4];
  59.                        
  60.                         if ($nz == 1) {
  61.                                 if ($ym < 0) {
  62.                                         $utrise = $hour + $z1;
  63.                                         $rise = true;
  64.                                 } else {
  65.                                         $utset = $hour + $z1;
  66.                                         $set = true;
  67.                                 }
  68.                         }
  69.                        
  70.                         if ($nz == 2) {
  71.                                 if ($ye < 0) {
  72.                                         $utrise = $hour + $z2;
  73.                                         $utset = $hour + $z1;
  74.                                 } else {
  75.                                         $utrise = $hour + $z1;
  76.                                         $utset = $hour + $z2;
  77.                                 }
  78.                         }
  79.                        
  80.                         $ym = $yp;
  81.                         $hour += 2.0;
  82.                
  83.                 }
  84.                 // Convert to unix timestamps and return as an object
  85.                 $retVal = new stdClass();
  86.                 $utrise = self::convertTime($utrise);
  87.                 $utset = self::convertTime($utset);
  88.                 $retVal->moonrise = $rise ? mktime($utrise['hrs'], $utrise['min'], 0, $month, $day, $year) : mktime(0, 0, 0, $month, $day + 1, $year);
  89.                 $retVal->moonset = $set ? mktime($utset['hrs'], $utset['min'], 0, $month, $day, $year) : mktime(0, 0, 0, $month, $day + 1, $year);
  90.                 return $retVal;
  91.        
  92.         }
  93.        
  94.         /**
  95.          *      finds the parabola throuh the three points (-1,ym), (0,yz), (1, yp)
  96.          *  and returns the coordinates of the max/min (if any) xe, ye
  97.          *  the values of x where the parabola crosses zero (roots of the self::quadratic)
  98.          *  and the number of roots (0, 1 or 2) within the interval [-1, 1]
  99.          *
  100.          *      well, this routine is producing sensible answers
  101.          *
  102.          *  results passed as array [nz, z1, z2, xe, ye]
  103.          */
  104.         private static function quad($ym, $yz, $yp) {
  105.  
  106.                 $nz = $z1 = $z2 = 0;
  107.                 $a = 0.5 * ($ym + $yp) - $yz;
  108.                 $b = 0.5 * ($yp - $ym);
  109.                 $c = $yz;
  110.                 $xe = -$b / (2 * $a);
  111.                 $ye = ($a * $xe + $b) * $xe + $c;
  112.                 $dis = $b * $b - 4 * $a * $c;
  113.                 if ($dis > 0) {
  114.                         $dx = 0.5 * sqrt($dis) / abs($a);
  115.                         $z1 = $xe - $dx;
  116.                         $z2 = $xe + $dx;
  117.                         $nz = abs($z1) < 1 ? $nz + 1 : $nz;
  118.                         $nz = abs($z2) < 1 ? $nz + 1 : $nz;
  119.                         $z1 = $z1 < -1 ? $z2 : $z1;
  120.                 }
  121.  
  122.                 return array($nz, $z1, $z2, $xe, $ye);
  123.                
  124.         }
  125.  
  126.         /**
  127.          *      this rather mickey mouse function takes a lot of
  128.          *  arguments and then returns the sine of the altitude of the moon
  129.          */
  130.         private static function sinAlt($mjd, $hour, $glon, $cglat, $sglat) {
  131.                
  132.                 $mjd += $hour / 24;
  133.                 $t = ($mjd - 51544.5) / 36525;
  134.                 $objpos = self::minimoon($t);
  135.  
  136.                 $ra = $objpos[1];
  137.                 $dec = $objpos[0];
  138.                 $decRad = deg2rad($dec);
  139.                 $tau = 15 * (self::lmst($mjd, $glon) - $ra);
  140.  
  141.                 return $sglat * sin($decRad) + $cglat * cos($decRad) * cos(deg2rad($tau));
  142.  
  143.         }
  144.  
  145.         /**
  146.          *      returns an angle in degrees in the range 0 to 360
  147.          */
  148.         private static function degRange($x) {
  149.                 $b = $x / 360;
  150.                 $a = 360 * ($b - (int)$b);
  151.                 $retVal = $a < 0 ? $a + 360 : $a;
  152.                 return $retVal;
  153.         }
  154.  
  155.         private static function lmst($mjd, $glon) {
  156.                 $d = $mjd - 51544.5;
  157.                 $t = $d / 36525;
  158.                 $lst = self::degRange(280.46061839 + 360.98564736629 * $d + 0.000387933 * $t * $t - $t * $t * $t / 38710000);
  159.                 return $lst / 15 + $glon / 15;
  160.         }
  161.  
  162.         /**
  163.          * takes t and returns the geocentric ra and dec in an array mooneq
  164.          * claimed good to 5' (angle) in ra and 1' in dec
  165.          * tallies with another approximate method and with ICE for a couple of dates
  166.          */
  167.         private static function minimoon($t) {
  168.                        
  169.                 $p2 = 6.283185307;
  170.                 $arc = 206264.8062;
  171.                 $coseps = 0.91748;
  172.                 $sineps = 0.39778;
  173.                
  174.                 $lo = self::frac(0.606433 + 1336.855225 * $t);
  175.                 $l = $p2 * self::frac(0.374897 + 1325.552410 * $t);
  176.                 $l2 = $l * 2;
  177.                 $ls = $p2 * self::frac(0.993133 + 99.997361 * $t);
  178.                 $d = $p2 * self::frac(0.827361 + 1236.853086 * $t);
  179.                 $d2 = $d * 2;
  180.                 $f = $p2 * self::frac(0.259086 + 1342.227825 * $t);
  181.                 $f2 = $f * 2;
  182.                
  183.                 $sinls = sin($ls);
  184.                 $sinf2 = sin($f2);
  185.                
  186.                 $dl = 22640 * sin($l);
  187.                 $dl += -4586 * sin($l - $d2);
  188.                 $dl += 2370 * sin($d2);
  189.                 $dl += 769 * sin($l2);
  190.                 $dl += -668 * $sinls;
  191.                 $dl += -412 * $sinf2;
  192.                 $dl += -212 * sin($l2 - $d2);
  193.                 $dl += -206 * sin ($l + $ls - $d2);
  194.                 $dl += 192 * sin($l + $d2);
  195.                 $dl += -165 * sin($ls - $d2);
  196.                 $dl += -125 * sin($d);
  197.                 $dl += -110 * sin($l + $ls);
  198.                 $dl += 148 * sin($l - $ls);
  199.                 $dl += -55 * sin($f2 - $d2);
  200.                
  201.                 $s = $f + ($dl + 412 * $sinf2 + 541 * $sinls) / $arc;
  202.                 $h = $f - $d2;
  203.                 $n = -526 * sin($h);
  204.                 $n += 44 * sin($l + $h);
  205.                 $n += -31 * sin(-$l + $h);
  206.                 $n += -23 * sin($ls + $h);
  207.                 $n += 11 * sin(-$ls + $h);
  208.                 $n += -25 * sin(-$l2 + $f);
  209.                 $n += 21 * sin(-$l + $f);
  210.                
  211.                 $L_moon = $p2 * self::frac($lo + $dl / 1296000);
  212.                 $B_moon = (18520.0 * sin($s) + $n) / $arc;
  213.                
  214.                 $cb = cos($B_moon);
  215.                 $x = $cb * cos($L_moon);
  216.                 $v = $cb * sin($L_moon);
  217.                 $w = sin($B_moon);
  218.                 $y = $coseps * $v - $sineps * $w;
  219.                 $z = $sineps * $v + $coseps * $w;
  220.                 $rho = sqrt(1 - $z * $z);
  221.                 $dec = (360 / $p2) * atan($z / $rho);
  222.                 $ra = (48 / $p2) * atan($y / ($x + $rho));
  223.                 $ra = $ra < 0 ? $ra + 24 : $ra;
  224.                
  225.                 return array($dec, $ra);
  226.                
  227.         }
  228.  
  229.         /**
  230.          *      returns the self::fractional part of x as used in self::minimoon and minisun
  231.          */
  232.         private static function frac($x) {
  233.                 $x -= (int)$x;
  234.                 return $x < 0 ? $x + 1 : $x;
  235.         }
  236.  
  237.         /**
  238.          * Takes the day, month, year and hours in the day and returns the
  239.          * modified julian day number defined as mjd = jd - 2400000.5
  240.          * checked OK for Greg era dates - 26th Dec 02
  241.          */
  242.         private static function modifiedJulianDate($month, $day, $year) {
  243.                
  244.                 if ($month <= 2) {
  245.                         $month += 12;
  246.                         $year--;
  247.                 }
  248.                
  249.                 $a = 10000 * $year + 100 * $month + $day;
  250.                 $b = 0;
  251.                 if ($a <= 15821004.1) {
  252.                         $b = -2 * (int)(($year + 4716) / 4) - 1179;
  253.                 } else {
  254.                         $b = (int)($year / 400) - (int)($year / 100) + (int)($year / 4);
  255.                 }
  256.  
  257.                 $a = 365 * $year - 679004;
  258.                 return $a + $b + (int)(30.6001 * ($month + 1)) + $day;
  259.                
  260.         }
  261.  
  262.         /**
  263.          * Converts an hours decimal to hours and minutes
  264.          */
  265.         private static function convertTime($hours) {
  266.  
  267.                 $hrs = (int)($hours * 60 + 0.5) / 60.0;
  268.                 $h = (int)($hrs);
  269.                 $m = (int)(60 * ($hrs - $h) + 0.5);
  270.                 return array('hrs'=>$h, 'min'=>$m);
  271.                
  272.         }
  273. }
RAW Paste Data
Top