Position of the Sun (azimuth) in Lua

Question

There is only one function in LUA I could find online, but it gives wrong values (measured with professional online tools).

It appears that from the sunrise till some time after the noon the math works, but after, the Sun's angle goes back to the sunrise position. Should be from 106° to 253°, currently it's from 106° to ~180° to 106°.

Function I'm using:

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

Example request:

lat, long = 45.327063, 14.442176 -- Rijeka, Croatia
time = {year=2016, month=2, day=17, hour=17, min=30} -- end of the day
altitude, azimuth = sunposition(lat, long, time)

Result is:

• -0.1 degrees in altitude
• 106 degrees in azimuth.

Result should be:

• -0.1 degrees in altitude
• 253 degrees in azimuth.

I have found multiple solutions in other programming languages and even tried to rewrite in Lua but without any success. Too complex math behind the solution.

I'm using it for my Corona SDK app that will show position of the Sun relative to the device. The only solution that currently works is a PHP or Javascript script that my app can ask via API call over the Internet but I would really like to avoid that.

I'm extremely grateful for any help from the community. Thank you and love you folks! :)

Answer 1

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.

Answer 2

Arccos of an angle equals to 360-angle as their cosine values is equal. You can simply return 360-angle.

Answer 3

thanks for sharing this.

I have seen that the modified function does not work well in the time slot between UTC day change 00:00 time and the local time day change, of the desired location.

To get the displayed azimuth value to match the navigation standard with a linear displacement of 15º/h from 0º to 360º:

• The Solar Azimuth will be to the South (180º) at solar noon in the Northern Hemisphere
• The Solar Azimuth will be to the North (0º) at solar noon in the Southern Hemisphere
• The sunrise always the azimuth will show the east (90º +/15º), 90º in equinoxes
• The Sunrise always the azimuth will show the West (270º +/15º),270º in equinoxes

The calculation also needs a different correction for E and W longitudes and N and S latitudes.

The calculation function makes 2 direction changes in the calculated azimuth, one at solar noon and another 12 hours later at night. Instead of calculating 0º to 360º, it calculates 0º to 180º and returns from 180º to 0º.

This is the code for the function, which only needs send it the latitude, longitude, and UTC time, with os.time()

The azimuth value is calculated 2 steps:

• The first with the actual UTC time and the values returned altitude and azimuth are saved
• The second step is sent to the function the UTC time - 60 sec
• The returned azimuth value confirm the direction of the azimuth increase (increasing, UP or decreasing DOWN).
• Depending on whether the trend is UP or DOWN, the azimuth of the UTC time is corrected by subtracting it from 360º depending on whether the location is in the NE, NW, SE or SW quadrant.

I have tested and checked it with 4 locations in the four quadrants, Beijing (NE), Melbourne (SE), New York (NW) and Rio de Janeiro (SW)

It works with a precision of +/- 0.5º with respect to an online calculator: enter link description here

This is the final code

    local 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 utcdate = os.date('*t', time)
    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)))

    print("<<<<<<<<<<<< utcdate.hour",utcdate.hour,"utcdate.min",utcdate.min)
    print("<<<<<< sza",sza,"<<<<<< saa",saa)
    print("<<<<< longitude",longitude,"Altitude",math.deg(sza),"Azimuth",math.deg(saa))

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


local function getSunPos(lat, long, time)

    -- calculate UTC altitude and azimuth
    local altitude, azimuth = sunposition(lat, long, time)

    -- Check previous minute values to detect azimuth direction change 0º and 180º
    time = time - 60
    local previuos_altitude, previuos_azimuth = sunposition(lat, long, time)

    print("previuos_azimuth",previuos_azimuth,"<<<< azimuth",azimuth)
    if previuos_azimuth < azimuth then
       print("<<<<<<<<<<<<< TREND:", "UP")
    else
       print("<<<<<<<<<<<<< TREND:", "DOWN")
    end
    if previuos_azimuth < azimuth then -- calcualted azimuth trend UP >>>>>>>>>>>>>>>>>>>
      if long < 0 then
        if lat > 0 then
            return altitude, azimuth
        else
            return altitude, 360 - azimuth
        end
      else ------------------- long > 0
        if lat > 0 then
            return altitude, azimuth
        else
            return altitude, 360 - azimuth
        end
      end
    else -- calcualted azimuth trend DOWN >>>>>>>>>>>>>>>>>>>
      if long < 0 then
        if lat > 0 then
            return altitude,360 - azimuth
        else
            return altitude, azimuth
        end
      else ------------------- long > 0
        if lat > 0 then
            return altitude, 360 - azimuth
        else
            return altitude, azimuth
        end
      end
    end
end


local function calculate_altitude_azimuth(lat, long, time)
    local altitude, azimuth = getSunPos(lat, long, os.time())
end

Answer 4

Based on the NOAA constants

lat, long are the format from google maps... datetime is UTC time (nil will grab from the system) tzoffset is hours from UTC without daylight time) e.g. -5 is EST

function noaa_sun_position(lat, long, datetime, tzoffset)
--https://gml.noaa.gov/grad/solcalc/solareqns.PDF
--
datetime = datetime("!*t")
tzoffset = tzoffset or -5
local datetime= os.date("*t",os.time(datetime) - 5 * 60 * 60)
local day_of_year = os.date("%j", os.time(datetime))
local fract_year = (2 * math.pi / 365) * ( day_of_year - 1 + (os.date("*t").hour - 5 - 12) / 24)

local eqtime = 229.18*
                 (0.000075+
                 0.001868*math.cos(fract_year)
               - 0.032077 * math.sin(fract_year)
               - 0.014615 * math.cos(2*fract_year)
               - 0.040849 * math.sin(2*fract_year))
local decl = 0.006918
        - 0.399912*math.cos(fract_year )
        + 0.070257*math.sin(fract_year )
        - 0.006758*math.cos(2*fract_year )
        + 0.000907*math.sin(2*fract_year )
        -0.002697*math.cos(3*fract_year )
        + 0.00148*math.sin (3*fract_year )
local time_offset = eqtime + 4 * long - 60 * -5
local tst = datetime.hour*60 + datetime.min  + datetime.sec/60 + time_offset
local ha = (tst / 4) -180
local cos_zenith =  math.sin(math.rad(lat))*math.sin(decl) + math.cos(math.rad(lat))* math.cos(decl) *math.cos(math.rad(ha))
--local azimuth = (math.sin(math.rad(lat) * cos_zenith - math.sin(decl))) / (math.cos(math.rad(lat) * math.sin(math.acos(cos_zenith))) )
local azimuth = math.atan2(-math.cos(math.rad(decl)) * math.sin(math.rad(ha)), math.sin(math.rad(decl)) * math.cos(math.rad(lat)) - math.cos(math.rad(decl)) * math.sin(math.rad(lat)) * math.cos(math.rad(ha)))
return ((math.deg(azimuth) > 0 and math.deg(azimuth)) or ( math.deg(azimuth) + 360 )) , 90-math.deg(math.acos(cos_zenith)) , ha, decl
end