3 Utility to let you figure out if we should've daylight now.
4 Usage: ./suntime.lua [debug]
5 Exit code is either 0, we have daylight - or 1, no daylight.
6 Adjust lon, lat, tolerance for your location and the required
9 Copyright (C) 2019 SatAgro Sp. z o.o. and contributors
10 Copyright (C) 2022 Sven Hoexter <sven@stormbind.net>
12 This program is free software: you can redistribute it and/or modify
13 it under the terms of the GNU Lesser General Public License as published
14 by the Free Software Foundation, either version 3 of the License, or
15 (at your option) any later version.
17 This program is distributed in the hope that it will be useful, but
18 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
19 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
22 You should have received a copy of the GNU Lesser General Public License
23 along with this program. If not, see <https://www.gnu.org/licenses/>.
27 lon,lat = "3.1234","51.1234"
29 dateTable = os.date("*t")
30 curEpoch = os.time(dateTable)
31 year = dateTable["year"]
32 month = dateTable["month"]
33 day = dateTable["day"]
35 function force_range(v, max)
38 elseif( v >= max) then
44 function daysInMonth(year, month)
45 month2days = {"31", "28", "31", "30", "31", "30", "31", "31", "30", "31", "30", "31"}
47 -- special case February http://www.henk-reints.nl/cal/gregcal.htm
49 -- no century change and can be devided by 4
50 if(((year % 100) ~= 0) and ((year % 4) == 0)) then
52 -- change of a century and can be devided by 400
53 elseif(((year % 100) == 0) and ((year % 400) == 0)) then
57 return month2days[month]
60 --[[ Ported from https://pypi.org/project/suntime/
61 which is based on https://stackoverflow.com/questions/19615350/
62 which took the algorithm from https://web.archive.org/web/20161022214335/https://williams.best.vwh.net/sunrise_sunset_algorithm.htm
67 lon: string, longitude of your location
68 lat: string, latitude of your location
69 getRiseTime: bool, true -> sunrise time, false -> sunset time
71 return value: time in Unix epoch
73 function calc_sun_time(year, month, day, lon, lat, getRiseTime)
75 TO_RAD = (math.pi/180)
78 -- 1. first calculate the day of the year
79 N1 = math.floor(275 * month / 9)
80 N2 = math.floor((month + 9 ) / 12)
81 N3 = ( 1 + math.floor((year - 4 * math.floor(year / 4) + 2) / 3))
82 N = N1 - (N2 * N3) + day - 30
84 -- 2. convert the longitude to hour value and calculate an approximate time
88 t = N + ((6 - lonHour) / 24)
90 t = N + ((18 - lonHour) / 24)
93 -- 3. calculate the Sun's mean anomaly
94 M = (0.9856 * t) - 3.289
96 -- 4. calculate the Sun's true longitude
97 L = M + (1.916 * math.sin(TO_RAD*M)) + (0.020 * math.sin(TO_RAD * 2 * M)) + 282.634
98 L = force_range(L, 360 ) -- NOTE: L adjusted into the range [0,360)
100 -- 5a. calculate the Sun's right ascension
101 RA = (1/TO_RAD) * math.atan(0.91764 * math.tan(TO_RAD*L))
102 RA = force_range(RA, 360 ) -- NOTE: RA adjusted into the range [0,360)
104 -- 5b. right ascension value needs to be in the same quadrant as L
105 Lquadrant = (math.floor( L/90)) * 90
106 RAquadrant = (math.floor(RA/90)) * 90
107 RA = RA + (Lquadrant - RAquadrant)
109 -- 5c. right ascension value needs to be converted into hours
112 -- 6. calculate the Sun's declination
113 sinDec = 0.39782 * math.sin(TO_RAD * L)
114 cosDec = math.cos(math.asin(sinDec))
116 -- 7a. calculate the Sun's local hour angle
117 cosH = (math.cos(TO_RAD * zenith) - (sinDec * math.sin(TO_RAD * lat))) / (cosDec * math.cos(TO_RAD * lat))
119 -- return None -- The sun never rises on this location (on the specified date)
120 print("Sun will never rise here and now.")
122 elseif(cosH < -1) then
123 -- return None -- The sun never sets on this location (on the specified date)
124 print("Sun will never rise here and now.")
128 -- 7b. finish calculating H and convert into hours
130 H = 360 - (1/TO_RAD) * math.acos(cosH)
132 H = (1/TO_RAD) * math.acos(cosH)
136 -- 8. calculate local mean time of rising/setting
137 T = H + RA - (0.06571 * t) - 6.622
139 -- 9. adjust back to UTC
141 UT = force_range(UT, 24) -- UTC time in decimal format (e.g. 23.23)
144 hr = force_range(math.floor(UT), 24)
145 min = tonumber(string.format("%.0f", (UT - math.floor(UT))*60))
151 -- 10. check corner case https://github.com/SatAgro/suntime/issues/1 and issue 8
156 if(day > daysInMonth(year, month)) then
167 return os.time{year=year, month=month, day=day, hour=hr, min=math.floor(min)}
170 -- calculate sunrise and sunset time adjusted by tolerance seconds
171 sunrise = (calc_sun_time(year, month, day, lon, lat, true) + tolerance)
172 sunset = (calc_sun_time(year, month, day, lon, lat, false) - tolerance)
175 if(arg[1] == "debug") then
176 print("Expected sunrise for today in UTC seconds adjusted +", tolerance)
177 print(os.date("%c", sunrise))
178 print("Expected sunset for today in UTC seconds adjusted -", tolerance)
179 print(os.date("%c", sunset))
182 -- return 0 in case we should've daylight, otherwise return 1
183 if((curEpoch >= sunrise) and (curEpoch <= sunset)) then