Code posted by: elf on 2013-09-29 17:30:33 GMT+10 - type PHP
  1. // various PHP functions related to determining positions of stars/dsos, planets, sun + moon based off a given time / lat / lon. Used for creationova.com. Some functions have been written from scratch, many by 3rd parties, some transcribed from other languages. Use at your own risk, ymmv!
  2.  
  3.  
  4. bcscale(20);
  5.  
  6. $DEGS = bcdiv('180',pi());
  7. $RADS = bcdiv(pi(),'180');
  8. $EPS = 1.0e-12;
  9. $K1 = 15*$RADS*1.0027379;
  10. $pname = array('Mercury', 'Venus', 'Sun', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto');
  11.  
  12. function phpsun($time, $lat, $lon, $offset)
  13. {
  14.         $phpsun = '';
  15.         $tzh = substr($offset,1,2);
  16.         $tzm = substr($offset,-2);
  17.         $tzpm = substr($offset,0,1);
  18.         $tzmin = $tzh*60 + $tzm;
  19.  
  20.  
  21.         $offset = $tzmin/60;
  22.         if ($tzpm == '-')
  23.                 $offset = -$offset;
  24.  
  25.         $zenith=108;
  26.         $phpsun[] = date_sunrise($time, SUNFUNCS_RET_STRING, $lat, $lon, $zenith, $offset);
  27.         $phpsun[] = date_sunset($time, SUNFUNCS_RET_STRING, $lat, $lon, $zenith, $offset);
  28.  
  29.         return $phpsun;
  30. }
  31.  
  32.  
  33. function coord_to_horizon($time, $ra, $dec, $lat, $lon, $h)
  34. {
  35.         //compute horizon (altaz) coordinates from ra/dec/lat/long + the given utc.
  36.  
  37.         $ha = bcsub(mean_sidereal_time($time, $lon), $ra);
  38.         if ($ha < 0)
  39.                 $ha = bcadd($ha, '360');
  40.  
  41.  
  42.         //convert from degrees to radians
  43.         $ha = bcdiv(bcmul($ha, pi()),180);
  44.         $dec = bcdiv(bcmul($dec, pi()),180);
  45.         $lat = bcdiv(bcmul($lat, pi()),180);
  46.  
  47.         //compute alt in radians
  48.         $sin_alt = bcadd(bcmul(sin($dec), sin($lat)), bcmul(bcmul(cos($dec), cos($lat)), cos($ha)));
  49.         $alt = asin($sin_alt);
  50.  
  51.         //compute azimuth in radians // does not cater for div by 0 errors (pole or alt90)
  52.         $cos_az = bcdiv(bcsub(sin($dec), bcmul(sin($alt), sin($lat))), bcmul(cos($alt), cos($lat)));
  53.         $az = acos($cos_az);
  54.  
  55.         $h['alt'] = bcdiv(bcmul($alt, '180'), pi());
  56.         $h['az'] = bcdiv(bcmul($az, '180'), pi());
  57.  
  58.         if (sin($ha) > 0)
  59.                 $h['az'] = bcsub('360', $h['az']);
  60.  
  61.         return $h;
  62. }
  63.  
  64.  
  65. function mean_sidereal_time($time, $longitude)
  66. {
  67.         //returns the Mean Sidereal Time in units of degrees.
  68.         //longitude of 0 = UTC aka GMT meridian.
  69.         //east is positive, west is negative.
  70.         $utc_year = gmdate('Y', $time);
  71.         $utc_month = gmdate('m', $time);
  72.         $utc_day = gmdate('d', $time);
  73.         $utc_hour = gmdate('H', $time);
  74.         $utc_minute = gmdate('i', $time);
  75.         $utc_second = gmdate('s', $time);
  76.  
  77.         //echo "$utc_year $utc_month $utc_day $utc_hour $utc_minute $utc_second";
  78.  
  79.         if ($utc_month == 1 || $utc_month == 2)
  80.         {
  81.                 $utc_year = bcsub($utc_year, '1');
  82.                 $utc_month = bcadd($utc_month, '12');
  83.         }
  84.  
  85.         $a = floor($utc_year/100);
  86.         $b = 2 - $a + floor($a/4);
  87.         $c = floor(365.25 * $utc_year);
  88.         $d = floor(30.6001 * ($utc_month + 1));
  89.  
  90.         $jd = $b + $c + $d - 730550.5 + $utc_day + ($utc_hour + $utc_minute/60 + $utc_second/3600)/24;
  91.  
  92.         $jt = $jd/36525;
  93.  
  94.         //mean sidereal time
  95.         $mst = 280.46061837 + 360.98564736629*$jd + 0.000387933*$jt*$jt - $jt*$jt*$jt/38710000 + $longitude;
  96.  
  97.         if ($mst > 0)
  98.         {
  99.                 while ($mst > 360)
  100.                         $mst -= 360;
  101.         }
  102.         else
  103.         {
  104.                 while ($mst < 0)
  105.                         $mst += 360;
  106.         }
  107.         return $mst;
  108. }
  109.  
  110.  
  111. function ra2real($hours, $minutes)
  112. {
  113.         //converts RA (hours, minutes) to degrees as real/float. 1 hour corresponds to 1/24th of a circle which equals 15 degrees. (useful for right ascension / longitude)
  114.         if (($hours >= 0 && $hours <= 24) && ($minutes >= 0 && $minutes <= 60))
  115.         {
  116.                 if ($minutes !== 0)
  117.                         $output = bcmul(15, bcadd($hours, bcdiv($minutes, '60')));
  118.                 else
  119.                         $output = bcmul(15, bcadd($hours, '0'));
  120.  
  121.                 return $output;
  122.         }
  123.         else
  124.                 return false;
  125. }
  126. function dms2real($degrees, $minutes, $seconds)
  127. {
  128.         //converts angle (degrees, minutes seconds) to degrees as real/float (useful for declination aka latitude)
  129.         if (($seconds >= 0 && $seconds <= 60) && ($minutes >= 0 && $minutes <= 60))
  130.         {
  131.                 if ($seconds === 0)
  132.                         $secondssum = 0;
  133.                 else
  134.                         $secondssum = $seconds / 3600;
  135.                 if ($minutes === 0)
  136.                         $minutessum = 0;
  137.                 else
  138.                         $minutessum = $minutes / 3600;
  139.  
  140.                 if ($degrees < 0)
  141.                 {
  142.  
  143.                         $output = $degrees - $minutessum - $secondssum;
  144.                 }
  145.                 else
  146.                 {
  147.                         //$output = bcadd($degrees, bcadd(bcdiv($minutes, '60'), bcdiv($seconds, '3600')));
  148.                         $output = $degrees + $minutessum + $secondssum;
  149.                 }
  150.                 //echo $output.'<br>';
  151.                 return $output;
  152.         }
  153.         else
  154.                 return false;
  155. }
  156.  
  157. function dms2real2($degrees, $minutes, $seconds)
  158. {
  159.         //echo "dms2real2 working with $degrees $minutes $seconds<br>";
  160.         if ($degrees < 0)
  161.                 $output = $degrees - $minutes/60 - $seconds/3600;
  162.         else
  163.                 $output = $degrees + $minutes/60 + $seconds/3600;
  164.  
  165.         return $output;
  166. }
  167.  
  168. // compute the true anomaly from mean anomaly using iteration
  169. //  M - mean anomaly in radians
  170. //  e - orbit eccentricity
  171. function true_anomaly($M, $e)
  172. {
  173.         $E = $M + $e*sin($M)*(1.0 + $e*cos($M));
  174.         do
  175.         {
  176.                 $E1 = $E;
  177.                 $E = bcsub($E1, bcdiv((bcsub($E1, bcmul($e,sin($E1))) - $M),(bcsub('1', bcmul($e,cos($E1))))));
  178.  
  179.         } while (abs(bcsub($E, $E1)) > $EPS);
  180.  
  181.         $V = 2*atan(sqrt((1 + $e)/(1 - $e))*tan(0.5*$E));
  182.  
  183.         if ($V < 0)
  184.                 $V = $V + (2*pi());
  185.  
  186.         return $V;
  187.  
  188. }
  189.  
  190. // converts hour angle in degrees into hour angle string
  191. function ha2str($x)
  192. {
  193.         if (($x < 0) || (360 < $x))
  194.                 echo 'ha2str range error!';
  195.  
  196.         $ra = $x / 15;
  197.         $h = floor($ra);
  198.         $m = 60*($ra - $h);
  199.         $output = cintstr($h, 3).'h '.frealstr($m, 4, 1) .'m';
  200.         return $output;
  201. }
  202.  
  203. // converts declination angle in degrees into string
  204. function dec2str($x)
  205. {
  206.         if (($x < -90) || (90 < $x))
  207.                 echo "dec2str range error! ($x)";
  208.  
  209.         $dec = abs($x);
  210.         $sgn = ($x < 0) ? -1 : +1;
  211.         $d = floor($dec);
  212.         $m = 60*($dec - $d);
  213.         $output = cintstr($sgn*$d, 3).'&deg; '.frealstr($m, 4, 1)."'";
  214.         return $output;
  215. }
  216.  
  217. // return the integer part of a number
  218. function abs_floor($x)
  219. {
  220.         $r = round($x,0);
  221.  
  222.         return $r;
  223. }
  224.  
  225. // return an angle in the range 0 to 2pi radians
  226. function mod2pi($x)
  227. {
  228.         $b = bcdiv($x,(bcmul('2',pi())));
  229.         $a = bcmul(bcmul('2',pi()),bcsub($b, abs_floor($b)));
  230.  
  231.         if ($a < 0)
  232.                 $a = bcadd(bcmul('2',pi()), $a);
  233.  
  234.         return $a;
  235. }
  236.  
  237. // converts angle in degrees into string
  238. function degr2str($x)
  239. {
  240.         $dec = abs($x);
  241.         $sgn = ($x < 0) ? -1 : +1;
  242.         $d = floor($dec);
  243.         $m = 60 * ($dec - $d);
  244.         $output = cintstr($sgn*$d, 3) .'&deg; '.frealstr($m, 4, 1)."'";
  245.         return $output;
  246. }
  247.  
  248. // converts latitude in signed degrees into string
  249. function lat2str($x, $method)
  250. {
  251.         $dec = abs($x);
  252.         $sgn = ($x < 0) ? 'S' : 'N';
  253.         $d = floor($dec);
  254.         $mfloor = 60 * ($dec - $d);
  255.         $m = floor($mfloor);
  256.         $s = floor(60 * ($mfloor - $m));
  257.         if ($method == 'string')
  258.         {
  259.                 $output = cintstr($d, 3) .'&deg; '.$m."' ".$s.'" '.$sgn;
  260.                 return $output;
  261.         }
  262.         else
  263.         {
  264.                 $output = array('d' => $d, 'm' => $m, 's' => $s, 'sgn' => $sgn);
  265.                 return $output;
  266.         }
  267.  
  268. }
  269.  
  270. // converts longitude in signed degrees into string
  271. function lon2str($x, $method)
  272. {
  273.         $dec = abs($x); //take -153.123 to 153.123 (removes negativity)
  274.         $sgn = ($x < 0) ? 'W' : 'E'; //West or East depending on whether the string is positive or negative
  275.         $d = floor($dec); // floor the degrees to get an integer
  276.         $mfloor = 60 * ($dec - $d); // take the remainder, multiply it by 60 to get the minutes (real)
  277.         $m = floor($mfloor); //floor the minutes to get an integer
  278.         $s = floor(60 * ($mfloor - $m));
  279.         if ($method == 'string')
  280.         {
  281.                 $output = cintstr($d, 3) .'&deg; '.$m."' ".$s.'" '.$sgn;
  282.                 return $output;
  283.         }
  284.         else
  285.         {
  286.                 $output = array('d' => $d, 'm' => $m, 's' => $s, 'sgn' => $sgn);
  287.                 return $output;
  288.         }
  289. }
  290.  
  291. // format two digits with leading zero if needed
  292. function d2($n)
  293. {
  294.         if ($n < 10)
  295.                 $n = "0$n";
  296.  
  297.         return $n;
  298. }
  299.  
  300. // format an integer
  301. function cintstr($num, $width)
  302. {
  303.         $numlen = strlen($num);
  304.         $diff = $width - $numlen;
  305.         if ($diff > 0)
  306.         {
  307.                 for ($i=0; $i < $diff; $i++)
  308.                         $num = ' '.$num;
  309.         }
  310.         return $num;
  311. }
  312.  
  313. function frealstr($num, $width, $fract)
  314. {
  315.         $num = number_format($num, $fract);
  316.  
  317.         $numlen = strlen($num);
  318.         $diff = $width - $numlen;
  319.         if ($diff > 0)
  320.         {
  321.                 for ($i=0; $i < $diff; $i++)
  322.                         $num = ' '.$num;
  323.         }
  324.  
  325.         return $num;
  326.  
  327. }
  328.  
  329. //day number to/from J2000
  330. function day_number($utc_year, $utc_month, $utc_day, $utc_hour, $utc_minute)
  331. {
  332.         global $pname, $DEGS, $RADS, $EPS;
  333.         if ($utc_minute !== 0)
  334.                 $h = bcadd($utc_hour, bcdiv($utc_minute,'60'));
  335.         else
  336.                 $h = bcadd($utc_hour, $utc_minute);
  337.  
  338.         $rv = bcadd(bcsub(bcadd(bcadd(bcsub(bcmul('367',$utc_year), floor(bcmul('7', bcdiv((bcadd($utc_year, floor(bcdiv(bcadd($utc_month,'9'),'12')))),'4')))), floor(bcmul('275',bcdiv($utc_month,'9')))), $utc_day), '730531.5'), bcdiv($h,'24'));
  339.  
  340.         return $rv;
  341. }
  342.  
  343. // compute RA, DEC, and distance of planet-p for day number-d
  344. // result returned in structure obj in degrees and astronomical units
  345. function get_coord($obj, $p, $d)
  346. {
  347.         global $pname, $DEGS, $RADS, $EPS;
  348.         $planet = array('a' => '', 'e' => '', 'i' => '', 'O' => '', 'w' => '', 'L' => '');
  349.         $planet = mean_elements($planet, $p, $d);
  350.         $ap = $planet['a'];
  351.         $ep = $planet['e'];
  352.         $ip = $planet['i'];
  353.         $op = $planet['O'];
  354.         $pp = $planet['w'];
  355.         $lp = $planet['L'];
  356.  
  357.         $earth = array('a' => '', 'e' => '', 'i' => '', 'O' => '', 'w' => '', 'L' => '');
  358.         $earth = mean_elements($earth, 2, $d);
  359.         $ae = $earth['a'];
  360.         $ee = $earth['e'];
  361.         $ie = $earth['i'];
  362.         $oe = $earth['O'];
  363.         $pe = $earth['w'];
  364.         $le = $earth['L'];
  365.  
  366.         //position of Earth in its orbit
  367.         $me = mod2pi($le - $pe);
  368.         $ve = true_anomaly($me, $ee);
  369.         $re = $ae*(1 - $ee*$ee)/(1 + $ee*cos($ve));
  370.  
  371.         // heliocentric rectangular coordinates of Earth
  372.         $xe = $re*cos($ve+$pe);
  373.         $ye = $re*sin($ve+$pe);
  374.         $ze = '0.0';
  375.  
  376.         // position of planet in its orbit
  377.         $mp = mod2pi(bcsub($lp, $pp));
  378.         $vp = true_anomaly($mp, $planet['e']);
  379.         $rp = $ap*(1 - $ep*$ep)/(1 + $ep*cos($vp));
  380.  
  381.         // heliocentric rectangular coordinates of planet
  382.         $xh = $rp*(cos($op)*cos($vp+$pp-$op) - sin($op)*sin($vp+$pp-$op)*cos($ip));
  383.         $yh = $rp*(sin($op)*cos($vp+$pp-$op) + cos($op)*sin($vp+$pp-$op)*cos($ip));
  384.         $zh = $rp*(sin($vp+$pp-$op)*sin($ip));
  385.  
  386.         //earth --> compute sun
  387.         if($p == 2)
  388.         {
  389.                 $xh = 0;
  390.                 $yh = 0;
  391.                 $zh = 0;
  392.         }
  393.  
  394.         // convert to geocentric rectangular coordinates
  395.         $xg = $xh - $xe;
  396.         $yg = $yh - $ye;
  397.         $zg = $zh - $ze;
  398.  
  399.         // rotate around x axis from ecliptic to equatorial coords
  400.         $ecl = 23.439281*$RADS;
  401.         $xeq = $xg;
  402.         $yeq = $yg*cos($ecl) - $zg*sin($ecl);
  403.         $zeq = $yg*sin($ecl) + $zg*cos($ecl);
  404.  
  405.         $obj['ra'] = mod2pi(atan2($yeq, $xeq))*$DEGS;
  406.         $obj['dec'] = atan($zeq/sqrt($xeq*$xeq + $yeq*$yeq))*$DEGS;
  407.         $obj['rvec'] = sqrt($xeq*$xeq + $yeq*$yeq + $zeq*$zeq);
  408.  
  409.         return $obj;
  410.  
  411. }
  412.  
  413. // Compute the elements of the orbit for planet-i at day number-d
  414. // result is returned in structure p
  415. function mean_elements($p, $i, $d)
  416. {
  417.  
  418.         global $pname, $DEGS, $RADS, $EPS;
  419.         $cy = bcdiv($d,'36525');
  420.  
  421.         switch($i)
  422.         {
  423.                 case 0: //mercury
  424.                         $p['a'] = bcadd('0.38709893', bcmul('0.00000066',$cy));
  425.                         $p['e'] = bcadd('0.20563069', bcmul('0.00002527',$cy));
  426.                         $p['i'] = bcmul(bcsub('7.00487', bcmul('23.51',bcdiv($cy,'3600'))),$RADS);
  427.                         $p['O'] = bcmul((bcsub('48.33167', bcmul('446.30',bcdiv($cy,'3600')))),$RADS);
  428.                         $p['w'] = bcmul(bcadd('77.45645', bcmul('573.57',bcdiv($cy,'3600'))),$RADS);
  429.                         $p['L'] = mod2pi(bcmul(bcadd('252.25084',bcmul('538101628.29',bcdiv($cy,'3600'))),$RADS));
  430.                         break;
  431.                 case 1: //venus
  432.                         $p['a'] = bcadd('0.72333199', bcmul('0.00000092',$cy));
  433.                         $p['e'] = bcsub('0.00677323', bcmul('0.00004938',$cy));
  434.                         $p['i'] = bcmul(bcsub('3.39471', bcmul('2.86',bcdiv($cy,'3600'))),$RADS);
  435.                         $p['O'] = bcmul(bcsub('76.68069', bcmul('996.89',bcdiv($cy,'3600'))),$RADS);
  436.                         $p['w'] = bcmul(bcsub('131.53298', bcmul('108.80',bcdiv($cy,'3600'))),$RADS);
  437.                         $p['L'] = mod2pi(bcmul(bcadd('181.97973', bcmul('210664136.06',bcdiv($cy,'3600'))),$RADS));
  438.                         break;
  439.                 case 2: //earth/sun
  440.                         $p['a'] = bcsub('1.00000011', bcmul('0.00000005',$cy));
  441.                         $p['e'] = bcsub('0.01671022', bcmul('0.00003804',$cy));
  442.                         $p['i'] = bcmul(bcsub('0.00005', bcmul('46.94',bcdiv($cy,'3600'))),$RADS);
  443.                         $p['O'] = bcmul(bcsub('-11.26064', bcmul('18228.25',bcdiv($cy,'3600'))),$RADS);
  444.                         $p['w'] = bcmul(bcadd('102.94719', bcmul('1198.28',bcdiv($cy,'3600'))),$RADS);
  445.                         $p['L'] = mod2pi(bcmul(bcadd('100.46435', bcmul('129597740.63',bcdiv($cy,'3600'))),$RADS));
  446.                         break;
  447.                 case 3: //mars
  448.                                         $p['a'] = bcsub('1.52366231', bcmul('0.00007221',$cy));
  449.                         $p['e'] = bcadd('0.09341233', bcmul('0.00011902',$cy));
  450.                         $p['i'] = bcmul(bcsub('1.85061', bcmul('25.47',bcdiv($cy,'3600'))),$RADS);
  451.                         $p['O'] = bcmul(bcsub('49.57854', bcmul('1020.19',bcdiv($cy,'3600'))),$RADS);
  452.                         $p['w'] = bcmul(bcadd('336.04084', bcmul('1560.78',bcdiv($cy,'3600'))),$RADS);
  453.                         $p['L'] = mod2pi(bcmul(bcadd('355.45332', bcmul('68905103.78',bcdiv($cy,'3600'))),$RADS));
  454.                         break;
  455.                 case 4: //jupiter
  456.                         $p['a'] = bcadd('5.20336301', bcmul('0.00060737',$cy));
  457.                         $p['e'] = bcsub('0.04839266', bcmul('0.00012880',$cy));
  458.                         $p['i'] = bcmul(bcsub('1.30530', bcmul('4.15',bcdiv($cy,'3600'))),$RADS);
  459.                         $p['O'] = bcmul(bcadd('100.55615', bcmul('1217.17',bcdiv($cy,'3600'))),$RADS);
  460.                         $p['w'] = bcmul(bcadd('14.75385', bcmul('839.93',bcdiv($cy,'3600'))),$RADS);
  461.                         $p['L'] = mod2pi(bcmul(bcadd('34.40438', bcmul('10925078.35',bcdiv($cy,'3600'))),$RADS));
  462.                         break;
  463.                 case 5: //saturn
  464.                         $p['a'] = bcsub('9.53707032', bcmul('0.00301530',$cy));
  465.                         $p['e'] = bcsub('0.05415060', bcmul('0.00036762',$cy));
  466.                         $p['i'] = bcmul(bcadd('2.48446', bcmul('6.11',bcdiv($cy,'3600'))),$RADS);
  467.                         $p['O'] = bcmul(bcsub('113.71504', bcmul('1591.05',bcdiv($cy,'3600'))),$RADS);
  468.                         $p['w'] = bcmul(bcsub('92.43194', bcmul('1948.89',bcdiv($cy,'3600'))),$RADS);
  469.                         $p['L'] = mod2pi(bcmul(bcadd('49.94432', bcmul('4401052.95',bcdiv($cy,'3600'))),$RADS));
  470.                         break;
  471.                 case 6: //uranus
  472.                         $p['a'] = bcadd('19.19126393', bcmul('0.00152025',$cy));
  473.                         $p['e'] = bcsub('0.04716771', bcmul('0.00019150',$cy));
  474.                         $p['i'] = bcmul(bcsub('0.76986', bcmul('2.09',bcdiv($cy,'3600'))),$RADS);
  475.                         $p['O'] = bcmul(bcsub('74.22988', bcmul('1681.40',bcdiv($cy,'3600'))),$RADS);
  476.                         $p['w'] = bcmul(bcadd('170.96424', bcmul('1312.56',bcdiv($cy,'3600'))),$RADS);
  477.                         $p['L'] = mod2pi(bcmul(bcadd('313.23218', bcmul('1542547.79',bcdiv($cy,'3600'))),$RADS));
  478.                         break;
  479.                 case 7: //neptune
  480.                         $p['a'] = bcsub('30.06896348', bcmul('0.00125196',$cy));
  481.                         $p['e'] = bcadd('0.00858587', bcmul('0.00002510',$cy));
  482.                         $p['i'] = bcmul(bcsub('1.76917', bcmul('3.64',bcdiv($cy,'3600'))),$RADS);
  483.                         $p['O'] = bcmul(bcsub('131.72169', bcmul('151.25',bcdiv($cy,'3600'))),$RADS);
  484.                         $p['w'] = bcmul(bcsub('44.97135', bcmul('844.43',bcdiv($cy,'3600'))),$RADS);
  485.                         $p['L'] = mod2pi(bcmul(bcadd('304.88003', bcmul('786449.21',bcdiv($cy,'3600'))),$RADS));
  486.                         break;
  487.                 case 8: //pluto
  488.                         $p['a'] = bcsub('39.48168677', bcmul('0.00076912',$cy));
  489.                         $p['e'] = bcadd('0.24880766', bcmul('0.00006465',$cy));
  490.                         $p['i'] = bcmul(bcadd('17.14175', bcmul('11.07',bcdiv($cy,'3600'))),$RADS);
  491.                         $p['O'] = bcmul(bcsub('110.30347', bcmul('37.33',bcdiv($cy,'3600'))),$RADS);
  492.                         $p['w'] = bcmul(bcsub('224.06676', bcmul('132.25',bcdiv($cy,'3600'))),$RADS);
  493.                         $p['L'] = mod2pi(bcmul(bcadd('238.92881', bcmul('522747.90',bcdiv($cy,'3600'))),$RADS));
  494.                         break;
  495.                 default:
  496.                         echo "mean_elements failed!";
  497.         }
  498.  
  499.         return $p;
  500. }
  501. function compute($time, $lat, $lon)
  502. {
  503.         global $pname, $DEGS, $RADS, $EPS;
  504.  
  505.         //$lat = dms2real2($lat_degrees, $lat_minutes, 0);
  506.         //$lon = dms2real2($lon_degrees, $lon_minutes, 0);
  507.  
  508.         $utc_year = gmdate('Y', $time);
  509.         $utc_month = gmdate('m', $time);
  510.         $utc_day = gmdate('d', $time);
  511.         $utc_hour = gmdate('H', $time);
  512.         $utc_minute = gmdate('i', $time);
  513.         $utc_second = gmdate('s', $time);
  514.  
  515.         $dn = day_number($utc_year, $utc_month, $utc_day, $utc_hour, $utc_minute);
  516.  
  517.     // compute location of objects
  518.     foreach ($pname AS $key => $value)
  519.         {
  520.                 if ($key != 2)
  521.                 {
  522.                         $h = array('alt' => '', 'az' => '');
  523.                         $obj = array('ra' => '', 'dec' => '', 'rvec' => '');
  524.                         $obj1 = get_coord($obj, $key, $dn);
  525.  
  526.                         $obj2 = coord_to_horizon($time, $obj1['ra'], $obj1['dec'], $lat, $lon, $h);
  527.                  
  528.                         $output[] = array('object' => $value, 'ra' => $obj1['ra'], 'dec' => $obj1['dec'], 'alt' => $obj2['alt'], 'az' => $obj2['az']);
  529.                 }
  530.  
  531.         }
  532.         return $output;
  533. }
  534.  
  535. class Moon
  536. {
  537. function phase($Year, $Month, $Day, $Hour, $Minutes, $Seconds)
  538.         {
  539.         $DateSec = mktime($Hour, $Minutes, $Seconds, $Month, $Day, $Year, 0);
  540.  
  541.         ini_set(precision, "20");       //Defini la precision des calcules
  542.  
  543.         # Astronomical constants.
  544.        $Epoch                                  = 2444238.5;            # 1980 January 0.0
  545.  
  546.         # Constants defining the Sun's apparent orbit.
  547.        $Elonge                                 = 278.833540;           # ecliptic longitude of the Sun at epoch 1980.0
  548.        $Elongp                                 = 282.596403;           # ecliptic longitude of the Sun at perigee
  549.        $Eccent                                 = 0.016718;                     # eccentricity of Earth's orbit
  550.        $Sunsmax                                = 1.495985e8;           # semi-major axis of Earth's orbit, km
  551.        $Sunangsiz                              = 0.533128;                     # sun's angular size, degrees, at semi-major axis distance
  552.  
  553.         # Elements of the Moon's orbit, epoch 1980.0.
  554.        $Mmlong                                 = 64.975464;            # moon's mean longitude at the epoch
  555.        $Mmlongp                                = 349.383063;           # mean longitude of the perigee at the epoch
  556.        $Mlnode                                 = 151.950429;           # mean longitude of the node at the epoch
  557.        $Minc                                   = 5.145396;                     # inclination of the Moon's orbit
  558.        $Mecc                                   = 0.054900;                     # eccentricity of the Moon's orbit
  559.        $Mangsiz                                = 0.5181;                       # moon's angular size at distance a from Earth
  560.        $Msmax                                  = 384401.0;                     # semi-major axis of Moon's orbit in km
  561.        $Mparallax                              = 0.9507;                       # parallax at distance a from Earth
  562.        $Synmonth                               = 29.53058868;          # synodic month (new Moon to new Moon)
  563.  
  564.         $pdate = Moon::jtime($DateSec);
  565.  
  566.         $pphase;                                # illuminated fraction
  567.        $mage;                                  # age of moon in days
  568.        $dist;                                  # distance in kilometres
  569.        $angdia;                                # angular diameter in degrees
  570.        $sudist;                                # distance to Sun
  571.        $suangdia;                              # sun's angular diameter
  572.  
  573.  
  574.         # Calculation of the Sun's position.
  575.  
  576.         $Day = $pdate - $Epoch;                                                                         # date within epoch
  577.        $N = Moon::fixangle((360 / 365.2422) * $Day);                           # mean anomaly of the Sun
  578.        $M = Moon::fixangle($N + $Elonge - $Elongp);                            # convert from perigee
  579.                                                                        # co-ordinates to epoch 1980.0
  580.        $Ec = Moon::kepler($M, $Eccent);                                                        # solve equation of Kepler
  581.        $Ec = sqrt((1 + $Eccent) / (1 - $Eccent)) * tan($Ec / 2);
  582.         $Ec = 2 * Moon::todeg(atan($Ec));                                                       # true anomaly
  583.        $Lambdasun = Moon::fixangle($Ec + $Elongp);                                     # Sun's geocentric ecliptic
  584.                                                                        # longitude
  585.        # Orbital distance factor.
  586.        $F = ((1 + $Eccent * cos(Moon::torad($Ec))) / (1 - $Eccent * $Eccent));
  587.         $SunDist = $Sunsmax / $F;                                                                       # distance to Sun in km
  588.        $SunAng = $F * $Sunangsiz;                                                                      # Sun's angular size in degrees
  589.  
  590.  
  591.         # Calculation of the Moon's position.
  592.  
  593.         # Moon's mean longitude.
  594.        $ml = Moon::fixangle(13.1763966 * $Day + $Mmlong);
  595.  
  596.         # Moon's mean anomaly.
  597.        $MM = Moon::fixangle($ml - 0.1114041 * $Day - $Mmlongp);
  598.  
  599.         # Moon's ascending node mean longitude.
  600.        $MN = Moon::fixangle($Mlnode - 0.0529539 * $Day);
  601.  
  602.         # Evection.
  603.        $Ev = 1.2739 * sin(Moon::torad(2 * ($ml - $Lambdasun) - $MM));
  604.  
  605.         # Annual equation.
  606.        $Ae = 0.1858 * sin(Moon::torad($M));
  607.  
  608.         # Correction term.
  609.        $A3 = 0.37 * sin(Moon::torad($M));
  610.  
  611.         # Corrected anomaly.
  612.        $MmP = $MM + $Ev - $Ae - $A3;
  613.  
  614.         # Correction for the equation of the centre.
  615.        $mEc = 6.2886 * sin(Moon::torad($MmP));
  616.  
  617.         # Another correction term.
  618.        $A4 = 0.214 * sin(Moon::torad(2 * $MmP));
  619.  
  620.         # Corrected longitude.
  621.        $lP = $ml + $Ev + $mEc - $Ae + $A4;
  622.  
  623.         # Variation.
  624.        $V = 0.6583 * sin(Moon::torad(2 * ($lP - $Lambdasun)));
  625.  
  626.         # True longitude.
  627.        $lPP = $lP + $V;
  628.  
  629.         # Corrected longitude of the node.
  630.        $NP = $MN - 0.16 * sin(Moon::torad($M));
  631.  
  632.         # Y inclination coordinate.
  633.        $y = sin(Moon::torad($lPP - $NP)) * cos(Moon::torad($Minc));
  634.  
  635.         # X inclination coordinate.
  636.        $x = cos(Moon::torad($lPP - $NP));
  637.  
  638.         # Ecliptic longitude.
  639.        $Lambdamoon = Moon::todeg(atan2($y, $x));
  640.         $Lambdamoon += $NP;
  641.  
  642.         # Ecliptic latitude.
  643.        $BetaM = Moon::todeg(asin(sin(Moon::torad($lPP - $NP)) * sin(Moon::torad($Minc))));
  644.  
  645.         # Calculation of the phase of the Moon.
  646.  
  647.         # Age of the Moon in degrees.
  648.        $MoonAge = $lPP - $Lambdasun;
  649.  
  650.         # Phase of the Moon.
  651.        $MoonPhase = (1 - cos(Moon::torad($MoonAge))) / 2;
  652.  
  653.         # Calculate distance of moon from the centre of the Earth.
  654.  
  655.         $MoonDist = ($Msmax * (1 - $Mecc * $Mecc)) /
  656.                 (1 + $Mecc * cos(Moon::torad($MmP + $mEc)));
  657.  
  658.         # Calculate Moon's angular diameter.
  659.  
  660.         $MoonDFrac = $MoonDist / $Msmax;
  661.         $MoonAng = $Mangsiz / $MoonDFrac;
  662.  
  663.         # Calculate Moon's parallax.
  664.  
  665.         $MoonPar = $Mparallax / $MoonDFrac;
  666.  
  667.         $pphase = $MoonPhase;                                                                   # illuminated fraction
  668.        $mage = $Synmonth * (Moon::fixangle($MoonAge) / 360.0); # age of moon in days
  669.        $dist = $MoonDist;                                                                              # distance in kilometres
  670.        $angdia = $MoonAng;                                                                             # angular diameter in degrees
  671.        $sudist = $SunDist;                                                                             # distance to Sun
  672.        $suangdia = $SunAng;                                                                    # sun's angular diameter
  673.        $mpfrac = Moon::fixangle($MoonAge) / 360.0;
  674.         return array( $pphase, $mage, $dist, $angdia, $sudist, $suangdia, $mpfrac, $mpfrac );
  675.         }
  676.  
  677. function fixangle($x)   { return ($x - 360.0 * (floor($x / 360.0))); }  # fix angle
  678. function torad($x)      { return ($x * (M_PI / 180.0)); }                               # deg->rad
  679. function todeg($x)      { return ($x * (180.0 / M_PI)); }                               # rad->deg
  680.  
  681. function jtime($t)
  682. {
  683.         $julian = ($t / 86400) + 2440587.5;     # (seconds /(seconds per day)) + julian date of epoch           2440587.5 / 86400 = 28,24753472222 Days
  684.        return ($julian);
  685. }
  686.  
  687. function kepler($m, $ecc)
  688. {
  689.         $EPSILON = 1e-6;
  690.  
  691.         $m = Moon::torad($m);
  692.         $e = $m;
  693.         while (abs($delta) > $EPSILON)
  694.                 {
  695.                 $delta = $e - $ecc * sin($e) - $m;
  696.                 $e -= $delta / (1 - $ecc * cos($e));
  697.                 }
  698.         return ($e);
  699.         }
  700.  
  701. }
  702.  
  703. function computesun($time, $timezone, $lat, $lon)
  704. {
  705.         $utc_year = gmdate('Y', $time);
  706.         $utc_month = gmdate('m', $time);
  707.         $utc_day = gmdate('d', $time);
  708.         $utc_hour = gmdate('H', $time);
  709.         $utc_minute = gmdate('i', $time);
  710.         $utc_second = gmdate('s', $time);
  711.  
  712.         $output = sunriseset($time, $timezone, $lat, $lon);
  713.  
  714.         return $output;
  715. }
  716.  
  717. function computemoon($time, $timezone, $lat, $lon)
  718. {
  719.         $utc_year = gmdate('Y', $time);
  720.         $utc_month = gmdate('m', $time);
  721.         $utc_day = gmdate('d', $time);
  722.         $utc_hour = gmdate('H', $time);
  723.         $utc_minute = gmdate('i', $time);
  724.         $utc_second = gmdate('s', $time);
  725.  
  726.     $output = riseset($time, $timezone, $lat, $lon);
  727.  
  728.         return $output;
  729. }
  730.  
  731. function sunriseset($time, $timezone, $lat, $lon)
  732. {
  733.  
  734.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Sunrise, $Sunset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  735.         $tzh = substr($timezone,1,2);
  736.         $tzm = substr($timezone,-2);
  737.         $tzpm = substr($timezone,0,1);
  738.         $tzmin = $tzh*60 + $tzm;
  739.  
  740.         $zone = $tzmin/60;
  741.         if ($tzpm == '+')
  742.                 $zone = -$zone;
  743.  
  744.         $jd = julian_day($time) - 2451545;           // Julian day relative to Jan 1.5, 2000
  745.  
  746.         if ((sgn($zone) == sgn($lon))&&($zone != 0))
  747.         echo "<br><font color='red'>WARNING: time zone and longitude are incompatible!</font><br>";
  748.  
  749.         $Sky = array('0.0', '0.0', '0.0');
  750.         $RAn = array('0.0', '0.0', '0.0');
  751.         $Dec = array('0.0', '0.0', '0.0');
  752.         $VHz = array('0.0', '0.0', '0.0');
  753.  
  754.         $mp = array(array('0.0','0.0','0.0'), array('0.0','0.0','0.0'), array('0.0','0.0','0.0'));
  755.  
  756.         $lon = $lon/360;
  757.         $tz  = $zone/24;
  758.         $ct  = $jd/36525 + 1;                    // centuries since 1900.0
  759.         $t0 = lst($lon, $jd, $tz);                 // local sidereal time
  760.  
  761.         $jd = $jd + $tz;                              // get sun position at start of day
  762.         sun($jd, $ct);
  763.         $ra0  = $Sky[0];
  764.         $dec0 = $Sky[1];
  765.  
  766.         $jd = $jd + 1;                               // get sun position at end of day
  767.         sun($jd, $ct);
  768.         $ra1  = $Sky[0];
  769.         $dec1 = $Sky[1];
  770.  
  771.         if ($ra1 <$ $ra0)                             // make continuous
  772.                 $ra1 = $ra1 + 2*pi();
  773.  
  774.         $Sunrise = false;                           // initialize
  775.         $Sunset  = false;
  776.         $RAn[0]   = $ra0;
  777.         $Dec[0]  = $dec0;
  778.  
  779.         for ($k = 0; $k < 24; $k++)                   // check each hour of this day
  780.         {
  781.                 $ph = ($k + 1)/24;
  782.  
  783.                 $RAn[2] = $ra0  + ($k + 1)*($ra1  - $ra0)/24;
  784.                 $Dec[2] = $dec0 + ($k + 1)*($dec1 - $dec0)/24;
  785.                 $VHz[2] = test_sun($k, $zone, $t0, $lat);
  786.  
  787.                 $RAn[0] = $RAn[2];                       // advance to next hour
  788.                 $Dec[0] = $Dec[2];
  789.                 $VHz[0] = $VHz[2];
  790.             }
  791.  
  792.         // results
  793.         $output = array('upallday' => 0, 'downallday' => 0, 'rises' => 0, 'sets' => 0, 'risetime' => 0, 'riseaz' => 0, 'settime' => 0, 'setaz' => 0);
  794.  
  795.  
  796.         // neither rise nor set
  797.         if ((!$Sunrise)&&(!$Sunset))
  798.         {
  799.                 if ($VHz[2] < 0)
  800.                         $output['downallday'] = 1;
  801.                 else
  802.                         $output['upallday'] = 1;
  803.  
  804.         }
  805.  
  806.         if($Sunrise)
  807.         {
  808.                 $output['rises'] = 1;
  809.                 $output['risetime'] = d2($Rise_time[0]).':'.d2($Rise_time[1]);
  810.                 $output['riseaz'] = frealstr($Rise_az, 5, 1);
  811.         }
  812.         if($Sunset)
  813.         {
  814.                 $output['sets'] = 1;
  815.                 $output['settime'] = d2($Set_time[0]).':'.d2($Set_time[1]);
  816.                 $output['setaz'] = frealstr($Set_az, 5, 1);
  817.         }
  818.  
  819.         return $output;
  820. }
  821.  
  822. // calculate moonrise and moonset times
  823. function riseset($time, $timezone, $lat, $lon )
  824. {
  825.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Sunrise, $Sunset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  826.         //need number of hours between UTC and user's timezone, expressed as UTC in relation to the tz, eg if in gmt+10 (+1000), then we need "-10".
  827.         $tzh = substr($timezone,1,2);
  828.         $tzm = substr($timezone,-2);
  829.         $tzpm = substr($timezone,0,1);
  830.         $tzmin = $tzh*60 + $tzm;
  831.  
  832.         $zone = $tzmin/60;
  833.         if ($tzpm == '+')
  834.                 $zone = -$zone;
  835.  
  836.         // Julian day relative to Jan 1.5, 2000
  837.         $jd = julian_day($time) - 2451545;
  838.  
  839.         if ((sgn($zone) == sgn($lon)) && ($zone !== 0))
  840.                 echo "<br><font color='red'>WARNING: time zone and longitude are incompatible!</font><br>";
  841.  
  842.         $Sky = array('0.0', '0.0', '0.0');
  843.         $RAn = array('0.0', '0.0', '0.0');
  844.         $Dec = array('0.0', '0.0', '0.0');
  845.         $VHz = array('0.0', '0.0', '0.0');
  846.  
  847.         $mp = array(array('0.0','0.0','0.0'), array('0.0','0.0','0.0'), array('0.0','0.0','0.0'));
  848.  
  849.         $lon = $lon/360;
  850.         //echo "<br>lon $lon";
  851.  
  852.         $tz = $zone/24;
  853.         // local sidereal time
  854.         $t0 = lst($lon, $jd, $tz);
  855.  
  856.         // get moon position at start of day
  857.         $jd = $jd + $tz;
  858.  
  859.         //echo "<br>t0 $t0 jd $jd";
  860.         for ($k = 0; $k <3; $k++)
  861.         {
  862.                 moon($jd);
  863.                 $mp[$k][0] = $Sky[0];
  864.                 $mp[$k][1] = $Sky[1];
  865.                 $mp[$k][2] = $Sky[2];
  866.                 $jd = $jd + 0.5;
  867.  
  868.         }
  869.  
  870.         if($mp[1][0] <= $mp[0][0])
  871.                 $mp[1][0] = $mp[1][0] + 2 * pi();
  872.  
  873.         if($mp[2][0] <= $mp[1][0])
  874.                 $mp[2][0] = $mp[2][0] + 2 * pi();
  875.  
  876.         $RAn[0] = $mp[0][0];
  877.         $Dec[0] = $mp[0][1];
  878.  
  879.         // check each hour of this day
  880.         for ($k = 0; $k < 24; $k++)
  881.         {
  882.                 $ph = ($k + 1) / 24;
  883.  
  884.                 $RAn[2] = interpolate($mp[0][0], $mp[1][0], $mp[2][0], $ph);
  885.                 $Dec[2] = interpolate($mp[0][1], $mp[1][1], $mp[2][1], $ph);
  886.  
  887.                 $VHz[2] = test_moon($k, $zone, $t0, $lat, $mp[1][2]);
  888.  
  889.                 // advance to next hour
  890.                 $RAn[0] = $RAn[2];
  891.                 $Dec[0] = $Dec[2];
  892.                 $VHz[0] = $VHz[2];
  893.         }
  894.  
  895.         // results
  896.         $output = array('upallday' => 0, 'downallday' => 0, 'rises' => 0, 'sets' => 0, 'risetime' => 0, 'riseaz' => 0, 'settime' => 0, 'setaz' => 0);
  897.  
  898.         // neither rise nor set
  899.         if ((!$Moonrise)&&(!$Moonset))
  900.         {
  901.                 if ($VHz[2] < 0)
  902.                         $output['downallday'] = 1;
  903.                 else
  904.                         $output['upallday'] = 1;
  905.  
  906.         }
  907.  
  908.         if($Moonrise)
  909.         {
  910.                 $output['rises'] = 1;
  911.                 $output['risetime'] = d2($Rise_time[0]).':'.d2($Rise_time[1]);
  912.                 $output['riseaz'] = frealstr($Rise_az, 5, 1);
  913.         }
  914.         if($Moonset)
  915.         {
  916.                 $output['sets'] = 1;
  917.                 $output['settime'] = d2($Set_time[0]).':'.d2($Set_time[1]);
  918.                 $output['setaz'] = frealstr($Set_az, 5, 1);
  919.         }
  920.  
  921.         return $output;
  922.  
  923.  
  924. }
  925.  
  926. // Local Sidereal Time for zone
  927. function lst($lon, $jd, $z)
  928. {
  929.     global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Sunrise, $Sunset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  930.         $s = 24110.5 + 8640184.812999999*$jd/36525 + 86636.6*$z + 86400*$lon;
  931.     $s = $s/86400;
  932.     $s = $s - floor($s);
  933.     $s = $s*360*$RADS;
  934.     return $s;
  935. }
  936.  
  937. // 3-point interpolation
  938. function interpolate( $f0, $f1, $f2, $p )
  939. {
  940.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  941.         $a = $f1 - $f0;
  942.         $b = $f2 - $f1 - $a;
  943.         $f = $f0 + $p*(2*$a + $b *(2*$p -1));
  944.  
  945.         //var f = f0 + p*(2*a + b*(2*p - 1));
  946.  
  947.         return $f;
  948. }
  949.  
  950.  
  951. function test_sun( $k, $zone, $t0, $lat)
  952. {
  953.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Sunrise, $Sunset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  954.    
  955.  
  956.         $ha[0] = $t0 - $RAn[0] + $k*$K1;
  957.         $ha[2] = $t0 - $RAn[2] + $k*$K1 + $K1;
  958.  
  959.         $ha[1]  = ($ha[2]  + $ha[0])/2;               // hour angle at half hour
  960.         $Dec[1] = ($Dec[2] + $Dec[0])/2 ;             // declination at half hour
  961.  
  962.         $s = sin($lat*$RADS);
  963.         $c = cos($lat*$RADS);
  964.         $z = cos(90.833*$RADS);                   // refraction + sun semidiameter at horizon
  965.  
  966.         if ($k <= 0)
  967.                 $VHz[0] = $s*sin($Dec[0]) + $c*cos($Dec[0])*cos($ha[0]) - $z;
  968.  
  969.         $VHz[2] = $s*sin($Dec[2]) + $c*cos($Dec[2])*cos($ha[2]) - $z;
  970.  
  971.         if (sgn($VHz[0]) == sgn($VHz[2]))
  972.                 return $VHz[2];                         // no event this hour
  973.  
  974.         $VHz[1] = $s*sin($Dec[1]) + $c*cos($Dec[1])*cos($ha[1]) - $z;
  975.  
  976.         $a =  2* $VHz[0] - 4*$VHz[1] + 2*$VHz[2];
  977.         $b = -3* $VHz[0] + 4*$VHz[1] - $VHz[2];
  978.         $d = $b*$b - 4*$a*$VHz[0];
  979.  
  980.         if ($d < 0)
  981.                 return $VHz[2];                         // no event this hour
  982.  
  983.         $d = sqrt($d);
  984.         $e = (-$b + $d)/(2 * $a);
  985.  
  986.         if (($e > 1)||($e < 0))
  987.                 $e = (-$b - $d)/(2*$a);
  988.  
  989.         $newtime = $k + $e + 1/120;                      // time of an event
  990.  
  991.         $hr = floor($newtime);
  992.         $min = floor(($newtime - $hr)*60);
  993.  
  994.         $hz = $ha[0] + $e*($ha[2] - $ha[0]);            // azimuth of the sun at the event
  995.         $nz = -cos($Dec[1])*sin($hz);
  996.         $dz = $c*sin($Dec[1]) - $s*cos($Dec[1])*cos($hz);
  997.         $az = atan2($nz, $dz)/$RADS;
  998.         if ($az < 0) $az = $az + 360;
  999.  
  1000.         if (($VHz[0] < 0)&&($VHz[2] > 0))
  1001.         {
  1002.                 $Rise_time[0] = $hr;
  1003.                 $Rise_time[1] = $min;
  1004.                 $Rise_az = $az;
  1005.                 $Sunrise = true;
  1006.         }
  1007.  
  1008.         if (($VHz[0] > 0)&&($VHz[2] < 0))
  1009.         {
  1010.                 $Set_time[0] = $hr;
  1011.                 $Set_time[1] = $min;
  1012.                 $Set_az = $az;
  1013.                 $Sunset = true;
  1014.         }
  1015.  
  1016.         return $VHz[2];
  1017. }
  1018.  
  1019. // test an hour for an event
  1020. function test_moon( $k, $zone, $t0, $lat, $plx )
  1021. {
  1022.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  1023.         $ha = array('0.0', '0.0', '0.0');
  1024.  
  1025.         if ($RAn[2] < $RAn[0])
  1026.                 $RAn[2] = $RAn[2] + 2* pi();
  1027.  
  1028.         $ha[0] = $t0 - $RAn[0] + $k*$K1;
  1029.         $ha[2] = $t0 - $RAn[2] + $k*$K1 + $K1;
  1030.  
  1031.         $ha[1]  = ($ha[2] + $ha[0])/2;                // hour angle at half hour
  1032.         $Dec[1] = ($Dec[2] + $Dec[0])/2;              // declination at half hour
  1033.  
  1034.         $s = sin($RADS * $lat);
  1035.         $c = cos($RADS * $lat);
  1036.  
  1037.         // refraction + sun semidiameter at horizon + parallax correction
  1038.         $z = cos($RADS*(90.567 - 41.685/$plx));
  1039.  
  1040.         // first call of function
  1041.         if ($k <= 0)
  1042.                 $VHz[0] = $s*sin($Dec[0]) + $c*cos($Dec[0])*cos($ha[0]) - $z;
  1043.  
  1044.         $VHz[2] = $s*sin($Dec[2]) + $c*cos($Dec[2])*cos($ha[2]) - $z;
  1045.  
  1046.         if (sgn($VHz[0]) == sgn($VHz[2]))
  1047.                 return $VHz[2]; // no event this hour
  1048.  
  1049.         $VHz[1] = $s*sin($Dec[1]) + $c*cos($Dec[1])*cos($ha[1]) - $z;
  1050.  
  1051.         $a = 2*$VHz[2] - 4*$VHz[1] + 2*$VHz[0];
  1052.         $b = 4*$VHz[1] - 3*$VHz[0] - $VHz[2];
  1053.         $d = $b*$b - 4*$a*$VHz[0];
  1054.  
  1055.         if ($d < 0)
  1056.                 return $VHz[2]; //no event this hour
  1057.  
  1058.         $d = sqrt($d);
  1059.         $e = (-$b + $d)/(2*$a);
  1060.  
  1061.         if (($e > 1) || ($e < 0))
  1062.                 $e = (-$b - $d)/(2*$a);
  1063.  
  1064.         // time of an event + round up
  1065.         $newtime = $k + $e + 1/120;
  1066.         $hr = floor($newtime);
  1067.         $min = floor(($newtime - $hr)*60);
  1068.  
  1069.         $hz = $ha[0] + $e*($ha[2] - $ha[0]);            // azimuth of the moon at the event
  1070.         $nz = -cos($Dec[1])*sin($hz);
  1071.         $dz = $c*sin($Dec[1]) - $s*cos($Dec[1])*cos($hz);
  1072.         $az = atan2($nz, $dz)/$RADS;
  1073.         if ($az < 0) $az = $az + 360;
  1074.  
  1075.         if (($VHz[0] < 0)&&($VHz[2] > 0))
  1076.         {
  1077.                 $Rise_time[0] = $hr;
  1078.                 $Rise_time[1] = $min;
  1079.                 $Rise_az = $az;
  1080.                 $Moonrise = true;
  1081.         }
  1082.  
  1083.         if (($VHz[0] > 0)&&($VHz[2] < 0))
  1084.         {
  1085.                 $Set_time[0] = $hr;
  1086.                 $Set_time[1] = $min;
  1087.                 $Set_az = $az;
  1088.                 $Moonset = true;
  1089.         }
  1090.  
  1091.         return $VHz[2];
  1092. }
  1093.  
  1094. // moon's position using fundamental arguments
  1095. // (Van Flandern & Pulkkinen, 1979)
  1096. function moon($jd )
  1097. {
  1098.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  1099.  
  1100.         $h = 0.606434 + 0.03660110129*$jd;
  1101.         $m = 0.374897 + 0.03629164709*$jd;
  1102.         $f = 0.259091 + 0.0367481952 *$jd;
  1103.         $d = 0.827362 + 0.03386319198*$jd;
  1104.         $n = 0.347343 - 0.00014709391*$jd;
  1105.         $g = 0.993126 + 0.0027377785 *$jd;
  1106.  
  1107.  
  1108.         $h = $h - floor($h);
  1109.         $m = $m - floor($m);
  1110.         $f = $f - floor($f);
  1111.         $d = $d - floor($d);
  1112.         $n = $n - floor($n);
  1113.         $g = $g - floor($g);
  1114.  
  1115.  
  1116.         $h = $h*2*pi();
  1117.         $m = $m*2*pi();
  1118.         $f = $f*2*pi();
  1119.         $d = $d*2*pi();
  1120.         $n = $n*2*pi();
  1121.         $g = $g*2*pi();
  1122.  
  1123.  
  1124.         $v = 0.39558*sin($f + $n);
  1125.         $v = $v + 0.082  *sin($f);
  1126.         $v = $v + 0.03257*sin($m - $f - $n);
  1127.         $v = $v + 0.01092*sin($m + $f + $n);
  1128.         $v = $v + 0.00666*sin($m - $f);
  1129.         $v = $v - 0.00644*sin($m + $f - 2*$d + $n);
  1130.         $v = $v - 0.00331*sin($f - 2*$d + $n);
  1131.         $v = $v - 0.00304*sin($f - 2*$d);
  1132.         $v = $v - 0.0024 *sin($m - $f - 2*$d - $n);
  1133.         $v = $v + 0.00226*sin($m + $f);
  1134.         $v = $v - 0.00108*sin($m + $f - 2*$d);
  1135.         $v = $v - 0.00079*sin($f - $n);
  1136.         $v = $v + 0.00078*sin($f + 2*$d + $n);
  1137.  
  1138.  
  1139.         $u = 1 - 0.10828*cos($m);
  1140.         $u = $u - 0.0188 *cos($m - 2*$d);
  1141.         $u = $u - 0.01479*cos(2*$d);
  1142.         $u = $u + 0.00181*cos(2*$m - 2*$d);
  1143.         $u = $u - 0.00147*cos(2*$m);
  1144.         $u = $u - 0.00105*cos(2*$d - $g);
  1145.         $u = $u - 0.00075*cos($m - 2*$d + $g);
  1146.  
  1147.  
  1148.         $w = 0.10478*sin($m);
  1149.         $w = $w - 0.04105*sin(2*$f + 2*$n);
  1150.         $w = $w - 0.0213 *sin($m - 2*$d);
  1151.         $w = $w - 0.01779*sin(2*$f + $n);
  1152.         $w = $w + 0.01774*sin($n);
  1153.         $w = $w + 0.00987*sin(2*$d);
  1154.         $w = $w - 0.00338*sin($m - 2*$f - 2*$n);
  1155.         $w = $w - 0.00309*sin($g);
  1156.         $w = $w - 0.0019 *sin(2*$f);
  1157.         $w = $w - 0.00144*sin($m + $n);
  1158.         $w = $w - 0.00144*sin($m - 2*$f - $n);
  1159.         $w = $w - 0.00113*sin($m + 2*$f + 2*$n);
  1160.         $w = $w - 0.00094*sin($m - 2*$d + $g);
  1161.         $w = $w - 0.00092*sin(2*$m - 2*$d);
  1162.  
  1163.         $s = $w/sqrt($u - $v*$v);                  // compute moon's right ascension ...
  1164.         $Sky[0] = $h + atan($s/sqrt(1 - $s*$s));
  1165.  
  1166.         $s = $v/sqrt($u);                        // declination ...
  1167.         $Sky[1] = atan($s/sqrt(1 - $s*$s));
  1168.  
  1169.         $Sky[2] = 60.40974*sqrt( $u );          // and parallax
  1170. }
  1171.  
  1172.  
  1173. function sun( $jd, $ct )
  1174. {
  1175.         global $RADS, $K1, $Sky, $RAn, $Dec, $VHz, $Moonrise, $Moonset, $Sunrise, $Sunset, $Rise_time, $Set_time, $Rise_az, $Set_az;
  1176.  
  1177.         $lo = 0.779072 + 0.00273790931*$jd;
  1178.         $lo = $lo - floor($lo);
  1179.         $lo = $lo*2*pi();
  1180.  
  1181.         $g = 0.993126 + 0.0027377785*$jd;
  1182.         $g = $g - floor($g);
  1183.         $g = $g*2*pi();
  1184.  
  1185.         $v = 0.39785*sin($lo);
  1186.         $v = $v - 0.01*sin($lo - $g);
  1187.         $v = $v + 0.00333*sin($lo + $g);
  1188.         $v = $v - 0.00021*$ct * sin($lo);
  1189.  
  1190.         $u = 1 - 0.03349*cos($g);
  1191.         $u = $u - 0.00014*cos(2*$lo);
  1192.         $u = $u + 0.00008*cos($lo);
  1193.  
  1194.         $w = -0.0001 - 0.04129*sin(2*$lo);
  1195.         $w = $w + 0.03211*sin($g );
  1196.         $w = $w + 0.00104*sin(2*$lo - $g);
  1197.         $w = $w - 0.00035*sin(2*$lo + $g);
  1198.         $w = $w - 0.00008*$ct*sin($g);
  1199.  
  1200.         $s = $w/sqrt($u - $v*$v);                  // compute sun's right ascension
  1201.         $Sky[0] = $lo + atan($s/sqrt(1 - $s*$s));
  1202.  
  1203.         $s = $v/sqrt($u);                        // ...and declination
  1204.         $Sky[1] = atan($s/sqrt(1 - $s*$s));
  1205. }
  1206.  
  1207. // determine Julian day from calendar date
  1208. // (Jean Meeus, "Astronomical Algorithms", Willmann-Bell, 1991)
  1209. function julian_day($time)
  1210. {
  1211.  
  1212.         $month = date('m', $time);
  1213.         $day = date('d', $time);
  1214.         $year = date('Y', $time);
  1215.  
  1216.         $gregorian = ($year < 1583) ? false : true;
  1217.  
  1218.         if (($month == 1)||($month == 2))
  1219.         {
  1220.                 $year  = $year  - 1;
  1221.                 $month = $month + 12;
  1222.         }
  1223.  
  1224.         $a = floor($year/100);
  1225.         if ($gregorian) $b = 2 - $a + floor($a/4);
  1226.         else           $b = 0.0;
  1227.  
  1228.         $jd = floor(365.25*($year + 4716)) + floor(30.6001*($month + 1)) + $day + $b - 1524.5;
  1229.  
  1230.         return $jd;
  1231. }
  1232.  
  1233. // returns value for sign of argument
  1234. function sgn( $x )
  1235. {
  1236.         if ($x > 0.0)      $rv =  1;
  1237.         else if ($x < 0.0) $rv = -1;
  1238.         else              $rv =  0;
  1239.         return $rv;
  1240. }
  1241.  
  1242. function direction($x)
  1243. {
  1244.         if ($x > 337.5 || $x <= 22.5)
  1245.                 $direction = 'north';
  1246.         if ($x > 22.5 && $x <= 67.5)
  1247.                 $direction = 'northeast';
  1248.         if ($x > 67.5 && $x <= 112.5)
  1249.                 $direction = 'east';
  1250.         if ($x > 112.5 && $x <= 157.5)
  1251.                 $direction = 'southeast';
  1252.         if ($x > 157.5 && $x <= 202.5)
  1253.                 $direction = 'south';
  1254.         if ($x > 202.5 && $x <= 247.5)
  1255.                 $direction = 'southwest';
  1256.         if ($x > 247.5 && $x <= 292.5)
  1257.                 $direction = 'west';
  1258.         if ($x > 292.5 && $x <= 337.5)
  1259.                 $direction = 'northwest';
  1260.  
  1261.         return $direction;
  1262. }
  1263.  
  1264. function timezoneoffsetarray()
  1265. {
  1266.         $timezone_identifiers = DateTimeZone::listIdentifiers();
  1267.         foreach ($timezone_identifiers AS $value) {
  1268.                 $thistimezone = new DateTimeZone($value);
  1269.                 $thistime = new DateTime("now", $thistimezone);
  1270.                 $tzarray[$value] = $thistime->getOffset();
  1271.         }
  1272.         return $tzarray;
  1273. }
  1274.  
  1275. function timezoneoffset($value)
  1276. {
  1277.         $thistimezone = new DateTimeZone($value);
  1278.         $thistime = new DateTime("now", $thistimezone);
  1279.         $tzarray[$value] = $thistime->getOffset();
  1280.         return $tzarray;
  1281. }
  1282.  
  1283. function timezoneoffsetfromme($value, $dateTimeZoneMe, $dateTimeMe)
  1284. {
  1285.         $thistimezone = new DateTimeZone($value);
  1286.         $thistime = new DateTime("now", $thistimezone);
  1287.         $offset = $thistimezone->getOffset($dateTimeMe) - $dateTimeZoneMe->getOffset($dateTimeMe);
  1288.         return $offset;
  1289. }
  1290.  
  1291. function isdst($value)
  1292. {
  1293.         $thistimezone = new DateTimeZone($value);
  1294.         $thistime = new DateTime("now", $thistimezone);
  1295.         if ($thistime->format('I') == 1)
  1296.                 return true;
  1297.         else
  1298.                 return false;
  1299. }
  1300.  
  1301. function getoffset($tzarray)
  1302. {
  1303.         //expects (key)'Australia/Brisbane', (value)'36000'
  1304.         foreach($tzarray as $key => $value)
  1305.         {
  1306.                 $isnegative = 0;
  1307.                 if ($value <0)
  1308.                         $isnegative = 1;
  1309.                 $value2 = abs($value);
  1310.  
  1311.                 if ($value2 % 3600 != 0)
  1312.                         $minutc = ($value2 % 3600) / 60;
  1313.                 else
  1314.                         $minutc = '00';
  1315.  
  1316.                         $hourutc = floor($value2/3600);
  1317.                 //echo "hourutc = $hourutc";
  1318.  
  1319.                 $utc = $hourutc.$minutc;
  1320.                 if ($utc > -999 && $utc < 999)
  1321.                 {
  1322.                         //echo "matches padding";
  1323.                         $utc = str_pad($utc, 4, "0", STR_PAD_LEFT);
  1324.                 }
  1325.  
  1326.                 if ($isnegative == 0)
  1327.                         $utc = '+'.$utc;
  1328.                 else
  1329.                         $utc = '-'.$utc;
  1330.  
  1331.                 return $utc;
  1332.         }
  1333.  
  1334. }
  1335.  
  1336.  

New post:
Respond to the original code above (post ID vxdxv24)? Yes (click here to paste the code from above), or click here to make a new post.

Syntax highlighting?
Name (blank for Anonymous)
How long to keep it for?
 
Type the above security code:


© 2006-2008 CreatioNova.com and CreationExchange.com.
Site design and code mostly developed by elf (with the exception of GeSHi of course).
Hits: 67041