How to calculate sun position?

How to calculate sun position? - Printable Version

+- LogicMachine Forum (https://forum.logicmachine.net)
+-- Forum: LogicMachine eco-system (https://forum.logicmachine.net/forumdisplay.php?fid=1)
+--- Forum: Scripting (https://forum.logicmachine.net/forumdisplay.php?fid=8)
+--- Thread: How to calculate sun position? (/showthread.php?tid=161)

Pages: 1 2

RE: How to calculate sun position? - LeoSantaCruz - 17.04.2018

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

RE: How to calculate sun position? - Erwin van der Zwart - 17.04.2018

Hi,

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

BR,

Erwin

RE: How to calculate sun position? - LeoSantaCruz - 18.04.2018

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

RE: How to calculate sun position? - admin - 18.04.2018

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

RE: How to calculate sun position? - Trond Hoyem - 18.03.2020

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

RE: How to calculate sun position? - Daniel - 18.03.2020

do you have time set correctly in LM ?

RE: How to calculate sun position? - Trond Hoyem - 18.03.2020

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

RE: How to calculate sun position? - Daniel - 18.03.2020

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

RE: How to calculate sun position? - Trond Hoyem - 18.03.2020

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

RE: How to calculate sun position? - Erwin van der Zwart - 02.07.2020

(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

RE: How to calculate sun position? - admin - 02.07.2020

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