This question has been asked before a little over three years ago. There was an answer given, however I've found a glitch in the solution.
Code below is in R. I've ported it to another language, however have tested the original code directly in R to ensure the issue wasn't with my porting.
sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
lat=46.5, long=6.5) {
twopi <- 2 * pi
deg2rad <- pi / 180
# Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
day <- day + cumsum(month.days)[month]
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
day[leapdays] <- day[leapdays] + 1
# Get Julian date - 2400000
hour <- hour + min / 60 + sec / 3600 # hour plus fraction
delta <- year - 1949
leap <- trunc(delta / 4) # former leapyears
jd <- 32916.5 + delta * 365 + leap + day + hour / 24
# The input to the Atronomer's almanach is the difference between
# the Julian date and JD 2451545.0 (noon, 1 January 2000)
time <- jd - 51545.
# Ecliptic coordinates
# Mean longitude
mnlong <- 280.460 + .9856474 * time
mnlong <- mnlong %% 360
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
# Mean anomaly
mnanom <- 357.528 + .9856003 * time
mnanom <- mnanom %% 360
mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
mnanom <- mnanom * deg2rad
# Ecliptic longitude and obliquity of ecliptic
eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
eclong <- eclong %% 360
eclong[eclong < 0] <- eclong[eclong < 0] + 360
oblqec <- 23.429 - 0.0000004 * time
eclong <- eclong * deg2rad
oblqec <- oblqec * deg2rad
# Celestial coordinates
# Right ascension and declination
num <- cos(oblqec) * sin(eclong)
den <- cos(eclong)
ra <- atan(num / den)
ra[den < 0] <- ra[den < 0] + pi
ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
dec <- asin(sin(oblqec) * sin(eclong))
# Local coordinates
# Greenwich mean sidereal time
gmst <- 6.697375 + .0657098242 * time + hour
gmst <- gmst %% 24
gmst[gmst < 0] <- gmst[gmst < 0] + 24.
# Local mean sidereal time
lmst <- gmst + long / 15.
lmst <- lmst %% 24.
lmst[lmst < 0] <- lmst[lmst < 0] + 24.
lmst <- lmst * 15. * deg2rad
# Hour angle
ha <- lmst - ra
ha[ha < -pi] <- ha[ha < -pi] + twopi
ha[ha > pi] <- ha[ha > pi] - twopi
# Latitude to radians
lat <- lat * deg2rad
# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
az <- asin(-cos(dec) * sin(ha) / cos(el))
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad
return(list(elevation=el, azimuth=az))
}
The problem I'm hitting is that the azimuth it returns seems wrong. For example, if I run the function on the (southern) summer solstice at 12:00 for locations 0ºE and 41ºS, 3ºS, 3ºN and 41ºN:
> sunPosition(2012,12,22,12,0,0,-41,0)
$elevation
[1] 72.42113
$azimuth
[1] 180.9211
> sunPosition(2012,12,22,12,0,0,-3,0)
$elevation
[1] 69.57493
$azimuth
[1] -0.79713
Warning message:
In asin(sin(dec)/sin(lat)) : NaNs produced
> sunPosition(2012,12,22,12,0,0,3,0)
$elevation
[1] 63.57538
$azimuth
[1] -0.6250971
Warning message:
In asin(sin(dec)/sin(lat)) : NaNs produced
> sunPosition(2012,12,22,12,0,0,41,0)
$elevation
[1] 25.57642
$azimuth
[1] 180.3084
These numbers just don't seem right. The elevation I'm happy with - the first two should be roughly the same, the third a touch lower, and the fourth much lower. However the first azimuth should be roughly due North, whereas the number it gives is the complete opposite. The remaining three should point roughly due South, however only the last one does. The two in the middle point just off North, again 180º out.
As you can see there are also a couple of errors triggered with the low latitudes (close the equator)
I believe the fault is in this section, with the error being triggered at the third line (starting with elc).
# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
az <- asin(-cos(dec) * sin(ha) / cos(el))
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
I googled around and found a similar chunk of code in C, converted to R the line it uses to calculate the azimuth would be something like
az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
The output here seems to be heading in the right direction, but I just can't get it to give me the right answer all the time when it's converted back to degrees.
A correction of the code (suspect it's just the few lines above) to make it calculate the correct azimuth would be fantastic.
This seems like an important topic, so I've posted a longer than typical answer: if this algorithm is to be used by others in the future, I think it's important that it be accompanied by references to the literature from which it has been derived.
As you've noted, your posted code does not work properly for locations near the equator, or in the southern hemisphere.
To fix it, simply replace these lines in your original code:
elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi
with these:
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]
It should now work for any location on the globe.
The code in your example is adapted almost verbatim from a 1988 article by J.J. Michalsky (Solar Energy. 40:227-235). That article in turn refined an algorithm presented in a 1978 article by R. Walraven (Solar Energy. 20:393-397). Walraven reported that the method had been used successfully for several years to precisely position a polarizing radiometer in Davis, CA (38° 33' 14" N, 121° 44' 17" W).
Both Michalsky's and Walraven's code contains important/fatal errors. In particular, while Michalsky's algorithm works just fine in most of the United States, it fails (as you've found) for areas near the equator, or in the southern hemisphere. In 1989, J.W. Spencer of Victoria, Australia, noted the same thing (Solar Energy. 42(4):353):
Dear Sir:
Michalsky's method for assigning the calculated azimuth to the correct quadrant, derived from Walraven, does not give correct values when applied for Southern (negative) latitudes. Further the calculation of the critical elevation (elc) will fail for a latitude of zero because of division by zero. Both these objections can be avoided simply by assigning the azimuth to the correct quadrant by considering the sign of cos(azimuth).
My edits to your code are based on the corrections suggested by Spencer in that published Comment. I have simply altered them somewhat to ensure that the R function sunPosition() remains 'vectorized' (i.e. working properly on vectors of point locations, rather than needing to be passed one point at a time).
sunPosition()To test that sunPosition() works correctly, I've compared its results with those calculated by the National Oceanic and Atmospheric Administration's Solar Calculator. In both cases, sun positions were calculated for midday (12:00 PM) on the southern summer solstice (December 22nd), 2012. All results were in agreement to within 0.02 degrees.
testPts <- data.frame(lat = c(-41,-3,3, 41),
long = c(0, 0, 0, 0))
# Sun's position as returned by the NOAA Solar Calculator,
NOAA <- data.frame(elevNOAA = c(72.44, 69.57, 63.57, 25.6),
azNOAA = c(359.09, 180.79, 180.62, 180.3))
# Sun's position as returned by sunPosition()
sunPos <- sunPosition(year = 2012,
month = 12,
day = 22,
hour = 12,
min = 0,
sec = 0,
lat = testPts$lat,
long = testPts$long)
cbind(testPts, NOAA, sunPos)
# lat long elevNOAA azNOAA elevation azimuth
# 1 -41 0 72.44 359.09 72.43112 359.0787
# 2 -3 0 69.57 180.79 69.56493 180.7965
# 3 3 0 63.57 180.62 63.56539 180.6247
# 4 41 0 25.60 180.30 25.56642 180.3083
There are at least two other (quite minor) errors in the posted code. The first causes February 29th and March 1st of leap years to both be tallied as day 61 of the year. The second error derives from a typo in the original article, which was corrected by Michalsky in a 1989 note (Solar Energy. 43(5):323).
This code block shows the offending lines, commented out and followed immediately by corrected versions:
# leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
day >= 60 & !(month==2 & day==60)
# oblqec <- 23.429 - 0.0000004 * time
oblqec <- 23.439 - 0.0000004 * time
sunPosition()Here is the corrected code that was verified above:
sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
lat=46.5, long=6.5) {
twopi <- 2 * pi
deg2rad <- pi / 180
# Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
day <- day + cumsum(month.days)[month]
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
day >= 60 & !(month==2 & day==60)
day[leapdays] <- day[leapdays] + 1
# Get Julian date - 2400000
hour <- hour + min / 60 + sec / 3600 # hour plus fraction
delta <- year - 1949
leap <- trunc(delta / 4) # former leapyears
jd <- 32916.5 + delta * 365 + leap + day + hour / 24
# The input to the Atronomer's almanach is the difference between
# the Julian date and JD 2451545.0 (noon, 1 January 2000)
time <- jd - 51545.
# Ecliptic coordinates
# Mean longitude
mnlong <- 280.460 + .9856474 * time
mnlong <- mnlong %% 360
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
# Mean anomaly
mnanom <- 357.528 + .9856003 * time
mnanom <- mnanom %% 360
mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
mnanom <- mnanom * deg2rad
# Ecliptic longitude and obliquity of ecliptic
eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
eclong <- eclong %% 360
eclong[eclong < 0] <- eclong[eclong < 0] + 360
oblqec <- 23.439 - 0.0000004 * time
eclong <- eclong * deg2rad
oblqec <- oblqec * deg2rad
# Celestial coordinates
# Right ascension and declination
num <- cos(oblqec) * sin(eclong)
den <- cos(eclong)
ra <- atan(num / den)
ra[den < 0] <- ra[den < 0] + pi
ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
dec <- asin(sin(oblqec) * sin(eclong))
# Local coordinates
# Greenwich mean sidereal time
gmst <- 6.697375 + .0657098242 * time + hour
gmst <- gmst %% 24
gmst[gmst < 0] <- gmst[gmst < 0] + 24.
# Local mean sidereal time
lmst <- gmst + long / 15.
lmst <- lmst %% 24.
lmst[lmst < 0] <- lmst[lmst < 0] + 24.
lmst <- lmst * 15. * deg2rad
# Hour angle
ha <- lmst - ra
ha[ha < -pi] <- ha[ha < -pi] + twopi
ha[ha > pi] <- ha[ha > pi] - twopi
# Latitude to radians
lat <- lat * deg2rad
# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
az <- asin(-cos(dec) * sin(ha) / cos(el))
# For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]
# if (0 < sin(dec) - sin(el) * sin(lat)) {
# if(sin(az) < 0) az <- az + twopi
# } else {
# az <- pi - az
# }
el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad
return(list(elevation=el, azimuth=az))
}
Michalsky, J.J. 1988. The Astronomical Almanac's algorithm for approximate solar position (1950-2050). Solar Energy. 40(3):227-235.
Michalsky, J.J. 1989. Errata. Solar Energy. 43(5):323.
Spencer, J.W. 1989. Comments on "The Astronomical Almanac's Algorithm for Approximate Solar Position (1950-2050)". Solar Energy. 42(4):353.
Walraven, R. 1978. Calculating the position of the sun. Solar Energy. 20:393-397.