The position of the Sun, taking into account the time of day, latitude and longitude

Question

This question was asked until a little over three years ago. An answer was given, however I found a failure in the solution.

The code below is in R. I ported it to another language, however I tested the source code directly in R to make sure that the problem is not related to 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, eg 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 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 am facing is that the azimuth that it returns seems wrong. For example, if I run the function at the southern summer solstice at 12:00 for places 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 height I'm pleased with is that the first two should be about the same, the third should be touched lower, and the fourth should be much lower. However, the first azimuth must be roughly calculated to the north, while the number it gives is the exact opposite. The remaining three should point roughly south, but only the last. Two in the middle point near the North, again at 180º.

As you can see, there are also a few errors caused by low latitudes (equator closure)

I believe that the error is in this section, while the error runs in 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 piece of code in C converted to R, which it uses to calculate the azimuth, would be something like

 az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat))) 

The result here seems to be going in the right direction, but I just can't get it to give me the correct answer all the time when it is converted back to degrees.

Correcting the code (suppose it's just a few lines above) to calculate the correct azimuth would be fantastic.

Answer 1

This seems to be an important topic, so I published a longer than typical answer: if this algorithm will be used by others in the future, I consider it important that it be accompanied by references to the literature from which it was derived.

Short answer

As you already noted, your published code does not work properly for locations near the equator or in the southern hemisphere.

To fix this, simply replace these lines in the source 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] 

Now he must work anywhere on the globe.

Discussion

The code in your example is adapted almost verbatim from a 1988 article by JJ Michalsky (Solar Energy, 40: 227-235). This article, in turn, refined the algorithm presented in R. Walraven (Solar Energy, 20: 393-397). Valraven reported that this method has been used successfully for several years to accurately position a polarizing radiometer in Davis, California (38 ° 33'14 "N, 121 ° 44'17" W).

Both Michalsky and Walraven codes contain important / fatal errors. In particular, although the Michalsky algorithm works very well in most of the United States, it fails (as you already found) for areas close to the equator or in the southern hemisphere. In 1989, JW Spencer from Victoria, Australia, noted the same (Solar Energy 42 (4): 353):

Dear Sir:

The Mikhalsky method for assigning the calculated azimuth to the correct quadrant obtained from Valraven does not give the correct values ​​when applied to the southern (negative) latitudes. Further, the calculation of the critical mark (elc) will fail for latitude zero due to division by zero. Both of these objections can be avoided simply by assigning the azimuth to the correct quadrant by considering the cos sign (azimuth).

My changes to your code are based on the corrections suggested by Spencer in the posted comment. I just modified them a bit to ensure that the R sunPosition() function remains "vectorized" (i.e., it works correctly on point location vectors, and should not be passed one point at a time).

The accuracy of the sunPosition() function

To verify that sunPosition() working correctly, I compared its results with those calculated by the National Oceanic and Atmospheric Administration Solar Calculator . In both cases, the solar positions were calculated at noon (12:00 PM) for the southern summer solstice (December 22), 2012. All results were consistent with an accuracy of 0.02 degrees.

 testPts <- data.frame(lat = c(-41,-3,3, 41), long = c(0, 0, 0, 0)) # Sun 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 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 code errors

There are at least two other (rather minor) errors in the published code. The first reasons that the February 29 and March 1 leap years coincide coincide with the 61st day. The second error results from a typo in the original article, which was corrected by Mikhalsky in a 1989 note (Solar Energy, 43 (5): 323).

This code block shows the lines of violation, commented out and immediately executed the 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 checked 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, eg 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 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, JW 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)) } 

Literature:

Michalsky, JJ 1988. Astronomical Almanac Algorithm for the Approximate Position of the Sun (1950–2050). Solar energy. 40 (3) :. 227-235

Michalsky, JJ 1989. Correction. Solar energy. 43 (5) :. 323

Spencer, JW 1989. Comments on "The Astronomical Almanac Algorithm for Approximate Solar Position (1950-2050)." Solar energy. 42 (4) :. 353

Walraven, R. 1978. Calculation of the position of the sun. Solar energy. twenty:. 393-397

Answer 2

Using the "NOAA Solar Calculations" from one of the above links, I slightly changed the final part of the function using an algorithm that, I hope, translated without errors. I commented now on the useless code and added a new algorithm immediately after converting the latitude to radians:

 # ----------------------------------------------- # 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 check the azimuthal trend in the four cases you mentioned, let's talk about it with the 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) 

Attached Image:

Answer 3

It corresponds here in what is more idiomatic for R, and easier to debug and maintain. This is essentially Josh, but with an azimuth calculated using Josh and Charlie algorithms for comparison. I also included date code simplifications from my other answer. The basic principle was to split the code into many smaller functions that you can more easily write for unit tests.

 astronomersAlmanacTime <- function(x) { # Astronomer 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) ) } 

Answer 4

This is the recommended update for Josh's excellent answer.

A significant part of the beginning of the function is a template code for calculating the number of days from noon on January 1, 2000. This does a much better job of using the existing date and time function R.

I also believe that instead of specifying the date and time on six different variables, it’s easier (and more consistent with other R functions) to specify an existing date object or a string of date + date 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 beginning of the function is now simplified 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) #... 

---

Another oddity is found in strings like

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

Since mnlong had %% calling its values, they should all be non-negative already, so this line is redundant.

Answer 5

I need the position of the sun in a Python project. I adapted the Josh O'Brien algorithm.

Thanks, Josh.

In case this can be useful to anyone, here is my adaptation.

Please note that my project only requires the instantaneous position of the sun, 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 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 

Answer 6

I am having a little problem with the data point and Richie Cotton functions above (in Charlie code implementation)

 longitude= 176.0433687000000020361767383292317390441894531250 latitude= -39.173830619999996827118593500927090644836425781250 event_time = as.POSIXct("2013-10-24 12:00:00", format="%Y-%m-%d %H:%M:%S", tz = "UTC") sunPosition(when=event_time, lat = latitude, long = longitude) elevation azimuthJ azimuthC 1 -38.92275 180 NaN Warning message: In acos((sin(lat) * cos(zenithAngle) - sin(dec))/(cos(lat) * sin(zenithAngle))) : NaNs produced 

because in the solarAzimuthRadiansCharlie function there was a floating-point expectation around an angle of 180, so (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)) is the smallest number greater than 1, 1.0000000000000004440892098, which generates NaN, since the input for acos should not be higher than 1 or lower - 1.

I suspect there could be similar cases for Josh calculation where floating point rounding effects cause the input for the asin step to be outside of -1: 1, but I did not hit them in my specific dataset.

In half a dozen cases, I hit this “truth” (mid day or night) when the problem arises so that the empirically true value should be 1 / -1. For this reason, it would be convenient for me to fix this by applying a rounding step within the solarAzimuthRadiansJosh and solarAzimuthRadiansCharlie . I'm not sure the theoretical accuracy of the NOAA algorithm (the point at which numerical accuracy ceases to matter anyway), but rounding to 12 decimal places fixed the data in my data set.