2
0
mirror of https://github.com/koreader/koreader synced 2024-11-13 19:11:25 +00:00
koreader/plugins/autowarmth.koplugin/suntime.lua
2021-12-05 22:14:00 +01:00

597 lines
20 KiB
Lua
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

-- Author: Martin Zwicknagl (zwim)
-- Date: 2021-10-29
-- 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.
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
* 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 <290s)
**) A few days after and befor midnight sun (error <1200s)
@usage
local SunTime = require("suntime")
time_zone = 0
altitude = 50
degree = true
SunTime:setPosition("Reykjavik", 64.14381, -21.92626, timezone, altitude, degree)
SunTime:setAdvanced()
SunTime:setDate()
SunTime:calculateTimes()
print(SunTime.rise, SunTime.set, SunTime.set_civil) -- or similar see calculateTime()
@module suntime
--]]--
-- math abbrevations
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
local tan = math.tan
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)
--------------------------------------------
-- 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
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 = {},
}
----------------------------------------------------------------
-- simple 'Equation of time' good for dates between 2008-2027
-- errors for latitude 20° are within 1min
-- 47° are within 1min 30sec
-- 65° are within 5min
-- https://www.astronomie.info/zeitgleichung/#Auf-_und_Untergang (German)
function SunTime:getZglSimple()
local T = self.date.yday
return -0.171 * sin(0.0337 * T + 0.465) - 0.1299 * sin(0.01787 * T - 0.168)
end
-- more advanced 'Equation of time' good for dates between 1800-2200
-- errors are better than with the simple method
-- https://de.wikipedia.org/wiki/Zeitgleichung (German) and
-- more infos on http://www.hlmths.de/Scilab/Zeitgleichung.pdf (German)
function SunTime:getZglAdvanced()
local e = self.num_ex
local e2 = e*e
local e3 = e2*e
local e4 = e3*e
local e5 = e4*e
-- https://de.wikibooks.org/wiki/Astronomische_Berechnungen_f%C3%BCr_Amateure/_Himmelsmechanik/_Sonne
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)
local tanLamb = tan(lamb)
local cosEps = cos(self.epsilon)
local zgl = atan( (tanL - tanLamb*cosEps) / (1 + tanL*tanLamb*cosEps) ) --rad
return zgl*toDeg/15 -- to hours *4'/60
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
if year and month and day then
self.date.year = year
self.date.month = month
self.date.day = day
if year % 4 == 0 and (year % 100 ~= 0 or year % 400 == 0) then
days_in_month[2] = 29
else
days_in_month[2] = 28
end
self.date.yday = day
for i = 1, month-1 do
self.date.yday = self.date.yday + days_in_month[i]
end
self.date.hour = hour or 12
self.date.min = min or 0
self.date.sec = sec or 0
if dst ~= nil then
self.date.isdst = dst
end
end
end
--[[--
Set position for later calculations
@param name Name of the location
@param latitude Geographical latitude, North is positive
@param longitude Geographical longitude, West is negative
@param time_zone Timezone e.g. CET = +1
@param altitude Altitude of the location above the sea level
@param degree if `nil` latitude and longitue are in radian, else in decimal degree
--]]--
function SunTime:setPosition(name, latitude, longitude, time_zone, altitude, degree)
altitude = altitude or 200
if degree then
latitude = latitude * toRad
longitude = longitude * toRad
end
-- check for sane values
-- latitudes are from -90° to +90°
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 > pi then
longitude = pi
elseif longitude < -pi then
longitude = -pi
end
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
--[[--
Use a simple equation of time (valid for the years 2008-2027)
--]]--
function SunTime:setSimple()
self.getZgl = self.getZglSimple
end
--[[--
Use an advanced equation of time (valid for the years 1800-2200 at least)
--]]--
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
local leap = floor((delta-1)/4)
return 365 * delta + leap + self.date.yday + (hour-12)/24 -- WMO No.8, rebased for 2000-01-01 12:00
end
-- more accurate parameters of earth orbit from
-- Title: Numerical expressions for precession formulae and mean elements for the Moon and the planets
-- Authors: Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touze, M., Francou, G., & Laskar, J., ,
-- Journal: Astronomy and Astrophysics (ISSN 0004-6361), vol. 282, no. 2, p. 663-683
-- Bibliographic Code: 1994A&A...282..663S
function SunTime:initVars(hour)
if not hour then
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
-- numerical eccentricity of earth's orbit
-- see wikipedia: https://de.wikipedia.org/wiki/Erdbahn-> Meeus
-- and Numerical expressions for preccession formulae
-- time is in Julian centuries
self.num_ex = 0.0167086342 + T*(-0.0004203654e-1
+ T*(-0.0000126734e-2 + T*( 0.0000001444e-3
+ T*(-0.0000000002e-4 + T* 0.0000000003e-5))))
-- self.epsilon = (23 + 26/60 + 21/3600 - 46.82/3600 * T) * toRad
-- earth's obliquity to the ecliptic
-- see Numerical expressions for precession formulae ...
-- Time is here in Julian centuries
-- local epsilon = 23 + 26/60 + (21.412 + T*(-468.0927E-1
-- + T*(-0.0152E-2 + T*(1.9989E-3
-- + T*(-0.0051E-4 - T*0.0025E-5)))))/3600 --°
-- Astronomical Almanac 2010, p. B52
-- Time is here in Julian centuries
local epsilon = 23 + 26/60 + (21.406 - T*(46.836769
- T*( 0.0001831 + T*(0.00200340
+ T*(-5.76E-7 - T* 4.34E-8)))))/3600 --°
self.epsilon = epsilon * toRad
-- see Numerical expressions for precession formulae ...
-- mean longitude
local nT = T * (36000.7690/35999.3720) -- convert from equinox to date
local L = 100.46645683 + (nT*(1295977422.83429E-1
+ 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
+ nT*(53.27577e-2 + nT*(-0.14095e-3
+ nT*( 0.11440e-4 + nT* 0.00478e-5))))/3600 --°
-- mean anomaly
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)
-- local decl = 0.40954 * sin(0.0172 * (self.date.yday - 79.349740))
-- Deklination WMO-No.8 page I-7-37
--local T = self.days_since_2000
--local L = 280.460 + 0.9856474 * T
--L = (L - floor(L/360)*360) * toRad
--local g = 357.528 + 0.9856003 * T -- mean anomaly
--g = (g - floor(g/360)*360) * toRad
--local l = L + (1.915 * sin (g) + 0.020 * sin (2*g))*toRad
--local ep = self.epsilon
-- -- sin(decl) = sin(ep)*sin(l)
--self.decl = asin(sin(ep)*sin(l))
-- Deklination WMO-No.8 page I-7-37
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)
local A = { 2.18243920 - 33.7570460 * T,
-2.77624462 + 1256.66393 * T,
7.62068856 + 16799.4182 * T,
4.36487839 - 67.5140919 * T}
local B = {92025e-4 + 8.9e-4 * T,
5736e-4 - 3.1e-4 * T,
977e-4 - 0.5e-4 * T,
-895e-4 + 0.5e-4 * T}
local delta_epsilon = 0 --"
for i = 1, #A do
delta_epsilon = delta_epsilon + B[i]*cos(A[i])
end
-- add nutation to declination
self.decl = self.decl + delta_epsilon/3600*toRad
-- https://de.wikipedia.org/wiki/Kepler-Gleichung#Wahre_Anomalie
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)
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) - self.sin_latitude*sin(self.decl))
/ (self.cos_latitude*cos(self.decl))
if abs(val) > 1 then
return
end
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 -- self.eod might be a bit too small
else
return asin(val)
end
end
-- Get time for a certain height
-- Set hour near to expected time
-- 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, 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/pi - self.zgl
if not after_noon then
hour = 12 - timeDiff + local_correction
else
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
-- 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)
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(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 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 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(hour)
-- hour:
-- 00 would be the beginning of the day
-- 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 pi_2 - self.pos.latitude - self.decl > self.eod then
if self:getHeight(hour) < 0 then
return hour + local_correction
end
end
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: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
self.rise
self.noon
self.set
self.set_civil
self.set_nautic
self.set_astronomic
self.midnight
Or as values in a table:
self.times[1] midnight_beginning
self.times[2] rise_astronomic
self.times[3] rise_nautic
self.times[4] rise_civil
self.times[5] rise
self.times[6] noon
self.times[7] set
self.times[8] set_civil
self.times[9] set_nautic
self.times[10] set_astronomic
self.times[11] midnight
--]]--
function SunTime:calculateTimes(fast_twilight)
-- All or some the times can be nil at great latitudes
-- 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)
self.rise_civil = self:calculateTimeIter(self.civil, 6)
self.set_civil = self:calculateTimeIter(self.civil, 18)
self.rise_nautic = self:calculateTimeIter(self.nautic, 6)
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()
-- 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
self.times[5] = self.rise
self.times[6] = self.noon
self.times[7] = self.set
self.times[8] = self.set_civil
self.times[9] = self.set_nautic
self.times[10] = self.set_astronomic
self.times[11] = self.midnight
end
-- Get time in seconds (either actual time in hours or date struct)
function SunTime:getTimeInSec(val)
if not val then
val = os.date("*t")
end
if type(val) == "table" then
return val.hour*3600 + val.min*60 + val.sec
end
return val*3600
end
return SunTime