# # Calculate solar position in equatorial coordinates (ascension, declination) # then convert to horizontal coordinates (altitude, azimuth). # Parameters L, g, lambda, and eps are appoximated using values chosen # for accuracy in the era on either side of Jan 2000. # # Input: # Date - string with format "dd-mm-YYYY HH:MM" # Latitude - in degrees # Longitude - in degrees # # Output: # Function definitions # Azimuth(t) = Solar azimuth at time t # Altitude(t) = Solar altitude at time t # where t is time in seconds relative to local solar noon on Date. # sunrise, sunset (string containing local standard time) # sunlight (string containing time between sunrise, sunset) # # Local values: Minute = 60. Hour = 3600. Day = 86400. TimeFormat = "%d-%m-%Y %H:%M" Phi = Latitude set angle degrees # # n = local time in days since noon on 1 Jan 2000 # L = mean solar longitude in ecliptic coordinates # = 280.460° + 0.9856474° * n # g = orbital anomaly # = 357.528° + 0.9856003° * n # lamba = solar longitude # = L + 1.915° * sin(g) + 0.020° * sin(2g) # eps = obliquity of the ecliptic # = 23.439° - 0.0000004° * n n = (strptime(TimeFormat, Date) - strptime(TimeFormat, "01-01-2000 12:00")) / Day L = n * 0.9856474 + 280.460 L = L - 360 * int(L / 360.) if (L < 0) { L += 360. } g = n * 0.9856003 + 357.528 g = g - 360 * int(g / 360.) if (g < 0) { g += 360. } lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g) eps = 23.439 - n * 0.0000004 # Equatorial coordinates # RAsc = right ascension # Dec = declination RAsc = atan2( cos(eps) * sin(lambda), cos(lambda) ) Dec = asin( sin(eps) * sin(lambda) ) # Length of day daylength = 2.0 * acos( -sin(Dec)*sin(Phi) / cos(Phi) ) # Correction for displacement of this location from center of timezone corr = Longitude - floor(Longitude/15.0) * 15.0 - 7.5 rise = (daylength/2.0 - corr) * Day/360. set = (daylength/2.0 + corr) * Day/360. sunlight = strftime("%tH h %tM m", daylength * Day/360.) sunrise = strftime("%tH:%tM", (Day/2)-rise) sunset = strftime("%tH:%tM", (Day/2)+set) # Horizontal coordinates h(t) = 360. * t / Day Altitude(t) = asin( sin(Dec) * sin(Phi) + cos(Dec) * cos(Phi) * cos(h(t)) ) cosAzi(t) = ( sin(Dec) * sin(Phi) - cos(Dec) * cos(h(t)) * sin(Phi) ) \ / cos(Altitude(t)) sinAzi(t) = ( -cos(Dec) * sin(h(t)) ) \ / cos(Altitude(t)) Azimuth(t) = atan2( sinAzi(t), cosAzi(t) )