From 28b33f9ee6d08b99e99b8158bb593c7dfd7e8365 Mon Sep 17 00:00:00 2001 From: Sven Hoexter Date: Wed, 2 Feb 2022 22:18:14 +0100 Subject: [PATCH] Add Lua script to check if we should have daylight or not --- home/suntime.lua | 171 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100755 home/suntime.lua diff --git a/home/suntime.lua b/home/suntime.lua new file mode 100755 index 0000000..3dd53f6 --- /dev/null +++ b/home/suntime.lua @@ -0,0 +1,171 @@ +#!/usr/bin/lua +--[[ + Utility to let you figure out if we should've daylight now. + Usage: ./suntime.lua [debug] + Exit code is either 0, we have daylight - or 1, no daylight. + Adjust lon, lat, tolerance for your location and the required + tolerance in seconds. +]] + +-- input values +lon,lat = "3.1234","51.1234" +tolerance = 4800 +dateTable = os.date("*t") +curEpoch = os.time(dateTable) +year = dateTable["year"] +month = dateTable["month"] +day = dateTable["day"] + +function force_range(v, max) + if( v < 0 ) then + return (v + max) + elseif( v >= max) then + return (v - max) + end + return v +end + +function daysInMonth(year, month) + month2days = {"31", "28", "31", "30", "31", "30", "31", "31", "30", "31", "30", "31"} + + -- special case February http://www.henk-reints.nl/cal/gregcal.htm + if(month == 2) then + -- no century change and can be devided by 4 + if(((year % 100) ~= 0) and ((year % 4) == 0)) then + month2days[2] = 29 + -- change of a century and can be devided by 400 + elseif(((year % 100) == 0) and ((year % 400) == 0)) then + month2days[2] = 29 + end + end + return month2days[month] +end + +--[[ Ported from https://pypi.org/project/suntime/ + which is based on https://stackoverflow.com/questions/19615350/calculate-sunrise-and-sunset-times-for-a-given-gps-coordinate-within-postgresql + which took the algorithm from https://web.archive.org/web/20161022214335/https://williams.best.vwh.net/sunrise_sunset_algorithm.htm + + year: number + month: number + day: number + lon: string, longitude of your location + lat: string, latitude of your location + getRiseTime: bool, true -> sunrise time, false -> sunset time + + return value: time in Unix epoch + ]] +function calc_sun_time(year, month, day, lon, lat, getRiseTime) + -- constants we need + TO_RAD = (math.pi/180) + zenith = 90.8 + + -- 1. first calculate the day of the year + N1 = math.floor(275 * month / 9) + N2 = math.floor((month + 9 ) / 12) + N3 = ( 1 + math.floor((year - 4 * math.floor(year / 4) + 2) / 3)) + N = N1 - (N2 * N3) + day - 30 + + -- 2. convert the longitude to hour value and calculate an approximate time + lonHour=lon/15 + + if(getRiseTime) then + t = N + ((6 - lonHour) / 24) + else + t = N + ((18 - lonHour) / 24) + end + + -- 3. calculate the Sun's mean anomaly + M = (0.9856 * t) - 3.289 + + -- 4. calculate the Sun's true longitude + L = M + (1.916 * math.sin(TO_RAD*M)) + (0.020 * math.sin(TO_RAD * 2 * M)) + 282.634 + L = force_range(L, 360 ) -- NOTE: L adjusted into the range [0,360) + + -- 5a. calculate the Sun's right ascension + RA = (1/TO_RAD) * math.atan(0.91764 * math.tan(TO_RAD*L)) + RA = force_range(RA, 360 ) -- NOTE: RA adjusted into the range [0,360) + + -- 5b. right ascension value needs to be in the same quadrant as L + Lquadrant = (math.floor( L/90)) * 90 + RAquadrant = (math.floor(RA/90)) * 90 + RA = RA + (Lquadrant - RAquadrant) + + -- 5c. right ascension value needs to be converted into hours + RA = RA / 15 + + -- 6. calculate the Sun's declination + sinDec = 0.39782 * math.sin(TO_RAD * L) + cosDec = math.cos(math.asin(sinDec)) + + -- 7a. calculate the Sun's local hour angle + cosH = (math.cos(TO_RAD * zenith) - (sinDec * math.sin(TO_RAD * lat))) / (cosDec * math.cos(TO_RAD * lat)) + if(cosH > 1) then + -- return None -- The sun never rises on this location (on the specified date) + print("Sun will never rise here and now.") + os.exit() + elseif(cosH < -1) then + -- return None -- The sun never sets on this location (on the specified date) + print("Sun will never rise here and now.") + os.exit() + end + + -- 7b. finish calculating H and convert into hours + if(getRiseTime) then + H = 360 - (1/TO_RAD) * math.acos(cosH) + else + H = (1/TO_RAD) * math.acos(cosH) + end + H = H / 15 + + -- 8. calculate local mean time of rising/setting + T = H + RA - (0.06571 * t) - 6.622 + + -- 9. adjust back to UTC + UT = T - lonHour + UT = force_range(UT, 24) -- UTC time in decimal format (e.g. 23.23) + + -- 10. Return + hr = force_range(math.floor(UT), 24) + min = tonumber(string.format("%.0f", (UT - math.floor(UT))*60)) + if(min == 60) then + hr = (hr + 1) + min = 0 + end + + -- 10. check corner case https://github.com/SatAgro/suntime/issues/1 and issue 8 + if(hr == 24) then + hr = 0 + day = (day + 1) + + if(day > daysInMonth(year, month)) then + day = 1 + month = (month + 1) + + if(month > 12) then + month = 1 + year = (year + 1) + end + end + end + + return os.time{year=year, month=month, day=day, hour=hr, min=math.floor(min)} +end + +-- calculate sunrise and sunset time adjusted by tolerance seconds +sunrise = (calc_sun_time(year, month, day, lon, lat, true) + tolerance) +sunset = (calc_sun_time(year, month, day, lon, lat, false) - tolerance) + +-- debugging aide +if(arg[1] == "debug") then + print("Expected sunrise for today in UTC seconds adjusted +", tolerance) + print(os.date("%c", sunrise)) + print("Expected sunset for today in UTC seconds adjusted -", tolerance) + print(os.date("%c", sunset)) +end + +-- return 0 in case we should've daylight, otherwise return 1 +if((curEpoch >= sunrise) and (curEpoch <= sunset)) then + os.exit(0) +else + os.exit(1) +end \ No newline at end of file -- 2.39.5