From c676aa6f0d259a9cb5803e5a9a5fc911ea5abb38 Mon Sep 17 00:00:00 2001
From: zwim <36999612+zwim@users.noreply.github.com>
Date: Sun, 5 Dec 2021 09:52:02 +0100
Subject: [PATCH] [autowarmth, plugin] use date time widget, optimize (#8502)

---
 plugins/autowarmth.koplugin/main.lua    |  28 +--
 plugins/autowarmth.koplugin/suntime.lua | 322 +++++++++++++++++-------
 2 files changed, 240 insertions(+), 110 deletions(-)

diff --git a/plugins/autowarmth.koplugin/main.lua b/plugins/autowarmth.koplugin/main.lua
index 7eaf62850..7820f89e5 100644
--- a/plugins/autowarmth.koplugin/main.lua
+++ b/plugins/autowarmth.koplugin/main.lua
@@ -7,6 +7,7 @@ Plugin for setting screen warmth based on the sun position and/or a time schedul
 local Device = require("device")
 
 local ConfirmBox = require("ui/widget/confirmbox")
+local DateTimeWidget = require("ui/widget/datetimewidget")
 local DoubleSpinWidget = require("/ui/widget/doublespinwidget")
 local DeviceListener = require("device/devicelistener")
 local Dispatcher = require("dispatcher")
@@ -515,6 +516,7 @@ function AutoWarmth:getLocationMenu()
         callback = function(touchmenu_instance)
             UIManager:show(SpinWidget:new{
                 title_text = _("Altitude"),
+                info_text = _("Enter the altitude in meters above sea level."),
                 value = self.altitude,
                 value_min = -100,
                 value_max = 15000, -- intercontinental flight
@@ -577,24 +579,15 @@ function AutoWarmth:getScheduleMenu()
                     hh = math.floor(self.scheduler_times[num])
                     mm = math.floor(frac(self.scheduler_times[num]) * 60 + 0.5)
                 end
-                UIManager:show(DoubleSpinWidget:new{
+                UIManager:show(DateTimeWidget:new{
                     title_text = _("Set time"),
-                    left_text = _("HH"),
-                    left_value = hh,
-                    left_min = 0,
-                    left_max = 23,
-                    left_step = 1,
-                    left_hold_step = 3,
-                    left_wrap = true,
-                    right_text = _("MM"),
-                    right_value = mm,
-                    right_min = 0,
-                    right_max = 59,
-                    right_step = 1,
-                    right_hold_step = 5,
-                    right_wrap = true,
-                    callback = function(left, right)
-                        local new_time = left + right / 60
+                    info_text = _("Enter time in hours and minutes."),
+                    is_date = false,
+                    hour = hh,
+                    min = mm,
+                    ok_text = _("Set time"),
+                    callback = function(time)
+                        local new_time = time.hour + time.min / 60
                         local function get_valid_time(n, dir)
                             for i = n+dir, dir > 0 and midnight_index or 1, dir do
                                 if self.scheduler_times[i] then
@@ -691,6 +684,7 @@ function AutoWarmth:getWarmthMenu()
                 if Device:hasNaturalLight() then
                     UIManager:show(SpinWidget:new{
                         title_text = text,
+                        info_text = _("Enter percentage of warmth."),
                         value = self.warmth[num],
                         value_min = 0,
                         value_max = 100,
diff --git a/plugins/autowarmth.koplugin/suntime.lua b/plugins/autowarmth.koplugin/suntime.lua
index 09aab9e62..a42c12836 100644
--- a/plugins/autowarmth.koplugin/suntime.lua
+++ b/plugins/autowarmth.koplugin/suntime.lua
@@ -4,20 +4,26 @@
 -- The current source code of this file can be found on https://github.com/zwim/suntime.
 
 --[[--
-Module to calculate ephemeris and other times depending on the sun position
+Module to calculate ephemeris and other times depending on the sun position.
 
-Maximal errors from 2020-2050 are:
+Maximal errors from 2020-2050 (compared to https://midcdmz.nrel.gov/spa/) are:
 
+* -43.52° Christchurch 66s
+* -20.16° Mauritius:   25s
+* 20.30° Honolulu:     47s
 * 33.58° Casablanca:   24s
-* 37.97° Athene:       25s
-* 41.91° Rome:         28s
-* 47.25° Innsbruck:    14s
-* 52.32° Berlin:       32s
-* 64.14° Reykjavik:   113s
-* 65.69° Akureyri:   <110s (except *)
+* 35.68° Tokio:        50s
+* 37.97° Athene:       24s
+* 38°    Sacramento:   67s
+* 41.91° Rome:         27s
+* 47.25° Innsbruck:    13s
+* 52.32° Berlin:       30s
+* 59.92° Oslo:         42s
+* 64.14° Reykjavik:    66s
+* 65.69° Akureyri:    <24s (except *)
 * 70.67° Hammerfest: <105s (except **)
 
-*) A few days around beginning of summer (error <530s)
+*) A few days around beginning of summer (error <290s)
 
 **) A few days after and befor midnight sun (error <1200s)
 
@@ -41,9 +47,10 @@ Maximal errors from 2020-2050 are:
 --]]--
 
 -- math abbrevations
-local toRad = math.pi/180
-local toDeg = 1/toRad
+local pi = math.pi
+local pi_2 = pi/2
 
+local abs = math.abs
 local floor = math.floor
 local sin = math.sin
 local cos = math.cos
@@ -52,20 +59,38 @@ local asin = math.asin
 local acos = math.acos
 local atan = math.atan
 
+local toRad = pi/180
+local toDeg = 1/toRad
+
 local function Rad(x)
     return x*toRad
 end
 
+--------------------------------------------
+local speed_of_light = 2.99792E8
+local sun_radius = 6.96342e8
+local average_earth_radius = 6371e3
+local semimajor_axis = 149598022.96E3 -- earth orbit's major semi-axis in meter
+local average_speed_earth = 29.7859e3
+local aberration = asin(average_speed_earth/speed_of_light) -- Aberration relativistic
+local average_speed_equator = (2*pi * average_earth_radius) / (24*3600)
 --------------------------------------------
 
-local SunTime = {}
+ -- minimal twillight times in hours
+local min_civil_twilight = 20/60
+local min_nautic_twilight = 45/60 - min_civil_twilight
+local min_astronomic_twilight = 20/60 - min_nautic_twilight
 
-SunTime.astronomic = Rad(-18)
-SunTime.nautic =  Rad(-12)
-SunTime.civil = Rad(-6)
--- SunTime.eod = Rad(-49/60) -- approx. end of day
-SunTime.earth_flatten = 1 / 298.257223563 -- WGS84
-SunTime.average_temperature = 10 -- °C
+local SunTime = {
+        astronomic = Rad(-18),
+        nautic =  Rad(-12),
+        civil = Rad(-6),
+        -- eod = Rad(-49/60), -- approx. end of day
+        earth_flatten = 1 / 298.257223563, -- WGS84
+        average_temperature = 10, -- °C
+
+        times = {},
+    }
 
 ----------------------------------------------------------------
 
@@ -90,13 +115,12 @@ function SunTime:getZglAdvanced()
     local e4 = e3*e
     local e5 = e4*e
 
-    local M = self.M
     -- https://de.wikibooks.org/wiki/Astronomische_Berechnungen_f%C3%BCr_Amateure/_Himmelsmechanik/_Sonne
-    local C = (2*e - e3/4 + 5/96*e5) * sin(M)
-            + (5/4*e2 + 11/24*e4) * sin(2*M)
-            + (13/12*e3 - 43/64*e5) * sin(3*M)
-            + 103/96*e4 * sin(4*M)
-            + 1097/960*e5 * sin(5*M) -- rad
+    local C = (2*e - e3/4 + 5/96*e5) * self.sin_M
+            + (5/4*e2 + 11/24*e4) * self.sin_2M
+            + (13/12*e3 - 43/64*e5) * self.sin_3M
+            + 103/96*e4 * self.sin_4M
+            + 1097/960*e5 * self.sin_5M -- rad
 
     local lamb = self.L + C
     local tanL = tan(self.L)
@@ -109,6 +133,8 @@ end
 
 -- set current date or year/month/day daylightsaving hh/mm/ss
 -- if dst == nil use curent daylight saving of the system
+local days_in_month = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}
+
 function SunTime:setDate(year, month, day, dst, hour, min, sec)
     self.date = os.date("*t") -- get current day
 
@@ -116,11 +142,11 @@ function SunTime:setDate(year, month, day, dst, hour, min, sec)
         self.date.year = year
         self.date.month = month
         self.date.day = day
-        local feb = 28
         if year % 4 == 0 and (year % 100 ~= 0 or year % 400 == 0) then
-            feb = 29
+            days_in_month[2] = 29
+        else
+            days_in_month[2] = 28
         end
-        local days_in_month = {31, feb, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}
         self.date.yday = day
         for i = 1, month-1 do
             self.date.yday = self.date.yday + days_in_month[i]
@@ -132,10 +158,6 @@ function SunTime:setDate(year, month, day, dst, hour, min, sec)
             self.date.isdst = dst
         end
     end
-
-    if not self.getZgl then
-        self.getZgl = self.getZglAdvanced
-    end
 end
 
 --[[--
@@ -157,22 +179,27 @@ function SunTime:setPosition(name, latitude, longitude, time_zone, altitude, deg
 
     -- check for sane values
     -- latitudes are from -90° to +90°
-    if latitude > math.pi/2 then
-        latitude = math.pi/2
-    elseif latitude < -math.pi/2 then
-        latitude = -math.pi/2
+    if latitude > pi_2 then
+        latitude = pi_2
+    elseif latitude < -pi_2 then
+        latitude = -pi_2
     end
     -- longitudes are from -180° to +180°
-    if longitude > math.pi then
-        longitude = math.pi
-    elseif longitude < -math.pi then
-        longitude = -math.pi
+    if longitude > pi then
+        longitude = pi
+    elseif longitude < -pi then
+        longitude = -pi
     end
 
-    self.pos = {name, latitude = latitude, longitude = longitude, altitude = altitude}
+    latitude = atan((1-self.earth_flatten)^2 * tan(latitude))
+
+    self.pos = {name = name, latitude = latitude, longitude = longitude, altitude = altitude}
     self.time_zone = time_zone
 --    self.refract = Rad(36.35/60 * .5 ^ (altitude / 5538)) -- constant temperature
     self.refract = Rad(36.20/60 * (1 - 0.0065*altitude/(273.15+self.average_temperature)) ^ 5.255 )
+
+    self.sin_latitude = sin(self.pos.latitude)
+    self.cos_latitude = cos(self.pos.latitude)
 end
 
 --[[--
@@ -188,6 +215,10 @@ function SunTime:setAdvanced()
     self.getZgl = self.getZglAdvanced
 end
 
+--[[--
+  Function to get the equation of time, can be set by setSimple() or setAdvanced()
+--]]--
+SunTime.getZgl = SunTime.getZglAdvanced
 
 function SunTime:daysSince2000(hour)
     local delta = self.date.year - 2000
@@ -205,6 +236,8 @@ function SunTime:initVars(hour)
         hour = 12
     end
 
+    local after_noon = hour > 12
+
     local T = self:daysSince2000(hour)/36525 -- in Julian centuries form 2000-01-01 12:00
 
     --    self.num_ex = 0.0167086342 - 0.000042 * T
@@ -239,7 +272,6 @@ function SunTime:initVars(hour)
             + nT*(-2.04411E-2 -  nT* 0.00523E-3)))/3600 --°
     self.L = (L - floor(L/360)*360) * toRad
 
-
     -- see Numerical expressions for precession formulae ...
     -- Time is here in Julian centuries
     local omega =     102.93734808 + nT*(11612.35290e-1
@@ -250,6 +282,21 @@ function SunTime:initVars(hour)
     local M = L - omega --°
     self.M = (M - floor(M/360)*360) * toRad
 
+    self.sin_M = sin(self.M)
+    self.cos_M = cos(self.M)
+
+    -- sin(2x)=2 sin(x) cos(x)
+    self.sin_2M = 2 * self.sin_M * self.cos_M
+
+    -- sin(3x) = 3 sin(x) − 4 sin(x)^3
+    self.sin_3M = 3 * self.sin_M - 4 * self.sin_M^3
+
+    -- sin(4x) = 8 sin(x) cos(x)^3 - 4 sin(x) cos(x)
+    self.sin_4M = 8 * self.sin_M * self.cos_M^3 - 4 * self.sin_M * self.cos_M
+
+    -- sin(5x) = 5 sin(x) - 20 sin(x)^3+ 16 sin(x)^5
+    self.sin_5M = 5 * self.sin_M - 20 * self.sin_M^3 + 16 * self.sin_M^5
+
     -- Deklination nach astronomie.info
     --  local decl = 0.4095 * sin(0.016906 * (self.date.yday - 80.086))
     --Deklination nach Brodbeck (2001)
@@ -267,7 +314,7 @@ function SunTime:initVars(hour)
     --self.decl = asin(sin(ep)*sin(l))
 
     -- Deklination WMO-No.8 page I-7-37
-    local l =  self.L + math.pi + (1.915 * sin (self.M) + 0.020 * sin (2*self.M))*toRad
+    local l =  self.L + pi + (1.915 * self.sin_M + 0.020 * self.sin_2M)*toRad
     self.decl = asin(sin(self.epsilon)*sin(l))
 
     -- Nutation see https://de.wikipedia.org/wiki/Nutation_(Astronomie)
@@ -288,105 +335,160 @@ function SunTime:initVars(hour)
     self.decl = self.decl + delta_epsilon/3600*toRad
 
     -- https://de.wikipedia.org/wiki/Kepler-Gleichung#Wahre_Anomalie
-    self.E = self.M + self.num_ex * sin(self.M) + self.num_ex^2 / 2 * sin(2*self.M)
-    self.a = 149598022.96E3 -- große Halbachse in meter
-    self.r = self.a * (1 - self.num_ex * cos(self.E))
+    self.E = self.M + self.num_ex * self.sin_M + self.num_ex^2 / 2 * self.sin_2M
+    self.r = semimajor_axis * (1 - self.num_ex * cos(self.E))
+
+    --    self.eod = -atan(sun_radius/self.r) - self.refract
+    --                                            ^- astronomical refraction (at altitude)
 
-    self.eod = -atan(6.96342e8/self.r) - self.refract
-    --                ^--sun radius            ^- astronomical refraction (at altitude)
+    if after_noon then
+        self.eod = -atan((sun_radius-average_earth_radius*self.cos_latitude)/self.r) - self.refract
+        self.eod = self.eod + aberration
+    else
+        self.eod = -atan((sun_radius+average_earth_radius*self.cos_latitude)/self.r) - self.refract
+        self.eod = self.eod - aberration
+    end
 
     self.zgl = self:getZgl()
 end
 
 function SunTime:getTimeDiff(height)
-    local val = (sin(height) - sin(self.pos.latitude)*sin(self.decl))
-                / (cos(self.pos.latitude)*cos(self.decl))
+    local val = (sin(height) - self.sin_latitude*sin(self.decl))
+                / (self.cos_latitude*cos(self.decl))
 
-    if math.abs(val) > 1 then
+    if abs(val) > 1 then
         return
     end
-    return 12/math.pi * acos(val)
+    return 12/pi * acos(val)
+end
+
+-- get the sun height for a given time
+-- eod for considering sun diameter and astronomic refraction
+function SunTime:getHeight(time, eod)
+    time = time - 12 -- subtrace 12, because JD starts at 12:00
+    local val = cos(self.decl)*self.cos_latitude*cos(pi/12*time)
+        + sin(self.decl)*self.sin_latitude
+
+    if abs(val) > 1 then
+        return
+    end
+
+    if eod then
+        return asin(val) - eod -- todo self.eod is a bit to small
+    else
+        return asin(val)
+    end
 end
 
 -- Get time for a certain height
 -- Set hour near to expected time
--- Sed after_noon to true, if sunset is wanted
+-- Set after_noon to true, if sunset is wanted
+-- Set no_correct_dst if no daylight saving correction is wanted
 -- Result rise or set time
 --        nil sun does not reach the height
-function SunTime:calculateTime(height, hour, after_noon)
-    local dst = self.date.isdst and 1 or 0
-    local timeDiff = self:getTimeDiff(height, hour)
+function SunTime:calculateTime(height, hour, after_noon, no_correct_dst)
+    if not no_correct_dst then
+        if self.date.isdst and hour then
+            hour = hour - 1
+        end
+    end
+    self:initVars(hour) -- calculate self.eod
+    local timeDiff = self:getTimeDiff(height or self.eod, hour)
     if not timeDiff then
         return
     end
 
-    local local_correction = self.time_zone - self.pos.longitude*12/math.pi + dst - self.zgl
+    local local_correction = self.time_zone - self.pos.longitude*12/pi - self.zgl
     if not after_noon then
-        return 12 - timeDiff + local_correction
+        hour = 12 - timeDiff + local_correction
     else
-        return 12 + timeDiff + local_correction
+        hour = 12 + timeDiff + local_correction
+    end
+    if not no_correct_dst then
+        if self.date.isdst and hour then
+            hour = hour + 1
+        end
     end
+    return hour
 end
 
+-- Calculates the hour, when the sun reaches height
 -- If height is nil, use newly calculated self.eod
-function SunTime:calculateTimeIter(height, hour)
-    local after_noon = hour > 12
+-- hour gives a start value, default is used when hour == nil
+function SunTime:calculateTimeIter(height, hour, default_hour)
+    local after_noon = (hour and hour > 12) or (default_hour and default_hour > 12)
 
-    self:initVars(hour) -- calculate self.eod
-    hour = self:calculateTime(height or self.eod, hour, after_noon)
-    if hour ~= nil then
-        self:initVars(hour)  -- calculate self.eod
+    if not hour then -- do the iteration with the default value
+        hour = self:calculateTime(height or self.eod, default_hour, after_noon, true)
+    elseif hour and not default_hour then -- do the full iteration with value
+        hour = self:calculateTime(height or self.eod, hour, after_noon, true)
+    end -- if hour and default_hour are given don't do the first step
+
+    if hour ~= nil then -- do the last calculation step
         hour = self:calculateTime(height or self.eod, hour, after_noon)
     end
     return hour
 end
 
-function SunTime:calculateNoon()
-    self:initVars(12)
+function SunTime:calculateNoon(hour)
+    hour = hour or 12
+    self:initVars(hour)
+    local aberration_time = aberration / pi * 12 -- aberration in hours (angle/(2pi)*24)
+    local dst = self.date.isdst and 1 or 0
+    local local_correction = self.time_zone - self.pos.longitude*12/pi + dst - self.zgl
     if self.pos.latitude >= 0 then -- northern hemisphere
-        if math.pi/2 - self.pos.latitude + self.decl > self.eod then
-            local dst = self.date.isdst and 1 or 0
-            local local_correction = self.time_zone - self.pos.longitude*12/math.pi + dst - self.zgl
-            return 12 + local_correction
+        if pi_2 - self.pos.latitude + self.decl > self.eod then
+            if self:getHeight(hour) > 0 then
+                return hour + local_correction + aberration_time
+            end
         end
     else -- sourthern hemisphere
-        if math.pi/2 + self.pos.latitude - self.decl > self.eod then
-            local dst = self.date.isdst and 1 or 0
-            local local_correction = self.time_zone - self.pos.longitude*12/math.pi + dst - self.zgl
-            return 12 + local_correction
+        if pi_2 + self.pos.latitude - self.decl > self.eod then
+            if self:getHeight(hour) > 0 then
+                return hour + local_correction + aberration_time
+            end
         end
     end
 end
 
-function SunTime:calculateMidnight()
-    -- 24 is the midnight at the end of the current day,
+function SunTime:calculateMidnight(hour)
+    -- hour:
     -- 00 would be the beginning of the day
-    self:initVars(24)
+    -- 24 is the midnight at the end of the current day,
+    hour = hour or 24
+    self:initVars(hour)
+    local dst = self.date.isdst and 1 or 0
+    -- no aberration correction here, as you can't see the sun on her nadir ;-)
+    local local_correction = self.time_zone - self.pos.longitude*12/pi + dst - self.zgl
     if self.pos.latitude >= 0 then -- northern hemisphere
-        if math.pi/2 - self.pos.latitude - self.decl > self.eod then
-            local dst = self.date.isdst and 1 or 0
-            local local_correction = self.time_zone - self.pos.longitude*12/math.pi + dst - self.zgl
-            return 24 + local_correction
+        if pi_2 - self.pos.latitude - self.decl > self.eod then
+            if self:getHeight(hour) < 0 then
+                return hour + local_correction
+            end
         end
-    else -- southern hemisphere
-        if math.pi/2 + self.pos.latitude + self.decl > self.eod then
-            local dst = self.date.isdst and 1 or 0
-            local local_correction = self.time_zone - self.pos.longitude*12/math.pi + dst - self.zgl
-            return 24 + local_correction
+    else -- sourthern hemisphere
+        if pi_2 + self.pos.latitude + self.decl > self.eod then
+            if self:getHeight(hour) < 0 then
+                return hour + local_correction
+            end
         end
     end
-
 end
 
 --[[--
 Calculates the ephemeris and twilight times
 
+@param exact_twilight If not nil, then exact twilight times will be calculated.
+
 @usage
-SunTime:calculateTime()
+SunTime:calculateTimes(exact_twilight)
+
 
 Times are in hours or `nil` if not applicable.
 
 You can then access:
+    self.midnight_beginning
+
     self.rise_astronomic
     self.rise_nautic
     self.rise_civil
@@ -402,7 +504,7 @@ You can then access:
     self.midnight
 
 Or as values in a table:
-    self.times[1]  midnight - 24h
+    self.times[1]  midnight_beginning
     self.times[2]  rise_astronomic
     self.times[3]  rise_nautic
     self.times[4]  rise_civil
@@ -414,9 +516,12 @@ Or as values in a table:
     self.times[10] set_astronomic
     self.times[11] midnight
 --]]--
-function SunTime:calculateTimes()
+function SunTime:calculateTimes(fast_twilight)
     -- All or some the times can be nil at great latitudes
-    -- but either noon or midnight is not nil!
+    -- but either noon or midnight is not nil
+
+if not fast_twilight then
+    -- The canonical way is to calculate everything from scratch
     self.rise = self:calculateTimeIter(nil, 6)
     self.set = self:calculateTimeIter(nil, 18)
 
@@ -426,12 +531,43 @@ function SunTime:calculateTimes()
     self.set_nautic = self:calculateTimeIter(self.nautic, 18)
     self.rise_astronomic = self:calculateTimeIter(self.astronomic, 6)
     self.set_astronomic = self:calculateTimeIter(self.astronomic, 18)
+else
+    -- Calculate rise and set from scratch, use these values for twilight times
+    self.rise = self:calculateTimeIter(nil, 6)
+    self.rise_civil = self:calculateTimeIter(self.civil, self.rise - min_civil_twilight, 6)
+    self.rise_nautic = self:calculateTimeIter(self.nautic, self.rise_civil - min_nautic_twilight, 6)
+    self.rise_astronomic = self:calculateTimeIter(self.astronomic, self.rise_nautic - min_astronomic_twilight, 6)
 
+    self.set = self:calculateTimeIter(nil, 18)
+    self.set_civil = self:calculateTimeIter(self.civil, self.set + min_civil_twilight, 18)
+    self.set_nautic = self:calculateTimeIter(self.nautic, self.set_civil + min_nautic_twilight, 18)
+    self.set_astronomic = self:calculateTimeIter(self.astronomic, self.set_nautic + min_astronomic_twilight, 18)
+end
+
+    self.midnight_beginning = self:calculateMidnight(0)
     self.noon = self:calculateNoon()
     self.midnight = self:calculateMidnight()
 
-    self.times = {}
-    self.times[1]  = self.midnight and (self.midnight - 24)
+    -- Sometimes at high latitudes noon or midnight does not get calculated.
+    -- Maybe there is a minor bug in the calculateNoon/calculateMidnight functions.
+    if self.rise and self.set then
+        if not self.noon then
+            self.noon = (self.rise + self.set) / 2
+        end
+        if not self.midnight then
+            self.midnight = self.noon + 12
+        end
+        if not self.midnight_beginning then
+            self.midnight_beginning = self.midnight - 24
+        end
+    elseif self.rise and not self.set then -- only sunrise on that day
+        self.midnight = nil
+        self.midnight_beginning = nil
+    elseif self.set and not self.rise then -- only sunset on that day
+        self.noon = nil
+    end
+
+    self.times[1]  = self.midnight_beginning
     self.times[2]  = self.rise_astronomic
     self.times[3]  = self.rise_nautic
     self.times[4]  = self.rise_civil