Posted By

tomas on 10/14/10


Tagged

php sun


Versions (?)

php sun position


 / Published in: PHP
 

  1. <?php
  2. // calculates the sun position and path throughout the day
  3. // input params: LAT&LON
  4. // output json
  5.  
  6. // not much commenting in code
  7. // ported from http://www.stjarnhimlen.se/comp/tutorial.html
  8. // MANY THANKS!
  9. // most var-names are identical to above tutorial
  10.  
  11.  
  12. $LAT = deg2rad($_GET["lat"]);
  13. $LON = deg2rad($_GET["lon"]);
  14.  
  15.  
  16. $year = gmdate("Y");
  17. $month = gmdate("m");
  18. $day = gmdate("d");
  19. $hour = gmdate("H") + (gmdate("i") / 60);
  20.  
  21. // get current position
  22. getsunpos($LAT, $LON, $year, $month, $day, $hour, $azimuth, $altitude, $sunstate);
  23.  
  24. // get 3 points during day to create circle
  25. getsunpos($LAT, $LON, $year, $month, $day, 1, $azm1, $alt1, $ignore);
  26. getsunpos($LAT, $LON, $year, $month, $day, 9, $azm2, $alt2, $ignore);
  27. getsunpos($LAT, $LON, $year, $month, $day, 13, $azm3, $alt3, $ignore);
  28.  
  29. // transform altitude from radians to number from 0 to 1
  30. $alt1 = ((pi()/2) - $alt1) / (pi()/2);
  31. $alt2 = ((pi()/2) - $alt2) / (pi()/2);
  32. $alt3 = ((pi()/2) - $alt3) / (pi()/2);
  33.  
  34. // polar to cartesian
  35. $x1 = cos($azm1) * $alt1; $y1 = sin($azm1) * $alt1;
  36. $x2 = cos($azm2) * $alt2; $y2 = sin($azm2) * $alt2;
  37. $x3 = cos($azm3) * $alt3; $y3 = sin($azm3) * $alt3;
  38.  
  39. findCentre($x1,$y1,$x2,$y2,$x3,$y3,$cx,$cy,$r);
  40.  
  41. echo "{\"result\":\"OK\",\"sunstate\":\"$sunstate\",\"azimuth\":$azimuth,\"altitude\":$altitude,";
  42. echo "\"centrex\":$cx,\"centrey\":$cy,\"radius\":$r}";
  43.  
  44.  
  45.  
  46. function getsunpos($LAT, $LON, $year, $month, $day, $hour, &$azimuth, &$altitude, &$sunstate) {
  47.  
  48. // julian date
  49. $d = 367*$year - floor((7*($year + floor(($month+9)/12)))/4)
  50. + floor((275*$month)/9) + $day - 730530;
  51.  
  52. $w = 4.9382415669097640822661983551248
  53. + .00000082193663128794959930855831205818* $d; // (longitude of perihelion)
  54. $a = 1.000000 ;// (mean distance, a.u.)
  55. $e = 0.016709 - .000000001151 * $d ;// (eccentricity)
  56. $M = 6.2141924418482506287494932704807
  57. + 0.017201969619332228715501852561964 * $d ;// (mean anomaly)
  58.  
  59.  
  60. $oblecl = 0.40909295936270689252387465029835
  61. - .0000000062186081248557962825791102081249 * $d ;// obliquity of the ecliptic
  62.  
  63. $L = $w + $M; // sun's mean longitude
  64.  
  65. $E = $M + $e * sin($M) * (1 + $e * cos($M));
  66.  
  67. $x = cos($E) - $e;
  68. $y = sin($E) * sqrt(1 - $e * $e);
  69.  
  70. $r = sqrt($x*$x + $y*$y);
  71. $v = atan2( $y, $x ) ;
  72.  
  73. $lon = $v + $w;
  74.  
  75. $x = $r * cos($lon);
  76. $y = $r * sin($lon);
  77. $z = 0.0;
  78.  
  79. $xequat = $x;
  80. $yequat = $y * cos($oblecl) + $z * sin($oblecl);
  81. $zequat = $y * sin($oblecl) + $z * cos($oblecl);
  82.  
  83. $RA = atan2( $yequat, $xequat );
  84. $Decl = asin( $zequat / $r );
  85.  
  86. $RAhours = r2d($RA)/15;
  87.  
  88. $GMST0 = r2d( $L + pi() ) / 15;//
  89. $SIDTIME = $GMST0 + $hour + rad2deg($LON)/15;
  90.  
  91. $HA = deg2rad(($SIDTIME - $RAhours) *15);
  92.  
  93. $x = cos($HA) * cos($Decl);
  94. $y = sin($HA) * cos($Decl);
  95. $z = sin($Decl);
  96.  
  97. $xhor = $x * cos(pi()/2 - $LAT) - $z * sin(pi()/2 - $LAT);
  98. $yhor = $y;
  99. $zhor = $x * sin(pi()/2 - $LAT) + $z * cos(pi()/2 - $LAT);
  100.  
  101. $azimuth = pi()*3/2 -atan2( $yhor, $xhor );
  102. $altitude = asin( $zhor );
  103.  
  104.  
  105. // check sunshine state
  106. $alt_d = rad2deg($altitude);
  107. if ($alt_d >= 0)
  108. $sunstate = "DAY";
  109. else if ($alt_d >= -6)
  110. $sunstate = "CIVIL-TWILIGHT";
  111. else if ($alt_d >= -12)
  112. $sunstate = "NAUTICAL-TWILIGHT";
  113. else if ($alt_d >= -18)
  114. $sunstate = "ASTRONOMICAL-TWILIGHT";
  115. else
  116. $sunstate = "NIGHT";
  117. }
  118.  
  119.  
  120. function r2d($r) {
  121. $d = rad2deg($r);
  122. while ($d<0) $d += 360;
  123. while ($d>360) $d -= 360;
  124. return $d;
  125. }
  126.  
  127.  
  128. // functions below are for the sun's path through the day
  129.  
  130. function findCentre($x1,$y1,$x2,$y2,$x3,$y3,&$cx,&$cy,&$r) {
  131.  
  132. if (!isPerpendicular($x1,$y1,$x2,$y2,$x3,$y3) ) CalcCircle($x1,$y1,$x2,$y2,$x3,$y3,$cx,$cy,$r);
  133. else if (!isPerpendicular($x1,$y1,$x3,$y3,$x2,$y2) ) CalcCircle($x1,$y1,$x3,$y3,$x2,$y2,$cx,$cy,$r);
  134. else if (!isPerpendicular($x2,$y2,$x1,$y1,$x3,$y3) ) CalcCircle($x2,$y2,$x1,$y1,$x3,$y3,$cx,$cy,$r);
  135. else if (!isPerpendicular($x2,$y2,$x3,$y3,$x1,$y1) ) CalcCircle($x2,$y2,$x3,$y3,$x1,$y1,$cx,$cy,$r);
  136. else if (!isPerpendicular($x3,$y3,$x2,$y2,$x1,$y1) ) CalcCircle($x3,$y3,$x2,$y2,$x1,$y1,$cx,$cy,$r);
  137. else if (!isPerpendicular($x3,$y3,$x1,$y1,$x2,$y2) ) CalcCircle($x3,$y3,$x1,$y1,$x2,$y2,$cx,$cy,$r);
  138. else {
  139. $cx=0;
  140. $cy=0;
  141. $r=0;
  142.  
  143. }
  144. }
  145.  
  146. function isPerpendicular($x1,$y1,$x2,$y2,$x3,$y3) {
  147. $m = 0.00000001;
  148. $dya = $y2 - $y1;
  149. $dxa = $x2 - $x1;
  150. $dyb = $y3 - $y2;
  151. $dxb = $x3 - $x2;
  152.  
  153. // checking whether the line of the two pts are vertical
  154. if (abs($dxa) <= $m && abs($dyb) <= $m){
  155. //TRACE("The points are pependicular and parallel to x-y axis\n");
  156. return false;
  157. }
  158.  
  159. if (abs($dxa) <= $m || abs($dxb) <= $m || abs($dya) < $m || abs($dyb) <= $m) {
  160. return true;
  161. }
  162.  
  163. }
  164.  
  165. function CalcCircle($x1,$y1,$x2,$y2,$x3,$y3,&$cx,&$cy,&$r) {
  166. $m = 0.00000001;
  167. $dya = $y2 - $y1;
  168. $dxa = $x2 - $x1;
  169. $dyb = $y3 - $y2;
  170. $dxb = $x3 - $x2;
  171.  
  172. if (abs($dxa) <= $m && abs($dxb) <= $m){
  173.  
  174. $cx= 0.5 * ($x2 + $x3);
  175. $cy= 0.5 * ($y1 + $y2);
  176. $r = sqrt(($x1-$cx)*($x1-$cx) + ($y1-$cy)*($y1-$cy));
  177. return;
  178. }
  179.  
  180. // IsPerpendicular() assure that xDelta(s) are not zero
  181. $aSlope = $dya/$dxa;
  182. $bSlope = $dyb/$dxb;
  183. if (abs($aSlope-$bSlope) <= 0.000000001){ // checking whether the given points are colinear.
  184.  
  185. return;
  186. }
  187.  
  188. // calc center
  189. $cx = ($aSlope*$bSlope*($y1 - $y3) + $bSlope*($x1 + $x2)
  190. - $aSlope*($x2+$x3) )/(2* ($bSlope-$aSlope) );
  191. $cy = -1*($cx - ($x1+$x2)/2)/$aSlope + ($y1+$y2)/2;
  192.  
  193. $r = sqrt(($x1-$cx)*($x1-$cx) + ($y1-$cy)*($y1-$cy));
  194. return;
  195.  
  196. }
  197.  
  198. ?>

Report this snippet  

You need to login to post a comment.