# Posted By

tomas on 10/14/10

# Statistics

Viewed 464 times
Favorited by 0 user(s)

# php sun position

/ Published in: PHP

Copy this code and paste it in your HTML
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.
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,";
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
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) {
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. ?>