R function for sun position gives unexpected results

Question

I would like to calculate the position of the sun at a given time, latitude and longitude. I found this wonderful question and answer here: The position of the sun, a given time of day, latitude and longitude . However, when I evaluate a function, I get incorrect results. Given the quality of the answer, I’m pretty sure that something is wrong on my part, but I ask this question as a record of an attempt to solve the problem.

Here's the function code reprinted below for convenience:

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

If I run:

 sunPosition(when = Sys.time(),lat = 43, long = -89) 

Result:

  elevation azimuthJ azimuthC 1 -24.56604 55.26111 55.26111 

Sys.time () gives:

 > Sys.time() [1] "2016-09-08 09:09:05 CDT" 

It's 9am and the sun is up. Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ I get an azimuth of 124 and a height of 38, which, in my opinion, is correct.

I thought maybe this is a problem with the code, but I also tried out the original SunPosition function from Josh from the answer above and got the same results. My next thought is that there is a problem with my time or time zone.

---

testing the winter solstice, as done in the above question, still gives the same results that they found and looks right:

 testPts <- data.frame(lat = c(-41,-3,3, 41), long = c(0, 0, 0, 0)) time <- as.POSIXct("2012-12-22 12:00:00") sunPosition(when = time, lat = testPts$lat, long = testPts$long) elevation azimuthJ azimuthC 1 72.43112 359.0787 359.0787 2 69.56493 180.7965 180.7965 3 63.56539 180.6247 180.6247 4 25.56642 180.3083 180.3083 

When I do the same test, but change the longitude (-89), I get a negative increase at noon.

 testPts <- data.frame(lat = c(-41,-3,3, 41), long = c(-89, -89, -89, -89)) time <- as.POSIXct("2012-12-22 12:00:00 CDT") sunPosition(when = time, lat = testPts$lat, long = testPts$long) elevation azimuthJ azimuthC 1 16.060136563 107.3420 107.3420 2 2.387033692 113.3522 113.3522 3 0.001378426 113.4671 113.4671 4 -14.190786786 108.8866 108.8866 

Answer 1

There is nothing wrong with the kernel code found in the linked message if you enter the when input in UTC. The confusion was that the OP entered the wrong Time Zone on the website for Sys.time() of 2016-09-08 09:09:05 CDT :

Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ I get an azimuth of 124 and a height of 38, which in my opinion is correct.

The correct Time Zone to enter the NOAA -5 site for the CDT ( see this site ), which gives

Calling sunPosition with time configured to UTC gives a similar result:

 sunPosition(when = "2016-09-08 14:09:05", format="%Y-%m-%d %H:%M:%S",lat = 43, long = -89) ## elevation azimuthJ azimuthC ##1 28.08683 110.4915 110.4915 

Now the code does not do this conversion to UTC. One way to do this is to replace the first line in sunPosition :

 if(is.character(when)) when <- strptime(when, format) 

from

 if(is.character(when)) when <- strptime(when, format, tz="UTC") else when <- as.POSIXlt(when, tz="UTC") 

Now we can call sunPosition with:

 sunPosition(when = "2016-09-08 09:09:05-0500", format="%Y-%m-%d %H:%M:%S%z",lat = 43, long = -89) ## elevation azimuthJ azimuthC ##1 28.08683 110.4915 110.4915 

to get the same result. Note that we NEED TO specify the offset from UTC in the string literal and in format ( %z ) when calling sunPosition in this way.

With this change, sunPosition can be called using Sys.time() (I'm on the east coast):

 Sys.time() ##[1] "2016-09-08 12:42:08 EDT" sunPosition(Sys.time(),lat = 43, long = -89) ## elevation azimuthJ azimuthC ##1 49.24068 152.1195 152.1195 

which corresponds to the NOAA website

for Time Zone = -4 for EDT .

Answer 2

I think the problem is related to longitude. If I set longitude to 0 and latitude to my latitude and time before my time, I get the correct values ​​for altitude and azimuth.

 > time <- Sys.time() > time [1] "2016-09-08 12:07:35 CDT" > sunPosition(when = time, lat = 43, long = 0) elevation azimuthJ azimuthC 1 52.36687 184.1056 184.1056 

It seems to me that longitude is longitude relative to your position. I am not an expert on this topic, but in a sense it makes sense that longitude will be that way, since it does not have much effect on the position of the Sun. A person at a given latitude at a given local time will see the sun in the same position in the sky as someone else, located in a different longitude, but with the same latitude and local time (ignoring the complications of time zones, which are cropped areas with discrete boundaries and earth moving around the Sun).

I may not have read the question or function well enough, but unexpectedly, that longitude behaves this way.

---

edit: Reading @aichao's answer shows why setting latitude to zero and the time spent in local time works approximately. However, I do not think it will be very accurate.