weblogpro/suntime.lua

   1 #!/usr/bin/lua
   2 --[[
   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
   7     tolerance in seconds.
   8
   9     Copyright (C) 2019 SatAgro Sp. z o.o. and contributors
  10     Copyright (C) 2022 Sven Hoexter <sven@stormbind.net>
  11
  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.
  16
  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
  20     for more details.
  21
  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/>.
  24 ]]
  25
  26 -- input values
  27 lon,lat = "3.1234","51.1234"
  28 tolerance = 4800
  29 dateTable = os.date("*t")
  30 curEpoch = os.time(dateTable)
  31 year = dateTable["year"]
  32 month = dateTable["month"]
  33 day = dateTable["day"]
  34
  35 function force_range(v, max)
  36     if( v < 0 ) then
  37         return (v + max)
  38     elseif( v >= max) then
  39         return (v - max)
  40     end
  41     return v
  42 end
  43
  44 function daysInMonth(year, month)
  45     month2days = {"31", "28", "31", "30", "31", "30", "31", "31", "30", "31", "30", "31"}
  46
  47     -- special case February http://www.henk-reints.nl/cal/gregcal.htm
  48     if(month == 2) then
  49         -- no century change and can be devided by 4
  50         if(((year % 100) ~= 0) and ((year % 4) == 0)) then
  51             month2days[2] = 29
  52         -- change of a century and can be devided by 400
  53         elseif(((year % 100) == 0) and ((year % 400) == 0)) then
  54             month2days[2] = 29
  55         end
  56     end
  57     return month2days[month]
  58 end
  59
  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
  63
  64      year: number
  65      month: number
  66      day: number
  67      lon: string, longitude of your location
  68      lat: string, latitude of your location
  69      getRiseTime: bool, true -> sunrise time, false -> sunset time
  70
  71      return value: time in Unix epoch
  72      ]]
  73 function calc_sun_time(year, month, day, lon, lat, getRiseTime)
  74     -- constants we need
  75     TO_RAD = (math.pi/180)
  76     zenith = 90.8
  77
  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
  83
  84     -- 2. convert the longitude to hour value and calculate an approximate time
  85     lonHour=lon/15
  86
  87     if(getRiseTime) then
  88         t = N + ((6 - lonHour) / 24)
  89     else
  90         t = N + ((18 - lonHour) / 24)
  91     end
  92
  93     -- 3. calculate the Sun's mean anomaly
  94     M = (0.9856 * t) - 3.289
  95
  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)
  99
 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)
 103
 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)
 108
 109     -- 5c. right ascension value needs to be converted into hours
 110     RA = RA / 15
 111
 112     -- 6. calculate the Sun's declination
 113     sinDec = 0.39782 * math.sin(TO_RAD * L)
 114     cosDec = math.cos(math.asin(sinDec))
 115
 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))
 118     if(cosH > 1) then
 119         -- return None -- The sun never rises on this location (on the specified date)
 120         print("Sun will never rise here and now.")
 121         os.exit(1)
 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.")
 125         os.exit(0)
 126     end
 127
 128     -- 7b. finish calculating H and convert into hours
 129     if(getRiseTime) then
 130         H = 360 - (1/TO_RAD) * math.acos(cosH)
 131     else
 132         H = (1/TO_RAD) * math.acos(cosH)
 133     end
 134     H = H / 15
 135
 136     -- 8. calculate local mean time of rising/setting
 137     T = H + RA - (0.06571 * t) - 6.622
 138
 139     -- 9. adjust back to UTC
 140     UT = T - lonHour
 141     UT = force_range(UT, 24)   -- UTC time in decimal format (e.g. 23.23)
 142
 143     -- 10. Return
 144     hr = force_range(math.floor(UT), 24)
 145     min = tonumber(string.format("%.0f", (UT - math.floor(UT))*60))
 146     if(min == 60) then
 147         hr = (hr + 1)
 148         min = 0
 149     end
 150
 151     -- 10. check corner case https://github.com/SatAgro/suntime/issues/1 and issue 8
 152     if(hr == 24) then
 153         hr = 0
 154         day = (day + 1)
 155
 156         if(day > daysInMonth(year, month)) then
 157             day = 1
 158             month = (month + 1)
 159
 160             if(month > 12) then
 161                 month = 1
 162                 year = (year + 1)
 163             end
 164         end
 165     end
 166
 167     return os.time{year=year, month=month, day=day, hour=hr, min=math.floor(min)}
 168 end
 169
 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)
 173
 174 -- debugging aide
 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))
 180 end
 181
 182 -- return 0 in case we should've daylight, otherwise return 1
 183 if((curEpoch >= sunrise) and (curEpoch <= sunset)) then
 184     os.exit(0)
 185 else
 186     os.exit(1)
 187 end