fork download
copy 
  1. # Code from http://stackoverflow.com/a/8764866/1468366 by Josh O'Brien,
  2. # augmented for for http://stackoverflow.com/q/19165498/1468366
  3.  
  4. sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
  5.                     lat=46.5, long=6.5) {
  6.  
  7.     twopi <- 2 * pi
  8.     cat("twopi =", twopi, "\n")
  9.     deg2rad <- pi / 180
  10.     cat("deg2rad =", deg2rad, "\n")
  11.  
  12.     # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
  13.     month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
  14.     cat("month.days =", month.days, "\n")
  15.     day <- day + cumsum(month.days)[month]
  16.     cat("day =", day, "\n")
  17.     leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & 
  18.                 day >= 60 & !(month==2 & day==60)
  19.     cat("leapdays =", leapdays, "\n")
  20.     day[leapdays] <- day[leapdays] + 1
  21.     cat("day[", leapdays, "] =", day[leapdays], "\n")
  22.     cat("day =", day, "\n")
  23.  
  24.     # Get Julian date - 2400000
  25.     hour <- hour + min / 60 + sec / 3600 # hour plus fraction
  26.     cat("hour =", hour, "\n")
  27.     delta <- year - 1949
  28.     cat("delta =", delta, "\n")
  29.     leap <- trunc(delta / 4) # former leapyears
  30.     cat("leap =", leap, "\n")
  31.     jd <- 32916.5 + delta * 365 + leap + day + hour / 24
  32.     cat("jd =", jd, "\n")
  33.  
  34.     # The input to the Atronomer's almanach is the difference between
  35.     # the Julian date and JD 2451545.0 (noon, 1 January 2000)
  36.     time <- jd - 51545.
  37.     cat("time =", time, "\n")
  38.  
  39.     # Ecliptic coordinates
  40.  
  41.     # Mean longitude
  42.     mnlong <- 280.460 + .9856474 * time
  43.     cat("mnlong =", mnlong, "\n")
  44.     mnlong <- mnlong %% 360
  45.     cat("mnlong =", mnlong, "\n")
  46.     cat("mnlong[", mnlong < 0, "] =", mnlong[mnlong < 0] + 360, "\n")
  47.     mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
  48.     cat("mnlong =", mnlong, "\n")
  49.  
  50.     # Mean anomaly
  51.     mnanom <- 357.528 + .9856003 * time
  52.     cat("mnanom =", mnanom, "\n")
  53.     mnanom <- mnanom %% 360
  54.     cat("mnanom =", mnanom, "\n")
  55.     cat("mnanom[", mnanom < 0, "] =", mnanom[mnanom < 0] + 360, "\n")
  56.     mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
  57.     cat("mnanom =", mnanom, "\n")
  58.     mnanom <- mnanom * deg2rad
  59.     cat("mnanom =", mnanom, "\n")
  60.  
  61.     # Ecliptic longitude and obliquity of ecliptic
  62.     eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
  63.     cat("eclong =", eclong, "\n")
  64.     eclong <- eclong %% 360
  65.     cat("eclong =", eclong, "\n")
  66.     cat("eclong[", eclong < 0, "] =", eclong[eclong < 0] + 360, "\n")
  67.     eclong[eclong < 0] <- eclong[eclong < 0] + 360
  68.     cat("eclong =", eclong, "\n")
  69.     oblqec <- 23.439 - 0.0000004 * time
  70.     cat("oblqec =", oblqec, "\n")
  71.     eclong <- eclong * deg2rad
  72.     cat("eclong =", eclong, "\n")
  73.     oblqec <- oblqec * deg2rad
  74.     cat("oblqec =", oblqec, "\n")
  75.  
  76.     # Celestial coordinates
  77.     # Right ascension and declination
  78.     num <- cos(oblqec) * sin(eclong)
  79.     cat("num =", num, "\n")
  80.     den <- cos(eclong)
  81.     cat("den =", den, "\n")
  82.     ra <- atan(num / den)
  83.     cat("ra =", ra, "\n")
  84.     ra[den < 0] <- ra[den < 0] + pi
  85.     cat("ra[", den < 0, "] =", ra[den < 0], "\n")
  86.     ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
  87.     cat("ra[", den >= 0 & num < 0, "] =", ra[den >= 0 & num < 0], "\n")
  88.     cat("ra =", ra, "\n")
  89.     dec <- asin(sin(oblqec) * sin(eclong))
  90.     cat("dec =", dec, "\n")
  91.  
  92.     # Local coordinates
  93.     # Greenwich mean sidereal time
  94.     gmst <- 6.697375 + .0657098242 * time + hour
  95.     cat("gmst =", gmst, "\n")
  96.     gmst <- gmst %% 24
  97.     cat("gmst =", gmst, "\n")
  98.     cat("gmst[", gmst < 0, "] =", gmst[gmst < 0] + 24., "\n")
  99.     gmst[gmst < 0] <- gmst[gmst < 0] + 24.
  100.  
  101.     # Local mean sidereal time
  102.     lmst <- gmst + long / 15.
  103.     cat("lmst =", lmst, "\n")
  104.     lmst <- lmst %% 24.
  105.     cat("lmst =", lmst, "\n")
  106.     cat("lmst[", lmst < 0, "] =", lmst[lmst < 0] + 24., "\n")
  107.     lmst[lmst < 0] <- lmst[lmst < 0] + 24.
  108.     lmst <- lmst * 15. * deg2rad
  109.     cat("lmst =", lmst, "\n")
  110.  
  111.     # Hour angle
  112.     ha <- lmst - ra
  113.     cat("ha =", ha, "\n")
  114.     cat("ha[", ha < -pi, "] =", ha[ha < -pi] + twopi, "\n")
  115.     ha[ha < -pi] <- ha[ha < -pi] + twopi
  116.     cat("ha[", ha > pi, "] =", ha[ha > pi] - twopi, "\n")
  117.     ha[ha > pi] <- ha[ha > pi] - twopi
  118.  
  119.     # Latitude to radians
  120.     lat <- lat * deg2rad
  121.     cat("lat =", lat, "\n")
  122.  
  123.     # Azimuth and elevation
  124.     el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
  125.     cat("el =", el, "\n")
  126.     az <- asin(-cos(dec) * sin(ha) / cos(el))
  127.     cat("az =", az, "\n")
  128.  
  129.     # For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
  130.     cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
  131.     cat("cosAzPos =", cosAzPos, "\n")
  132.     sinAzNeg <- (sin(az) < 0)
  133.     cat("sinAzNeg =", sinAzNeg, "\n")
  134.     az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
  135.     cat("az[", cosAzPos & sinAzNeg, "] =", az[cosAzPos & sinAzNeg], "\n")
  136.     az[!cosAzPos] <- pi - az[!cosAzPos]
  137.     cat("az[", !cosAzPos, "] =", az[!cosAzPos], "\n")
  138.  
  139.     el <- el / deg2rad
  140.     cat("el =", el, "\n")
  141.     az <- az / deg2rad
  142.     cat("az =", az, "\n")
  143.     lat <- lat / deg2rad
  144.     cat("lat =", lat, "\n")
  145.  
  146.     return(list(elevation=el, azimuth=az))
  147. }
  148.  
  149. sunPosition(2013, 9, 1, 10, 0, 0, -4.77867, 11.86364)
  150.  
Success #stdin #stdout 0.52s 22824KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
twopi = 6.283185 
deg2rad = 0.01745329 
month.days = 0 31 28 31 30 31 30 31 31 30 31 30 
day = 244 
leapdays = FALSE 
day[ FALSE ] =  
day = 244 
hour = 10 
delta = 64 
leap = 16 
jd = 56536.92 
time = 4991.917 
mnlong = 5200.73 
mnlong = 160.7297 
mnlong[ FALSE ] =  
mnlong = 160.7297 
mnanom = 5277.563 
mnanom = 237.5626 
mnanom[ FALSE ] =  
mnanom = 237.5626 
mnanom = 4.146249 
eclong = 159.1316 
eclong = 159.1316 
eclong[ FALSE ] =  
eclong = 159.1316 
oblqec = 23.43700 
eclong = 2.77737 
oblqec = 0.4090529 
num = 0.326834 
den = -0.934401 
ra = -0.3364781 
ra[ TRUE ] = 2.805115 
ra[ FALSE ] =  
ra = 2.805115 
dec = 0.1421627 
gmst = 344.7153 
gmst = 8.715342 
gmst[ FALSE ] =  
lmst = 9.50625 
lmst = 9.50625 
lmst[ FALSE ] =  
lmst = 2.488731 
ha = -0.3163839 
ha[ FALSE ] =  
ha[ FALSE ] =  
lat = -0.08340353 
el = 1.182897 
az = 0.9514718 
cosAzPos = TRUE 
sinAzNeg = FALSE 
az[ FALSE ] =  
az[ FALSE ] =  
el = 67.77503 
az = 54.51532 
lat = -4.77867 
$elevation
[1] 67.77503

$azimuth
[1] 54.51532
https://ideone.com/xGlopE
language:
R (R 3.5.2)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!