## Posted By

tomas on 10/14/10

# php sun position

/ Published in: PHP

`<?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; 	} ?>`