How to calculate sun position?

[buuuudzik]

Hello

Do you have some algorithm to calculate a sun position?

[admin]

This function should work, though it was not tested much:

Code:

12345678910111213141516171819202122232425262728

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

[buuuudzik]

This function should work, though it was not tested much:

Code:

12345678910111213141516171819202122232425262728

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

It works. Thank you

I create a library with this function in User library and in resident script I add this script:

Code:

1234567

latitude = 50.990 longtitude = -23.179 time = nil altitude, azimuth = sunposition(latitude, longitude, time) log(altitude, azimuth)

Result is very close to the results from online calculators.

[admin]

Yes, I think it was ported from JavaScript based on one of online calculators.

[oab_nor]

Hi, while we’re on this topic. 

Has anyone tried to work out an algorithm to calculate slate position of venetian blinds? Given the orientation of a window/blind, and the solar altitude and azimuth (given in post above). 

What is the optimal slate position to avoid direct sunlight? 

I’we been playing with some calculations. Even tried to “reverse engineer” observations from a major brand “KNX Wetterstation”, but never really got hold of it. Anyone with good experiences on this?

Any suggestions on how to approach this is welcome. 
E.g. Is there one formula? 
Should i use a lookup table with different altitude- intervals for different azimuth “ranges”? 
Should i use a combination?

[admin]

For starters, you can try with just using altitude since blinds only move on one axis. Blinds should be perpendicular to current sun position. Then you can tweak it by checking azimuth so you don't adjust blinds when there's no direct sunlight.

[oab_nor]

For starters, you can try with just using altitude since blinds only move on one axis. Blinds should be perpendicular to current sun position. Then you can tweak it by checking azimuth so you don't adjust blinds when there's no direct sunlight.

Tanks, maybe  i'm just making it too complex.  
There should, of course, be a 1:1 relation between sun altitude and  slate position (e.g. 0-100%) – regardless of azimuth and orientation of the blinds/window.

[buuuudzik]

Can I add this library to LM?

https://github.com/mourner/suncalc

[admin]

You cannot use this library directly, but rewriting it in Lua should be quite easy.

[Trond Hoyem]

This function should work, though it was not tested much:

Code:

12345678910111213141516171819202122232425262728

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

Hi

I used this today. What I see is that when the azimuth get to 180, it starts counting down rather than going on to the full 360... Any ideas why?

[admin]

Try this one instead:

Code:

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576

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

[Trond Hoyem]

Try this one instead:

Code:

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576

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

Thanks. That one worked.

Trond

[Auzzii]

Try this one instead:

Code:

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576

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

Thanks. That one worked.

Trond

Hi,

Maybe someone have example code?

I need to get sun position every 5 min.
1) Elevation degree (0-100%)
2) Azimuth (0-360*)


thx

[Erwin van der Zwart]

Hi,

Use function suncalc() like this:

Code:

1234567

lat = 52.5167747 -- Zwolle - Netherlands lng = 6.0830219 -- Zwolle - Netherlands time = os.time() alt, az = suncalc(lat, lng, time) grp.update('33/1/1', alt) grp.update('33/1/2', az)

BR,

Erwin

[Auzzii]

Hi,

Use function suncalc() like this:

Code:

1234567

lat = 52.5167747 -- Zwolle - Netherlands lng = 6.0830219 -- Zwolle - Netherlands time = os.time() alt, az = suncalc(lat, lng, time) grp.update('33/1/1', alt) grp.update('33/1/2', az)

BR,

Erwin

Hi,

I get error:

Common functions:212: attempt to index global 'cec' (a nil value)
stack traceback:
Common functions:212: in main chunk

[Erwin van der Zwart]

Hi,

Did you placed the function inside common functions? Do you do the call from a resident script?

I runned the function and the call from the same resident script and that worked perfect..

BR,

Erwin

[Auzzii]

Hi,

Did you placed the function inside common functions? Do you do the call from a resident script?

I runned the function and the call from the same resident script and that worked perfect..

BR,

Erwin

After fw upgrade and  factory reset all works fine.

thank you Erwin

[buuuudzik]

What about time zone and time shift? Please be carefull especially with the time shift because after it your calculations can be offset by 1 hour.

Code:

123

isdst = os.date('*t').isdst if isdst then d = d - 3600 end

[Erwin van der Zwart]

Hi,

I don't think that is needed, i send the raw unix timestamp and time extracted from the raw stamp already includes dst..

BR,

Erwin

[buuuudzik]

I've checked this and you're right, raw unix time is in UTC. I am using other algorithm and on it I must add the info about time zone and I must consider time shift