Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on May 11th, 2012  |  syntax: None  |  size: 12.91 KB  |  hits: 24  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
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.  
  72.   # Latitude to radians
  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.  
  240.     # Latitude to radians
  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.  
  348. degreesToRadians <- function(degrees)
  349. {
  350.   degrees * pi / 180
  351. }
  352.  
  353. radiansToDegrees <- function(radians)
  354. {
  355.   radians * 180 / pi
  356. }
  357.  
  358. meanLongitudeDegrees <- function(time)
  359. {
  360.   (280.460 + 0.9856474 * time) %% 360
  361. }
  362.  
  363. meanAnomalyRadians <- function(time)
  364. {
  365.   degreesToRadians((357.528 + 0.9856003 * time) %% 360)
  366. }
  367.  
  368. eclipticLongitudeRadians <- function(mnlong, mnanom)
  369. {
  370.   degreesToRadians(
  371.       (mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)) %% 360
  372.   )
  373. }
  374.  
  375. eclipticObliquityRadians <- function(time)
  376. {
  377.   degreesToRadians(23.439 - 0.0000004 * time)
  378. }
  379.  
  380. rightAscensionRadians <- function(oblqec, eclong)
  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.  
  390. rightDeclinationRadians <- function(oblqec, eclong)
  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.  
  400. localMeanSiderealTimeRadians <- function(gmst, long)
  401. {
  402.   degreesToRadians(15 * ((gmst + long / 15) %% 24))
  403. }
  404.  
  405. hourAngleRadians <- function(lmst, ra)
  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)  
  440.   mnanom <- meanAnomalyRadians(time)  
  441.   eclong <- eclipticLongitudeRadians(mnlong, mnanom)    
  442.   oblqec <- eclipticObliquityRadians(time)
  443.  
  444.   # Celestial coordinates
  445.   ra <- rightAscensionRadians(oblqec, eclong)
  446.   dec <- rightDeclinationRadians(oblqec, eclong)
  447.  
  448.   # Local coordinates
  449.   gmst <- greenwichMeanSiderealTimeHours(time, hour)  
  450.   lmst <- localMeanSiderealTimeRadians(gmst, long)
  451.  
  452.   # Hour angle
  453.   ha <- hourAngleRadians(lmst, ra)
  454.  
  455.   # Latitude to radians
  456.   lat <- degreesToRadians(lat)
  457.  
  458.   # Azimuth and elevation
  459.   el <- elevationRadians(lat, dec, ha)
  460.   azJ <- solarAzimuthRadiansJosh(lat, dec, ha, el)
  461.   azC <- solarAzimuthRadiansCharlie(lat, dec, ha)
  462.  
  463.   list(
  464.       elevation = radiansToDegrees(el),
  465.       azimuthJ  = radiansToDegrees(azJ),
  466.       azimuthC  = radiansToDegrees(azC)
  467.   )
  468. }
create a new version of this paste RAW Paste Data