How to calculate sun position?

[LeoSantaCruz]

Try this one instead:

Code:

function suncalc(lat, lng, time)
 -- sun calculations are based on http://aa.quae.nl/en/reken/zonpositie.html formulas
 local PI   = math.pi
 local sin  = math.sin
 local cos  = math.cos
 local tan  = math.tan
 local asin = math.asin
 local atan = math.atan2
 local acos = math.acos
 local deg  = math.deg
 local rad  = PI / 180
 local e = rad * 23.4397 -- obliquity of the Earth
 local daysec = 60 * 60 * 24
 local J1970 = 2440588
 local J2000 = 2451545

 local function toDays(time)
   return time / daysec - 0.5 + J1970 - J2000
 end

 local function rightAscension(l, b)
   return atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l))
 end

 local function declination(l, b)
   return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l))
 end

 local function azimuth(H, phi, dec)
   return atan(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi))
 end

 local function altitude(H, phi, dec)
   return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H))
 end

 local function siderealTime(d, lw)
   return rad * (280.16 + 360.9856235 * d) - lw
 end

 local function astroRefraction(h)
   if h < 0 then -- the following formula works for positive altitudes only.
     h = 0 -- if h = -0.08901179 a div/0 would occur.
   end

   -- formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
   -- 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad:
   return 0.0002967 / math.tan(h + 0.00312536 / (h + 0.08901179))
 end

 -- general sun calculations
 local function solarMeanAnomaly(d)
   return rad * (357.5291 + 0.98560028 * d)
 end

 local function eclipticLongitude(M)
   local C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)) -- equation of center
   local P = rad * 102.9372 -- perihelion of the Earth
   return M + C + P + PI
 end

 local function sunCoords(d)
   local M = solarMeanAnomaly(d)
   local L = eclipticLongitude(M)
   return declination(L, 0), rightAscension(L, 0)
 end

 local lw      = rad * -lng
 local phi     = rad * lat
 local d       = toDays(time or os.time())
 local dec, ra = sunCoords(d)
 local H       = siderealTime(d, lw) - ra

 local alt, az = altitude(H, phi, dec), azimuth(H, phi, dec)
 return deg(alt), 180 + deg(az)
end

Where can i get the License for this code?

[Erwin van der Zwart]

Hi,

You don't need a license, it's just a mathematical function..

BR,

Erwin

[LeoSantaCruz]

Hi,

You don't need a license, it's just a mathematical function..

BR,

Erwin

So I can just Reference to this site? I am asking because i am using this function for a Project.

[admin]

Original JavaScript code is under BSD license: https://github.com/mourner/suncalc

[Trond Hoyem]

Try this one instead:

Code:

function suncalc(lat, lng, time)
  -- sun calculations are based on http://aa.quae.nl/en/reken/zonpositie.html formulas
  local PI   = math.pi
  local sin  = math.sin
  local cos  = math.cos
  local tan  = math.tan
  local asin = math.asin
  local atan = math.atan2
  local acos = math.acos
  local deg  = math.deg
  local rad  = PI / 180
  local e = rad * 23.4397 -- obliquity of the Earth
  local daysec = 60 * 60 * 24
  local J1970 = 2440588
  local J2000 = 2451545

  local function toDays(time)
    return time / daysec - 0.5 + J1970 - J2000
  end

  local function rightAscension(l, b)
    return atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l))
  end

  local function declination(l, b)
    return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l))
  end

  local function azimuth(H, phi, dec)
    return atan(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi))
  end

  local function altitude(H, phi, dec)
    return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H))
  end

  local function siderealTime(d, lw)
    return rad * (280.16 + 360.9856235 * d) - lw
  end

  local function astroRefraction(h)
    if h < 0 then -- the following formula works for positive altitudes only.
      h = 0 -- if h = -0.08901179 a div/0 would occur.
    end

    -- formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
    -- 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad:
    return 0.0002967 / math.tan(h + 0.00312536 / (h + 0.08901179))
  end

  -- general sun calculations
  local function solarMeanAnomaly(d)
    return rad * (357.5291 + 0.98560028 * d)
  end

  local function eclipticLongitude(M)
    local C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)) -- equation of center
    local P = rad * 102.9372 -- perihelion of the Earth
    return M + C + P + PI
  end

  local function sunCoords(d)
    local M = solarMeanAnomaly(d)
    local L = eclipticLongitude(M)
    return declination(L, 0), rightAscension(L, 0)
  end

  local lw      = rad * -lng
  local phi     = rad * lat
  local d       = toDays(time or os.time())
  local dec, ra = sunCoords(d)
  local H       = siderealTime(d, lw) - ra

  local alt, az = altitude(H, phi, dec), azimuth(H, phi, dec)
  return deg(alt), 180 + deg(az)
end

I have used this for a few different applications on SL, and it works. But today I implementet this on a LM5, and i get all wrong values from it.

When I log the returned values I get; 

* arg: 1
  * number: 383
* arg: 2
  * number: 1112

Any ideas why? I even tested to enter it at the SL I use for testing at my office, and then I get the correct values.
The firmware of the LM5 is 20180822.

[Daniel]

do you have time set correctly in LM ?

[Trond Hoyem]

do you have time set correctly in LM ?

Yes, first thing I checked.

But anyway, the values should still be within 360 degrees, as that would be the highest possible Azimuth, and highest possible Elevation should be 90 degrees.

[Daniel]

Post your full script I will check on both devices but there can't be any difference

[Trond Hoyem]

Post your full script I will check on both devices but there can't be any difference

I disabled the script, and then enabled it again, and now it works....

[Erwin van der Zwart]

Try this one instead:

Code:

function suncalc(lat, lng, time)
  -- sun calculations are based on http://aa.quae.nl/en/reken/zonpositie.html formulas
  local PI   = math.pi
  local sin  = math.sin
  local cos  = math.cos
  local tan  = math.tan
  local asin = math.asin
  local atan = math.atan2
  local acos = math.acos
  local deg  = math.deg
  local rad  = PI / 180
  local e = rad * 23.4397 -- obliquity of the Earth
  local daysec = 60 * 60 * 24
  local J1970 = 2440588
  local J2000 = 2451545

  local function toDays(time)
    return time / daysec - 0.5 + J1970 - J2000
  end

  local function rightAscension(l, b)
    return atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l))
  end

  local function declination(l, b)
    return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l))
  end

  local function azimuth(H, phi, dec)
    return atan(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi))
  end

  local function altitude(H, phi, dec)
    return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H))
  end

  local function siderealTime(d, lw)
    return rad * (280.16 + 360.9856235 * d) - lw
  end

  local function astroRefraction(h)
    if h < 0 then -- the following formula works for positive altitudes only.
      h = 0 -- if h = -0.08901179 a div/0 would occur.
    end

    -- formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
    -- 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad:
    return 0.0002967 / math.tan(h + 0.00312536 / (h + 0.08901179))
  end

  -- general sun calculations
  local function solarMeanAnomaly(d)
    return rad * (357.5291 + 0.98560028 * d)
  end

  local function eclipticLongitude(M)
    local C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)) -- equation of center
    local P = rad * 102.9372 -- perihelion of the Earth
    return M + C + P + PI
  end

  local function sunCoords(d)
    local M = solarMeanAnomaly(d)
    local L = eclipticLongitude(M)
    return declination(L, 0), rightAscension(L, 0)
  end

  local lw      = rad * -lng
  local phi     = rad * lat
  local d       = toDays(time or os.time())
  local dec, ra = sunCoords(d)
  local H       = siderealTime(d, lw) - ra

  local alt, az = altitude(H, phi, dec), azimuth(H, phi, dec)
  return deg(alt), 180 + deg(az)
end

Hi Admin,

On what calculation is suncalc (above) based? The Azimuth is 28 degrees wrong, i just tested function suncalc vs sunposition and these are the results:

log(suncalc(52.5167747, 6.0830219, os.time()) )

log(sunposition(52.5167747, 6.0830219, os.time()) )

Suncalc 02.07.2020 14:21:25
* arg: 1 Altitude
  * number: 59.4576360078012
* arg: 2 Azimuth
  * number: 198.878804294899

Sunposition 02.07.2020 14:21:25
* arg: 1 Altitude
  * number: 59.3596373934031
* arg: 2 Azimuth
  * number: 160.948315437194 


Any idea what might be wrong in the function suncalc?

BR,

Erwin

[admin]

Hi Erwin, looks like the sunpositon function is wrong. I've checked several online calculators and they give a result similar to what suncalc produces. Suncalc formulas are taken from here: https://github.com/mourner/suncalc