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.  
Reply posted by: justhef@3drugs.com on 2017-10-27 22:23:16 GMT+10 (reply - 0 replies)
  1. Cialis 0 5  <a href=http://costofcial.com>online pharmacy</a> Newzealand Pharmacy Domperidone
Reply posted by: dannysilva121m@yahoo.com on 2017-11-13 15:58:04 GMT+10 (reply - 0 replies)
  1. Trying To Find The Best Diet Pill?
  2.  
  3. Trying to find the best diet pill may seem like an impossible task, especially with the multitude of diet pills available for purchase. Many people purchase a diet pill only to find out that the pill makes them feel jittery, nervous, or often has no effect at all.
  4.  
  5. Diet pills frequently contain the same or similar combination of ingredients and rarely contain anything new, innovative, or undiscovered to the supplement / weight loss industry. So, how can you find the best diet pill when most diet pills are made with similar ingredients?
  6.  
  7. One of the most common problems associated with taking diet pills is that the person taking the diet pill is uneducated about the dosage, effects, and promises offered as they relate to each diet pill. The research at website finds that there are three factors that should be taken into consideration when deciding to take a diet pill.
  8.  
  9. Dosage:
  10. It is important to take the pill exactly as recommended on the product label. Some people choose to increase the dosage thinking that the product will work faster or better. This is not the case, and many people become sick in response to the large dose. Reviewers at website often suggest that the recommended dosage be cut in half to give the body time to adjust to the stimulant in the diet pill. After the body has adjusted, it is fine to begin taking the regular dosage as recommended on the product label.
  11.  
  12. Effects:
  13. The effects listed on the product label are there because these are the effects that the product has had on 'some' of the test group. Some of the diet pill testers may be fine taking the product, while others may have adverse effects. The diet pill companies print this information to educate the buyer as well as to protect themselves from lawsuits. The consumer needs to read the label and educate themselves before taking the product. Many people who are sensitive to caffeine are surprised when the diet pill makes them feel nervous or nauseous, but this information is likely printed on the product, so with a little research these affects can be avoided.
  14.  
  15. Promises:
  16. If you read the fine print on product claims for diet pills and other weight loss supplements, you will see 'results not typical' printed very small somewhere where you are not expected to look. The diet pills advertised on television are responsible for some of the most outlandish claims. The results claimed in these advertisements are often unattainable within the given amount of time outlined in the ad. Don't expect to see results in two weeks like a lot of ads claim.
  17.  
  18. Wouldn't it be great if you could read reviews for diet pills from actual users of each diet pill? Diet Pill Reviews website has taken the trouble out of searching for the best diet pill. You can read reviews of over 150 of the most popular diet pills available.
  19.  
  20. Copyright 2006, Diet Pill Reviews  <a href=https://www.viagrapascherfr.com/>viagra pas cher</a>

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: 58518