1# 2# Calculate solar position in equatorial coordinates (ascension, declination) 3# then convert to horizontal coordinates (altitude, azimuth). 4# Parameters L, g, lambda, and eps are appoximated using values chosen 5# for accuracy in the era on either side of Jan 2000. 6# 7# Input: 8# Date - string with format "dd-mm-YYYY HH:MM" 9# Latitude - in degrees 10# Longitude - in degrees 11# 12# Output: 13# Function definitions 14# Azimuth(t) = Solar azimuth at time t 15# Altitude(t) = Solar altitude at time t 16# where t is time in seconds relative to local solar noon on Date. 17# sunrise, sunset (string containing local standard time) 18# sunlight (string containing time between sunrise, sunset) 19# 20# Local values: 21 Minute = 60. 22 Hour = 3600. 23 Day = 86400. 24 TimeFormat = "%d-%m-%Y %H:%M" 25 Phi = Latitude 26 set angle degrees 27# 28# n = local time in days since noon on 1 Jan 2000 29# L = mean solar longitude in ecliptic coordinates 30# = 280.460° + 0.9856474° * n 31# g = orbital anomaly 32# = 357.528° + 0.9856003° * n 33# lamba = solar longitude 34# = L + 1.915° * sin(g) + 0.020° * sin(2g) 35# eps = obliquity of the ecliptic 36# = 23.439° - 0.0000004° * n 37 38n = (strptime(TimeFormat, Date) - strptime(TimeFormat, "01-01-2000 12:00")) / Day 39 40L = n * 0.9856474 + 280.460 41L = L - 360 * int(L / 360.) 42if (L < 0) { L += 360. } 43 44g = n * 0.9856003 + 357.528 45g = g - 360 * int(g / 360.) 46if (g < 0) { g += 360. } 47 48lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g) 49 50eps = 23.439 - n * 0.0000004 51 52# Equatorial coordinates 53# RAsc = right ascension 54# Dec = declination 55 56RAsc = atan2( cos(eps) * sin(lambda), cos(lambda) ) 57Dec = asin( sin(eps) * sin(lambda) ) 58 59# Length of day 60daylength = 2.0 * acos( -sin(Dec)*sin(Phi) / cos(Phi) ) 61# Correction for displacement of this location from center of timezone 62corr = Longitude - floor(Longitude/15.0) * 15.0 - 7.5 63rise = (daylength/2.0 - corr) * Day/360. 64set = (daylength/2.0 + corr) * Day/360. 65 66sunlight = strftime("%tH h %tM m", daylength * Day/360.) 67sunrise = strftime("%tH:%tM", (Day/2)-rise) 68sunset = strftime("%tH:%tM", (Day/2)+set) 69 70# Horizontal coordinates 71 72h(t) = 360. * t / Day 73Altitude(t) = asin( sin(Dec) * sin(Phi) + cos(Dec) * cos(Phi) * cos(h(t)) ) 74 75cosAzi(t) = ( sin(Dec) * sin(Phi) - cos(Dec) * cos(h(t)) * sin(Phi) ) \ 76 / cos(Altitude(t)) 77sinAzi(t) = ( -cos(Dec) * sin(h(t)) ) \ 78 / cos(Altitude(t)) 79Azimuth(t) = atan2( sinAzi(t), cosAzi(t) ) 80 81