Revision: 33887
Initial Code
Initial URL
Initial Description
Initial Title
Initial Tags
Initial Language
at October 14, 2010 20:25 by tomas
Initial Code
<?php // calculates the sun position and path throughout the day // input params: LAT&LON // output json // not much commenting in code // ported from http://www.stjarnhimlen.se/comp/tutorial.html // MANY THANKS! // most var-names are identical to above tutorial $LAT = deg2rad($_GET["lat"]); $LON = deg2rad($_GET["lon"]); $year = gmdate("Y"); $month = gmdate("m"); $day = gmdate("d"); $hour = gmdate("H") + (gmdate("i") / 60); // get current position getsunpos($LAT, $LON, $year, $month, $day, $hour, $azimuth, $altitude, $sunstate); // get 3 points during day to create circle getsunpos($LAT, $LON, $year, $month, $day, 1, $azm1, $alt1, $ignore); getsunpos($LAT, $LON, $year, $month, $day, 9, $azm2, $alt2, $ignore); getsunpos($LAT, $LON, $year, $month, $day, 13, $azm3, $alt3, $ignore); // transform altitude from radians to number from 0 to 1 $alt1 = ((pi()/2) - $alt1) / (pi()/2); $alt2 = ((pi()/2) - $alt2) / (pi()/2); $alt3 = ((pi()/2) - $alt3) / (pi()/2); // polar to cartesian $x1 = cos($azm1) * $alt1; $y1 = sin($azm1) * $alt1; $x2 = cos($azm2) * $alt2; $y2 = sin($azm2) * $alt2; $x3 = cos($azm3) * $alt3; $y3 = sin($azm3) * $alt3; findCentre($x1,$y1,$x2,$y2,$x3,$y3,$cx,$cy,$r); echo "{\"result\":\"OK\",\"sunstate\":\"$sunstate\",\"azimuth\":$azimuth,\"altitude\":$altitude,"; echo "\"centrex\":$cx,\"centrey\":$cy,\"radius\":$r}"; function getsunpos($LAT, $LON, $year, $month, $day, $hour, &$azimuth, &$altitude, &$sunstate) { // julian date $d = 367*$year - floor((7*($year + floor(($month+9)/12)))/4) + floor((275*$month)/9) + $day - 730530; $w = 4.9382415669097640822661983551248 + .00000082193663128794959930855831205818* $d; // (longitude of perihelion) $a = 1.000000 ;// (mean distance, a.u.) $e = 0.016709 - .000000001151 * $d ;// (eccentricity) $M = 6.2141924418482506287494932704807 + 0.017201969619332228715501852561964 * $d ;// (mean anomaly) $oblecl = 0.40909295936270689252387465029835 - .0000000062186081248557962825791102081249 * $d ;// obliquity of the ecliptic $L = $w + $M; // sun's mean longitude $E = $M + $e * sin($M) * (1 + $e * cos($M)); $x = cos($E) - $e; $y = sin($E) * sqrt(1 - $e * $e); $r = sqrt($x*$x + $y*$y); $v = atan2( $y, $x ) ; $lon = $v + $w; $x = $r * cos($lon); $y = $r * sin($lon); $z = 0.0; $xequat = $x; $yequat = $y * cos($oblecl) + $z * sin($oblecl); $zequat = $y * sin($oblecl) + $z * cos($oblecl); $RA = atan2( $yequat, $xequat ); $Decl = asin( $zequat / $r ); $RAhours = r2d($RA)/15; $GMST0 = r2d( $L + pi() ) / 15;// $SIDTIME = $GMST0 + $hour + rad2deg($LON)/15; $HA = deg2rad(($SIDTIME - $RAhours) *15); $x = cos($HA) * cos($Decl); $y = sin($HA) * cos($Decl); $z = sin($Decl); $xhor = $x * cos(pi()/2 - $LAT) - $z * sin(pi()/2 - $LAT); $yhor = $y; $zhor = $x * sin(pi()/2 - $LAT) + $z * cos(pi()/2 - $LAT); $azimuth = pi()*3/2 -atan2( $yhor, $xhor ); $altitude = asin( $zhor ); // check sunshine state $alt_d = rad2deg($altitude); if ($alt_d >= 0) $sunstate = "DAY"; else if ($alt_d >= -6) $sunstate = "CIVIL-TWILIGHT"; else if ($alt_d >= -12) $sunstate = "NAUTICAL-TWILIGHT"; else if ($alt_d >= -18) $sunstate = "ASTRONOMICAL-TWILIGHT"; else $sunstate = "NIGHT"; } function r2d($r) { $d = rad2deg($r); while ($d<0) $d += 360; while ($d>360) $d -= 360; return $d; } // functions below are for the sun's path through the day function findCentre($x1,$y1,$x2,$y2,$x3,$y3,&$cx,&$cy,&$r) { if (!isPerpendicular($x1,$y1,$x2,$y2,$x3,$y3) ) CalcCircle($x1,$y1,$x2,$y2,$x3,$y3,$cx,$cy,$r); else if (!isPerpendicular($x1,$y1,$x3,$y3,$x2,$y2) ) CalcCircle($x1,$y1,$x3,$y3,$x2,$y2,$cx,$cy,$r); else if (!isPerpendicular($x2,$y2,$x1,$y1,$x3,$y3) ) CalcCircle($x2,$y2,$x1,$y1,$x3,$y3,$cx,$cy,$r); else if (!isPerpendicular($x2,$y2,$x3,$y3,$x1,$y1) ) CalcCircle($x2,$y2,$x3,$y3,$x1,$y1,$cx,$cy,$r); else if (!isPerpendicular($x3,$y3,$x2,$y2,$x1,$y1) ) CalcCircle($x3,$y3,$x2,$y2,$x1,$y1,$cx,$cy,$r); else if (!isPerpendicular($x3,$y3,$x1,$y1,$x2,$y2) ) CalcCircle($x3,$y3,$x1,$y1,$x2,$y2,$cx,$cy,$r); else { $cx=0; $cy=0; $r=0; } } function isPerpendicular($x1,$y1,$x2,$y2,$x3,$y3) { $m = 0.00000001; $dya = $y2 - $y1; $dxa = $x2 - $x1; $dyb = $y3 - $y2; $dxb = $x3 - $x2; // checking whether the line of the two pts are vertical if (abs($dxa) <= $m && abs($dyb) <= $m){ //TRACE("The points are pependicular and parallel to x-y axis\n"); return false; } if (abs($dxa) <= $m || abs($dxb) <= $m || abs($dya) < $m || abs($dyb) <= $m) { return true; } } function CalcCircle($x1,$y1,$x2,$y2,$x3,$y3,&$cx,&$cy,&$r) { $m = 0.00000001; $dya = $y2 - $y1; $dxa = $x2 - $x1; $dyb = $y3 - $y2; $dxb = $x3 - $x2; if (abs($dxa) <= $m && abs($dxb) <= $m){ $cx= 0.5 * ($x2 + $x3); $cy= 0.5 * ($y1 + $y2); $r = sqrt(($x1-$cx)*($x1-$cx) + ($y1-$cy)*($y1-$cy)); return; } // IsPerpendicular() assure that xDelta(s) are not zero $aSlope = $dya/$dxa; $bSlope = $dyb/$dxb; if (abs($aSlope-$bSlope) <= 0.000000001){ // checking whether the given points are colinear. return; } // calc center $cx = ($aSlope*$bSlope*($y1 - $y3) + $bSlope*($x1 + $x2) - $aSlope*($x2+$x3) )/(2* ($bSlope-$aSlope) ); $cy = -1*($cx - ($x1+$x2)/2)/$aSlope + ($y1+$y2)/2; $r = sqrt(($x1-$cx)*($x1-$cx) + ($y1-$cy)*($y1-$cy)); return; } ?>
Initial URL
Initial Description
Initial Title
php sun position
Initial Tags
php
Initial Language
PHP