This forum uses cookies
This forum makes use of cookies to store your login information if you are registered, and your last visit if you are not. Cookies are small text documents stored on your computer; the cookies set by this forum can only be used on this website and pose no security risk. Cookies on this forum also track the specific topics you have read and when you last read them. Please confirm whether you accept or reject these cookies being set.

How to calculate sun position?
#21
(11.05.2017, 07:49)admin Wrote: 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?
Reply
#22
Hi,

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

BR,

Erwin
Reply
#23
(17.04.2018, 16:29)Erwin van der Zwart Wrote: 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.
Reply
#24
Original JavaScript code is under BSD license: https://github.com/mourner/suncalc
Reply
#25
(11.05.2017, 07:49)admin Wrote: 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.
There are 10 kinds of people in the world; those who can read binary and those who don't  Cool
Reply
#26
do you have time set correctly in LM ?
Reply
#27
(18.03.2020, 13:27)Daniel. Wrote: 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.
There are 10 kinds of people in the world; those who can read binary and those who don't  Cool
Reply
#28
Post your full script I will check on both devices but there can't be any difference
Reply
#29
(18.03.2020, 13:33)Daniel. Wrote: 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....
There are 10 kinds of people in the world; those who can read binary and those who don't  Cool
Reply
#30
(11.05.2017, 07:49)admin Wrote: 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
Reply
#31
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
Reply


Forum Jump: