Position of the Sun Given Time of Day, Latitude and Longitude

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

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.

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

I believe the line

if ($dayfrac < 0) $dayfrac += 1;

is in error. If you are before noon, you don't want to refer to the same time one day later, but instead you want to specify a time before noon, i.e. subtract from the julian date which represents noon.

Removing that line, your example date corresponds to the one computed using http://www.imcce.fr/en/grandpublic/temps/jour_julien.php, namely 2456536.9166666665. The resulting

$el = 67.775028608168
$az = 54.515316112281

looks pretty good to me. In particular, it agrees with the R run

elevation = 67.77503
azimuth = 54.51532

and also with what Stellarium says (although I quoted this incorrectly in a comment above):

Alt = 67°46'30" = 67.775
Az  = 54°30'60" = 45.5167

It also (almost) agrees with sunearthtools.com, so I guess you made a mistake when first entering the data there:

So I'd say that solves the problem.

How can we compute solar position at a given place on a given day and time?

It is an interesting problem and I think I have a good answer - well, at least starting points.

Check out the awesome astropy package. I believe you need to use the coordinates module.

Something along these lines:

import astropy.coordinates as coord
from astropy.time import Time
import astropy.units as u


loc = coord.EarthLocation(lon=0.1 * u.deg,
                          lat=51.5 * u.deg)
now = Time.now()

altaz = coord.AltAz(location=loc, obstime=now)
sun = coord.get_sun(now)

print(sun.transform_to(altaz).alt)

Here, we are getting the angle of the sun above the horizon for the 0.1 degrees longitude and 51.5 latitude location at the current time.

FYI, .zen would give you the zenith angle.

R function for position of sun giving unexpected results

There is nothing wrong with the core code found in the linked post if the input when is given in UTC. The confusion was that the OP entered the wrong Time Zone in the website for the 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 elevation of 38, which I think is correct.

The correct Time Zone to input into the NOAA website is -5 for CDT (see this website), which gives:

Calling sunPosition with the time adjusted 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 that is to replace the first line in sunPosition:

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

with

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

We can now 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 the format (%z) when calling sunPosition this way.

With this change sunPosition can be called with 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 matches the NOAA website

for Time Zone = -4 for EDT.

Time from sun's position and observer's location

Here's one (slightly rough) method, as the position calculators I've seen don't seem particularly easy to reverse-engineer for "time given a known date".

I took an existing sun position calculator, Sun Position from the file exchange. I then wrote two wrapper functions to this function to output info in a slightly more useful format for searching:

sun_hour →　calculates hourly (zero minute) positions across a given date.

sun_min → given a date and hour, calculates minutely positions one hour either side of thate date/hour.

Since the original function is quick, both of these take fractions of a second even on my not particularly good computer.

Then, I wrote a quick time location function along the following lines, using measured azimuth alone (along with time/location inputs giving the date and position in the formats required by sun_position).

time.min = 0;

% find list of azimuths, find closest match
[azlist, zlist, hours] = sun_hour(time, location);
n = find(abs(azlist-azimuth)==min(abs(azlist-azimuth)));
time.hour = hours(n);

% ditto, but around the time found at previous step
[azlist, zlist, hours, minutes] = sun_min(time, location);
n = find(abs(azlist-azimuth)==min(abs(azlist-azimuth)));

% output time
time.hour = hours(n);
time.min = minutes(n);

I haven't thoroughly tested it but it seems to work fine for a rough estimate (of course, assuming the accuracy of the original position calculator is good). I'm assuming that you don't need any more accuracy than minutely (unless you're managing to measure the position very accurately). You will need to be careful about how whichever calculator you end up using handles time-zones, whether the output is local time or not, etc.

Position of the Sun (azimuth) in Lua

I have found a way/hack to fix this issue.

Function is still the same:

-- solar altitude, azimuth (degrees)
function sunposition(latitude, longitude, time)
    time = time or os.time()
    if type(time) == 'table' then time = os.time(time) end

    local date = os.date('*t', time)
    local timezone = (os.time(date) - os.time(os.date('!*t', time))) / 3600
    if date.isdst then timezone = timezone + 1 end

    local utcdate = os.date('*t', time - timezone * 3600)
    local latrad = math.rad(latitude)
    local fd = (utcdate.hour + utcdate.min / 60 + utcdate.sec / 3600) / 24
    local g = (2 * math.pi / 365.25) * (utcdate.yday + fd)
    local d = math.rad(0.396372 - 22.91327 * math.cos(g) + 4.02543 * math.sin(g) - 0.387205 * math.cos(2 * g)
      + 0.051967 * math.sin(2 * g) - 0.154527 * math.cos(3 * g) + 0.084798 * math.sin(3 * g))
    local t = math.rad(0.004297 + 0.107029 * math.cos(g) - 1.837877 * math.sin(g)
      - 0.837378 * math.cos(2 * g) - 2.340475 * math.sin(2 * g))
    local sha = 2 * math.pi * (fd - 0.5) + t + math.rad(longitude)

    local sza = math.acos(math.sin(latrad) * math.sin(d) + math.cos(latrad) * math.cos(d) * math.cos(sha))
    local saa = math.acos((math.sin(d) - math.sin(latrad) * math.cos(sza)) / (math.cos(latrad) * math.sin(sza)))

    return 90 - math.deg(sza), math.deg(saa)
end

The added code that fixed the issue:

function getSunPos(lat, long, time)
    findTime = {}
    findTime.hour, findTime.min = time.hour, time.min
    fixedAzimuthLast, fixedAzimuth = 0, 0
    for i=0,23 do
        for j=0,59 do
            time.hour, time.min = i, j
            local altitude, azimuth = sunposition(lat, long, time)
            -- fix azimuth
            if fixedAzimuthLast < azimuth then 
                fixedAzimuthLast = azimuth
                fixedAzimuth = fixedAzimuthLast
            else
                fixedAzimuth = fixedAzimuthLast + (180 - azimuth)
            end
            -- find azimuth at target time
            if findTime.hour == i and findTime.min == j then
                -- final result
                return altitude, fixedAzimuth
            end
        end
    end
end

And finally to get the correct result:

lat, long = 45.327063, 14.442176
altitude, azimuth = getSunPos(lat, long, os.date('*t', os.time()))

That is it.
I would be happier with the complete function that does the math correctly but this will suffice. It works and was tested on 3 locations globally.

---

Related Topics