Position of the sun given time of day, latitude and longitude

Question

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.

Solution 1

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.

The short answer

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.

Discussion

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).

Accuracy of the function 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

Other errors in the code

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

Corrected version of 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))
}

References:

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.

Solution 2

Using "NOAA Solar Calculations" from one of the links above I have changed a bit the final part of the function by using a slighly different algorithm that, I hope, have translated without errors. I have commented out the now-useless code and added the new algorithm just after the latitude to radians conversion:

# -----------------------------------------------
# New code
# Solar zenith angle
zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
# Solar azimuth
az <- acos(((sin(lat) * cos(zenithAngle)) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
rm(zenithAngle)
# -----------------------------------------------

# 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

# -----------------------------------------------
# New code
if (ha > 0) az <- az + 180 else az <- 540 - az
az <- az %% 360
# -----------------------------------------------

return(list(elevation=el, azimuth=az))

To verify azimuth trend in the four cases you mentioned let's plot it against time of day:

hour <- seq(from = 0, to = 23, by = 0.5)
azimuth <- data.frame(hour = hour)
az41S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-41,0)$azimuth)
az03S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-03,0)$azimuth)
az03N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,03,0)$azimuth)
az41N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,41,0)$azimuth)
azimuth <- cbind(azimuth, az41S, az03S, az41N, az03N)
rm(az41S, az03S, az41N, az03N)
library(ggplot2)
azimuth.plot <- melt(data = azimuth, id.vars = "hour")
ggplot(aes(x = hour, y = value, color = variable), data = azimuth.plot) + 
    geom_line(size = 2) + 
    geom_vline(xintercept = 12) + 
    facet_wrap(~ variable)

Image attached:

Solution 3

Here's a rewrite in that's more idiomatic to R, and easier to debug and maintain. It is essentially Josh's answer, but with azimuth calculated using both Josh and Charlie's algorithms for comparison. I've also included the simplifications to the date code from my other answer. The basic principle was to split the code up into lots of smaller functions that you can more easily write unit tests for.

astronomersAlmanacTime <- function(x)
{
  # Astronomer's almanach time is the number of 
  # days since (noon, 1 January 2000)
  origin <- as.POSIXct("2000-01-01 12:00:00")
  as.numeric(difftime(x, origin, units = "days"))
}

hourOfDay <- function(x)
{
  x <- as.POSIXlt(x)
  with(x, hour + min / 60 + sec / 3600)
}

degreesToRadians <- function(degrees)
{
  degrees * pi / 180
}

radiansToDegrees <- function(radians)
{
  radians * 180 / pi
}

meanLongitudeDegrees <- function(time)
{
  (280.460 + 0.9856474 * time) %% 360
}

meanAnomalyRadians <- function(time)
{
  degreesToRadians((357.528 + 0.9856003 * time) %% 360)
}

eclipticLongitudeRadians <- function(mnlong, mnanom)
{
  degreesToRadians(
      (mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)) %% 360
  )
}

eclipticObliquityRadians <- function(time)
{
  degreesToRadians(23.439 - 0.0000004 * time)
}

rightAscensionRadians <- function(oblqec, eclong)
{
  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] + 2 * pi 
  ra
}

rightDeclinationRadians <- function(oblqec, eclong)
{
  asin(sin(oblqec) * sin(eclong))
}

greenwichMeanSiderealTimeHours <- function(time, hour)
{
  (6.697375 + 0.0657098242 * time + hour) %% 24
}

localMeanSiderealTimeRadians <- function(gmst, long)
{
  degreesToRadians(15 * ((gmst + long / 15) %% 24))
}

hourAngleRadians <- function(lmst, ra)
{
  ((lmst - ra + pi) %% (2 * pi)) - pi
}

elevationRadians <- function(lat, dec, ha)
{
  asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
}

solarAzimuthRadiansJosh <- function(lat, dec, ha, el)
{
  az <- asin(-cos(dec) * sin(ha) / cos(el))
  cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
  sinAzNeg <- (sin(az) < 0)
  az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + 2 * pi
  az[!cosAzPos] <- pi - az[!cosAzPos]
  az
}

solarAzimuthRadiansCharlie <- function(lat, dec, ha)
{
  zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
  az <- acos((sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
  ifelse(ha > 0, az + pi, 3 * pi - az) %% (2 * pi)
}

sunPosition <- function(when = Sys.time(), format, lat = 46.5, long = 6.5) 
{    
  if(is.character(when)) when <- strptime(when, format)
  when <- lubridate::with_tz(when, "UTC")
  time <- astronomersAlmanacTime(when)
  hour <- hourOfDay(when)

  # Ecliptic coordinates  
  mnlong <- meanLongitudeDegrees(time)   
  mnanom <- meanAnomalyRadians(time)  
  eclong <- eclipticLongitudeRadians(mnlong, mnanom)     
  oblqec <- eclipticObliquityRadians(time)

  # Celestial coordinates
  ra <- rightAscensionRadians(oblqec, eclong)
  dec <- rightDeclinationRadians(oblqec, eclong)

  # Local coordinates
  gmst <- greenwichMeanSiderealTimeHours(time, hour)  
  lmst <- localMeanSiderealTimeRadians(gmst, long)

  # Hour angle
  ha <- hourAngleRadians(lmst, ra)

  # Latitude to radians
  lat <- degreesToRadians(lat)

  # Azimuth and elevation
  el <- elevationRadians(lat, dec, ha)
  azJ <- solarAzimuthRadiansJosh(lat, dec, ha, el)
  azC <- solarAzimuthRadiansCharlie(lat, dec, ha)

  data.frame(
      elevation = radiansToDegrees(el), 
      azimuthJ  = radiansToDegrees(azJ),
      azimuthC  = radiansToDegrees(azC)
  )
}

Solution 4

This is a suggested update to Josh's excellent answer.

Much of the start of the function is boilerplate code for calculating the number of days since midday on 1st Jan 2000. This is much better dealt with using R's existing date and time function.

I also think that rather than having six different variables to specify the date and time, it's easier (and more consistent with other R functions) to specify an existing date object or a date strings + format strings.

Here are two helper functions

astronomers_almanac_time <- function(x)
{
  origin <- as.POSIXct("2000-01-01 12:00:00")
  as.numeric(difftime(x, origin, units = "days"))
}

hour_of_day <- function(x)
{
  x <- as.POSIXlt(x)
  with(x, hour + min / 60 + sec / 3600)
}

And the start of the function now simplifies to

sunPosition <- function(when = Sys.time(), format, lat=46.5, long=6.5) {

  twopi <- 2 * pi
  deg2rad <- pi / 180

  if(is.character(when)) when <- strptime(when, format)
  time <- astronomers_almanac_time(when)
  hour <- hour_of_day(when)
  #...

---

The other oddity is in the lines like

mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

Since mnlong has had %% called on its values, they should all be non-negative already, so this line is superfluous.

Solution 5

I needed sun position in a Python project. I adapted Josh O'Brien's algorithm.

Thank you Josh.

In case it could be useful to anyone, here's my adaptation.

Note that my project only needed instant sun position so time is not a parameter.

def sunPosition(lat=46.5, long=6.5):

    # Latitude [rad]
    lat_rad = math.radians(lat)

    # Get Julian date - 2400000
    day = time.gmtime().tm_yday
    hour = time.gmtime().tm_hour + \
           time.gmtime().tm_min/60.0 + \
           time.gmtime().tm_sec/3600.0
    delta = time.gmtime().tm_year - 1949
    leap = delta / 4
    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)
    t = jd - 51545

    # Ecliptic coordinates

    # Mean longitude
    mnlong_deg = (280.460 + .9856474 * t) % 360

    # Mean anomaly
    mnanom_rad = math.radians((357.528 + .9856003 * t) % 360)

    # Ecliptic longitude and obliquity of ecliptic
    eclong = math.radians((mnlong_deg + 
                           1.915 * math.sin(mnanom_rad) + 
                           0.020 * math.sin(2 * mnanom_rad)
                          ) % 360)
    oblqec_rad = math.radians(23.439 - 0.0000004 * t)

    # Celestial coordinates
    # Right ascension and declination
    num = math.cos(oblqec_rad) * math.sin(eclong)
    den = math.cos(eclong)
    ra_rad = math.atan(num / den)
    if den < 0:
        ra_rad = ra_rad + math.pi
    elif num < 0:
        ra_rad = ra_rad + 2 * math.pi
    dec_rad = math.asin(math.sin(oblqec_rad) * math.sin(eclong))

    # Local coordinates
    # Greenwich mean sidereal time
    gmst = (6.697375 + .0657098242 * t + hour) % 24
    # Local mean sidereal time
    lmst = (gmst + long / 15) % 24
    lmst_rad = math.radians(15 * lmst)

    # Hour angle (rad)
    ha_rad = (lmst_rad - ra_rad) % (2 * math.pi)

    # Elevation
    el_rad = math.asin(
        math.sin(dec_rad) * math.sin(lat_rad) + \
        math.cos(dec_rad) * math.cos(lat_rad) * math.cos(ha_rad))

    # Azimuth
    az_rad = math.asin(
        - math.cos(dec_rad) * math.sin(ha_rad) / math.cos(el_rad))

    if (math.sin(dec_rad) - math.sin(el_rad) * math.sin(lat_rad) < 0):
        az_rad = math.pi - az_rad
    elif (math.sin(az_rad) < 0):
        az_rad += 2 * math.pi

    return el_rad, az_rad

View more solutions

Comments

[mbask]

@Josh O'Brien: Your very detailed answer is an excellent read. As a related note, our SunPosition functions yield exactly the same results.

[mdsumner]

I attached the image file, if you want it.

[Josh O'Brien]

@Charlie - Great answer, and the plots are an especially nice addition. Before seeing them, I hadn't appreciated how different the night-time azimuthal coordinates of the sun would be at 'equatorial' vs more 'temperate' locations. Truly cool.

[SpoonNZ]

Thanks for the fantastic answer! I haven't been on here over the weekend so missed it sorry. Won't get a chance to try this until at least tonight, but looks like it'll do the trick. Cheers!

[Josh O'Brien]

@SpoonNZ -- My pleasure. If you need pdf copies of any of those cited references, let me know at my email address, and I can send them to you.

[Richie Cotton]

@JoshO'Brien: Just added some suggestions in a separate answer. You might want to take a look and incorporate them in your own.

[SpoonNZ]

Great thanks! As mentioned, I've ported this to PHP (and will go to Javascript probably - just gotta decide where I want what functions handled) so that code isn't much help to me, but should be able to be ported (though with slightly more thinking involved than with the original code!). I need to tweak the code that handles the time zones a little, so might be able to integrate this change at the same time.

[mbask]

Nifty changes @Richie Cotton. Note that the assignment hour <- hour_of_day should actually be hour <- hour_of_day(when) and that variable time should hold the number of days, not an object of class "difftime". The second line of function astronomers_almanac_time should be changed to something like as.numeric(difftime(x, origin, units = "days"), units = "days").

[Josh O'Brien]

@RichieCotton -- Thanks for posting your suggestions. I'm not going to add them here, but only because they are R-specific and the OP was using the R-code to try to debug it before porting it to another language. (In fact, I've just now edited my post to correct a date-processing error in the original code, and it's exactly the sort of error that argues for using higher-level code like that you proposed.) Cheers!

[Josh O'Brien]

Thanks for the great suggestions. It might be nice (if you're interested) to include in your post an edited version of the entire sunPosition() function that is more R-ish in its construction.

[Richie Cotton]

@JoshO'Brien: Done. I've made the answer community wiki, since it's a combination of all our answers. It give the same answer as yours for the current time and default (Swiss?) coordinates but lots more testing is needed.

[Josh O'Brien]

@RichieCotton -- What a nice idea. I'll take a deeper look at what you've done, as soon as I get the chance.

[Neon22]

Note when testing against NOAA website here: esrl.noaa.gov/gmd/grad/solcalc/azel.html That NOAA uses Longitude West as +ve. This algorithm uses Longitude West as -ve.

[Mark Ireland]

This was really useful to me. Thanks. One thing I did do is add an adjustment for Daylight Savings. In case it's of use, it was simply: if (time.localtime().tm_isdst == 1): hour += 1

[Hawk]

One could also combine the julian dates to: time = 365*(year - 2000) + floor((year - 1949)/4) + day + hour - 13.5

[Tedward]

When I run "sunPosition(lat = 43, long = -89)", I get an elevation of 52 and an azimuth of 175. But using NOAA's web app esrl.noaa.gov/gmd/grad/solcalc, I get elevation of around 5 and azimuth of 272. Am I missing something? NOAA is correct, but I can't get sunPosition to give accurate results.

[Richie Cotton]

@Tedward sunPosition defaults to using the current time and date. Is that what you wanted?

[Tedward]

Yes. I also tested with some different times. This was late in the day, I'm going to try again today with a fresh start. I'm pretty certain I'm doing something wrong, but don't know what. I'll keep working at it.

[Tedward]

I needed to convert "when" to UTC to get accurate results. See stackoverflow.com/questions/39393514/…. @aichao suggests code for conversion.

[Richie Cotton]

@Tedward I added a call to lubridate::with_tz in sunPosition. Please can you check to see if this gives the correct answer now.