# Untitled

By: a guest on May 11th, 2012  |  syntax: None  |  size: 12.91 KB  |  hits: 24  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
1. Position of the sun given time of day, latitude and longitude
2. sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
3.                     lat=46.5, long=6.5) {
4.
5.
6.   twopi <- 2 * pi
7.   deg2rad <- pi / 180
8.
9.   # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
10.   month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
11.   day <- day + cumsum(month.days)[month]
12.   leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
13.   day[leapdays] <- day[leapdays] + 1
14.
15.   # Get Julian date - 2400000
16.   hour <- hour + min / 60 + sec / 3600 # hour plus fraction
17.   delta <- year - 1949
18.   leap <- trunc(delta / 4) # former leapyears
19.   jd <- 32916.5 + delta * 365 + leap + day + hour / 24
20.
21.   # The input to the Atronomer's almanach is the difference between
22.   # the Julian date and JD 2451545.0 (noon, 1 January 2000)
23.   time <- jd - 51545.
24.
25.   # Ecliptic coordinates
26.
27.   # Mean longitude
28.   mnlong <- 280.460 + .9856474 * time
29.   mnlong <- mnlong %% 360
30.   mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
31.
32.   # Mean anomaly
33.   mnanom <- 357.528 + .9856003 * time
34.   mnanom <- mnanom %% 360
35.   mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
36.   mnanom <- mnanom * deg2rad
37.
38.   # Ecliptic longitude and obliquity of ecliptic
39.   eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
40.   eclong <- eclong %% 360
41.   eclong[eclong < 0] <- eclong[eclong < 0] + 360
42.   oblqec <- 23.429 - 0.0000004 * time
43.   eclong <- eclong * deg2rad
44.   oblqec <- oblqec * deg2rad
45.
46.   # Celestial coordinates
47.   # Right ascension and declination
48.   num <- cos(oblqec) * sin(eclong)
49.   den <- cos(eclong)
50.   ra <- atan(num / den)
51.   ra[den < 0] <- ra[den < 0] + pi
52.   ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
53.   dec <- asin(sin(oblqec) * sin(eclong))
54.
55.   # Local coordinates
56.   # Greenwich mean sidereal time
57.   gmst <- 6.697375 + .0657098242 * time + hour
58.   gmst <- gmst %% 24
59.   gmst[gmst < 0] <- gmst[gmst < 0] + 24.
60.
61.   # Local mean sidereal time
62.   lmst <- gmst + long / 15.
63.   lmst <- lmst %% 24.
64.   lmst[lmst < 0] <- lmst[lmst < 0] + 24.
65.   lmst <- lmst * 15. * deg2rad
66.
67.   # Hour angle
68.   ha <- lmst - ra
69.   ha[ha < -pi] <- ha[ha < -pi] + twopi
70.   ha[ha > pi] <- ha[ha > pi] - twopi
71.
73.   lat <- lat * deg2rad
74.
75.   # Azimuth and elevation
76.   el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
77.   az <- asin(-cos(dec) * sin(ha) / cos(el))
78.   elc <- asin(sin(dec) / sin(lat))
79.   az[el >= elc] <- pi - az[el >= elc]
80.   az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
81.
82.   el <- el / deg2rad
83.   az <- az / deg2rad
84.   lat <- lat / deg2rad
85.
86.   return(list(elevation=el, azimuth=az))
87. }
88.
89. > sunPosition(2012,12,22,12,0,0,-41,0)
90. \$elevation
91. [1] 72.42113
92.
93. \$azimuth
94. [1] 180.9211
95.
96. > sunPosition(2012,12,22,12,0,0,-3,0)
97. \$elevation
98. [1] 69.57493
99.
100. \$azimuth
101. [1] -0.79713
102.
103. Warning message:
104. In asin(sin(dec)/sin(lat)) : NaNs produced
105. > sunPosition(2012,12,22,12,0,0,3,0)
106. \$elevation
107. [1] 63.57538
108.
109. \$azimuth
110. [1] -0.6250971
111.
112. Warning message:
113. In asin(sin(dec)/sin(lat)) : NaNs produced
114. > sunPosition(2012,12,22,12,0,0,41,0)
115. \$elevation
116. [1] 25.57642
117.
118. \$azimuth
119. [1] 180.3084
120.
121. # Azimuth and elevation
122.   el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
123.   az <- asin(-cos(dec) * sin(ha) / cos(el))
124.   elc <- asin(sin(dec) / sin(lat))
125.   az[el >= elc] <- pi - az[el >= elc]
126.   az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
127.
128. az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
129.
130. elc <- asin(sin(dec) / sin(lat))
131. az[el >= elc] <- pi - az[el >= elc]
132. az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
133.
134. cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
135. sinAzNeg <- (sin(az) < 0)
136. az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
137. az[!cosAzPos] <- pi - az[!cosAzPos]
138.
139. testPts <- data.frame(lat = c(-41,-3,3, 41),
140.                       long = c(0, 0, 0, 0))
141.
142. # Sun's position as returned by the NOAA Solar Calculator,
143. NOAA <- data.frame(elevNOAA = c(72.44, 69.57, 63.57, 25.6),
144.                    azNOAA = c(359.09, 180.79, 180.62, 180.3))
145.
146. # Sun's position as returned by sunPosition()
147. sunPos <- sunPosition(year = 2012,
148.                       month = 12,
149.                       day = 22,
150.                       hour = 12,
151.                       min = 0,
152.                       sec = 0,
153.                       lat = testPts\$lat,
154.                       long = testPts\$long)
155.
156. cbind(testPts, NOAA, sunPos)
157. #   lat long elevNOAA azNOAA elevation  azimuth
158. # 1 -41    0    72.44 359.09  72.43112 359.0787
159. # 2  -3    0    69.57 180.79  69.56493 180.7965
160. # 3   3    0    63.57 180.62  63.56539 180.6247
161. # 4  41    0    25.60 180.30  25.56642 180.3083
162.
163. # leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
164.   leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
165.               day >= 60 & !(month==2 & day==60)
166.
167. # oblqec <- 23.429 - 0.0000004 * time
168.   oblqec <- 23.439 - 0.0000004 * time
169.
170. sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
171.                     lat=46.5, long=6.5) {
172.
173.     twopi <- 2 * pi
174.     deg2rad <- pi / 180
175.
176.     # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
177.     month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
178.     day <- day + cumsum(month.days)[month]
179.     leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
180.                 day >= 60 & !(month==2 & day==60)
181.     day[leapdays] <- day[leapdays] + 1
182.
183.     # Get Julian date - 2400000
184.     hour <- hour + min / 60 + sec / 3600 # hour plus fraction
185.     delta <- year - 1949
186.     leap <- trunc(delta / 4) # former leapyears
187.     jd <- 32916.5 + delta * 365 + leap + day + hour / 24
188.
189.     # The input to the Atronomer's almanach is the difference between
190.     # the Julian date and JD 2451545.0 (noon, 1 January 2000)
191.     time <- jd - 51545.
192.
193.     # Ecliptic coordinates
194.
195.     # Mean longitude
196.     mnlong <- 280.460 + .9856474 * time
197.     mnlong <- mnlong %% 360
198.     mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
199.
200.     # Mean anomaly
201.     mnanom <- 357.528 + .9856003 * time
202.     mnanom <- mnanom %% 360
203.     mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
204.     mnanom <- mnanom * deg2rad
205.
206.     # Ecliptic longitude and obliquity of ecliptic
207.     eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
208.     eclong <- eclong %% 360
209.     eclong[eclong < 0] <- eclong[eclong < 0] + 360
210.     oblqec <- 23.439 - 0.0000004 * time
211.     eclong <- eclong * deg2rad
212.     oblqec <- oblqec * deg2rad
213.
214.     # Celestial coordinates
215.     # Right ascension and declination
216.     num <- cos(oblqec) * sin(eclong)
217.     den <- cos(eclong)
218.     ra <- atan(num / den)
219.     ra[den < 0] <- ra[den < 0] + pi
220.     ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
221.     dec <- asin(sin(oblqec) * sin(eclong))
222.
223.     # Local coordinates
224.     # Greenwich mean sidereal time
225.     gmst <- 6.697375 + .0657098242 * time + hour
226.     gmst <- gmst %% 24
227.     gmst[gmst < 0] <- gmst[gmst < 0] + 24.
228.
229.     # Local mean sidereal time
230.     lmst <- gmst + long / 15.
231.     lmst <- lmst %% 24.
232.     lmst[lmst < 0] <- lmst[lmst < 0] + 24.
233.     lmst <- lmst * 15. * deg2rad
234.
235.     # Hour angle
236.     ha <- lmst - ra
237.     ha[ha < -pi] <- ha[ha < -pi] + twopi
238.     ha[ha > pi] <- ha[ha > pi] - twopi
239.
241.     lat <- lat * deg2rad
242.
243.     # Azimuth and elevation
244.     el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
245.     az <- asin(-cos(dec) * sin(ha) / cos(el))
246.
247.     # For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
248.     cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
249.     sinAzNeg <- (sin(az) < 0)
250.     az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
251.     az[!cosAzPos] <- pi - az[!cosAzPos]
252.
253.     # if (0 < sin(dec) - sin(el) * sin(lat)) {
254.     #     if(sin(az) < 0) az <- az + twopi
255.     # } else {
256.     #     az <- pi - az
257.     # }
258.
259.
260.     el <- el / deg2rad
261.     az <- az / deg2rad
262.     lat <- lat / deg2rad
263.
264.     return(list(elevation=el, azimuth=az))
265. }
266.
267. # -----------------------------------------------
268. # New code
269. # Solar zenith angle
270. zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
271. # Solar azimuth
272. az <- acos(((sin(lat) * cos(zenithAngle)) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
273. rm(zenithAngle)
274. # -----------------------------------------------
275.
276. # Azimuth and elevation
277. el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
278. #az <- asin(-cos(dec) * sin(ha) / cos(el))
279. #elc <- asin(sin(dec) / sin(lat))
280. #az[el >= elc] <- pi - az[el >= elc]
281. #az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
282.
283. el <- el / deg2rad
284. az <- az / deg2rad
285. lat <- lat / deg2rad
286.
287. # -----------------------------------------------
288. # New code
289. if (ha > 0) az <- az + 180 else az <- 540 - az
290. az <- az %% 360
291. # -----------------------------------------------
292.
293. return(list(elevation=el, azimuth=az))
294.
295. hour <- seq(from = 0, to = 23, by = 0.5)
296. azimuth <- data.frame(hour = hour)
297. az41S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-41,0)\$azimuth)
298. az03S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-03,0)\$azimuth)
299. az03N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,03,0)\$azimuth)
300. az41N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,41,0)\$azimuth)
301. azimuth <- cbind(azimuth, az41S, az03S, az41N, az03N)
302. rm(az41S, az03S, az41N, az03N)
303. library(ggplot2)
304. azimuth.plot <- melt(data = azimuth, id.vars = "hour")
305. ggplot(aes(x = hour, y = value, color = variable), data = azimuth.plot) +
306.     geom_line(size = 2) +
307.     geom_vline(xintercept = 12) +
308.     facet_wrap(~ variable)
309.
310. astronomers_almanac_time <- function(x)
311. {
312.   origin <- as.POSIXct("2000-01-01 12:00:00")
313.   as.numeric(difftime(x, origin, units = "days"))
314. }
315.
316. hour_of_day <- function(x)
317. {
318.   x <- as.POSIXlt(x)
319.   with(x, hour + min / 60 + sec / 3600)
320. }
321.
322. sunPosition <- function(when = Sys.time(), format, lat=46.5, long=6.5) {
323.
324.   twopi <- 2 * pi
325.   deg2rad <- pi / 180
326.
327.   if(is.character(when)) when <- strptime(when, format)
328.   time <- astronomers_almanac_time(when)
329.   hour <- hour_of_day(when)
330.   #...
331.
332. mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
333.
334. astronomersAlmanacTime <- function(x)
335. {
336.   # Astronomer's almanach time is the number of
337.   # days since (noon, 1 January 2000)
338.   origin <- as.POSIXct("2000-01-01 12:00:00")
339.   as.numeric(difftime(x, origin, units = "days"))
340. }
341.
342. hourOfDay <- function(x)
343. {
344.   x <- as.POSIXlt(x)
345.   with(x, hour + min / 60 + sec / 3600)
346. }
347.
349. {
350.   degrees * pi / 180
351. }
352.
354. {
355.   radians * 180 / pi
356. }
357.
358. meanLongitudeDegrees <- function(time)
359. {
360.   (280.460 + 0.9856474 * time) %% 360
361. }
362.
364. {
365.   degreesToRadians((357.528 + 0.9856003 * time) %% 360)
366. }
367.
369. {
371.       (mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)) %% 360
372.   )
373. }
374.
376. {
377.   degreesToRadians(23.439 - 0.0000004 * time)
378. }
379.
381. {
382.   num <- cos(oblqec) * sin(eclong)
383.   den <- cos(eclong)
384.   ra <- atan(num / den)
385.   ra[den < 0] <- ra[den < 0] + pi
386.   ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + 2 * pi
387.   ra
388. }
389.
391. {
392.   asin(sin(oblqec) * sin(eclong))
393. }
394.
395. greenwichMeanSiderealTimeHours <- function(time, hour)
396. {
397.   (6.697375 + 0.0657098242 * time + hour) %% 24
398. }
399.
401. {
402.   degreesToRadians(15 * ((gmst + long / 15) %% 24))
403. }
404.
406. {
407.   ((lmst - ra + pi) %% (2 * pi)) - pi
408. }
409.
410. elevationRadians <- function(lat, dec, ha)
411. {
412.   asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
413. }
414.
415. solarAzimuthRadiansJosh <- function(lat, dec, ha, el)
416. {
417.   az <- asin(-cos(dec) * sin(ha) / cos(el))
418.   cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
419.   sinAzNeg <- (sin(az) < 0)
420.   az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + 2 * pi
421.   az[!cosAzPos] <- pi - az[!cosAzPos]
422.   az
423. }
424.
425. solarAzimuthRadiansCharlie <- function(lat, dec, ha)
426. {
427.   zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
428.   az <- acos((sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
429.   ifelse(ha > 0, az + pi, 3 * pi - az) %% (2 * pi)
430. }
431.
432. sunPosition <- function(when = Sys.time(), format, lat = 46.5, long = 6.5)
433. {
434.   if(is.character(when)) when <- strptime(when, format)
435.   time <- astronomersAlmanacTime(when)
436.   hour <- hourOfDay(when)
437.
438.   # Ecliptic coordinates
439.   mnlong <- meanLongitudeDegrees(time)
443.
444.   # Celestial coordinates
447.
448.   # Local coordinates
449.   gmst <- greenwichMeanSiderealTimeHours(time, hour)
451.
452.   # Hour angle
454.