Module:Lib astro

From Touhou Wiki
Jump to navigation Jump to search

Documentation for this module may be created at Module:Lib astro/doc

-- Helper library for sun and moon calculations
-- It's based on Chinese calendar calculations, but some of the functions were modified to calculate the Youkai calendar dates
-- References:
 -- "Calendrical Calculations" third edition by N. Dershowitz and E. M. Reingold
 -- "Astronomical Algorithms" errata 2 and 3 by J. Meeus

-- additional moon calculation tables
local common = require("Module:Common")
local moon = require("Module:Lib astro moon")

-- define global objects "astro" and "time"
if astro == nil then astro = {} end
if time == nil then time = {} end

-- == constants ==
-- general astrological constants
astro.J2000 = 2451545.0
astro.JulianCentury = 36525.0
astro.JulianMillenium = (astro.JulianCentury * 10)
astro.AstronomicalUnit = 149597870.0
astro.TropicalYear = 365.24219878
astro.MeanSynodicMonth = 29.530588853

-- solar longitudes for equinoxes and solstices
astro.spring = 0
astro.summer = 90
astro.autumn = 180
astro.winter = 270

-- time constants
time.wdays = {'Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday'}
time.GREGORIAN_EPOCH = 1721425.5
time.GENSOKYO_EPOCH = 2409676.5 -- to be calculated

-- for ease of use, we define some of the constants as local
local J2000 = astro.J2000
local JulianCentury = astro.JulianCentury
local JulianMillenium = astro.JulianMillenium
local AstronomicalUnit = astro.AstronomicalUnit
local TropicalYear = astro.TropicalYear
local MeanSynodicMonth = astro.MeanSynodicMonth

local GREGORIAN_EPOCH = time.GREGORIAN_EPOCH
local GENSOKYO_EPOCH = time.GENSOKYO_EPOCH

-- == local helper functions ==

-- Standard Japanese location for time calculation
-- We treat the location as function, since sometimes the reference point changes in time
-- For constant locations function needs only to return one table, regardless of time
-- {['lat'] = latitude, ['lon'] = longitude, ['dt'] = time offset from solar time (in hours), ['tz'] = timezone (in hours)}
-- lat and lon are positive for North and East
function time.location(t)
  -- TODO: should Mt. Fuji be used as reference point?
  -- 35°21′28.8″N 138°43′51.6″E

  local ret = {}
  -- The location changed starting with Gregorian year 1888,
  -- but since it's after Gensokyo was sealed off, we use only the old location
  -- if t < 2410637.5 then
    -- Before 1888
    ret = {['lat'] = 35.7, ['lon'] = 139.46, ['dt'] = 24/60, ['tz'] = (9 + 143 / 450)}
  -- else
    -- 1888 and later
    -- ret = {['lat'] = 35, ['lon'] = 135, ['dt'] = 0, ['tz'] = 9}
  -- end

  return ret
end

-- degrees to radians
local function d2r(d)
  return d * math.pi / 180.0
end

-- radians to degrees
local function r2d(r)
  return r * 180.0 / math.pi
end

-- reduce angle in degrees
local function fixangle(a)
  return a - 360.0 * math.floor(a / 360.0)
end

-- reduce angle in radians
local function fixangr(a)
  return a - 2 * math.pi * math.floor(a / (2 * math.pi))
end

-- sin of angle in degrees
local function dsin(d)
  return math.sin(d2r(d))
end

-- cos of angle in degrees
local function dcos(d)
  return math.cos(d2r(d))
end

-- mod - modulus for non-integers (in theory that's how lua's % operator works, but gonna define it explicitly)
local function mod(a, b)
  return a - b * math.floor(a / b)
end

-- amod - modulus function which returns numerator if modulus is zero
local function amod(a, b)
  return mod(a - 1, b) + 1
end

-- round - round to nearest
local function round(x)
  return math.floor(x + 0.5)
end

-- fulld - round to nearest full day
local function fulld(jd)
  return math.floor(jd - 0.5) + 0.5
end

-- == Gregorian time conversion functions ==

-- convert julian time to hour, minute, second
function time.jhms(jd)
  jd = jd + 0.5 -- astronomical to civil
  ij = ((jd - math.floor(jd)) * 86400.0) + 0.5
  return {math.floor(ij / 3600), math.floor((ij / 60) % 60), math.floor(ij % 60)}
end

-- calculate the day of week from julian day
function time.jwday(jd)
  -- +0.5 for 0 as Monday, +1.5 for 0 as Sunday
  return mod(math.floor(jd + 0.5), 7)
end

function time.leap_gregorian(year)
  return ((year % 4) == 0) and (not (((year % 100) == 0) and ((year % 400) ~= 0)))
end

function time.gregorian_to_jd(year, month, day, hour)
  return (GREGORIAN_EPOCH - 1) + 
    (365 * (year - 1)) +
    math.floor((year - 1) / 4) + 
    (-math.floor((year - 1) / 100)) +
    math.floor((year - 1) / 400) +
    math.floor((((367 * month) - 362) / 12) + common.cv((month <= 2), 0, common.cv(time.leap_gregorian(year), -1, -2)) + day) + hour / 24
end

function time.jd_to_gregorian(jd)
  local wjd, depoch, quadricent, qdc, cent, dcent, quad, dquad, yindex, dyindex, year, yearday, leapadj
  
  wjd = fulld(jd)
  depoch = wjd - GREGORIAN_EPOCH
  quadricent = math.floor(depoch / 146097)
  qdc = depoch % 146097
  cent = math.floor(qdc / 36524)
  dcent = qdc % 36524
  quad = math.floor(dcent / 1461)
  dquad = dcent % 1461
  yindex = math.floor(dquad / 365)
  year = (quadricent * 400) + (cent * 100) + (quad * 4) + yindex
  if not ((cent == 4) or (yindex == 4)) then
    year = year + 1
  end
  
  yearday = wjd - time.gregorian_to_jd(year, 1, 1, 0)
  leapadj = common.cv((wjd < time.gregorian_to_jd(year, 3, 1, 0)), 0, common.cv(time.leap_gregorian(year), 1, 2))
  
  month = math.floor((((yearday + leapadj) * 12) + 373) / 367)
  day = (wjd - time.gregorian_to_jd(year, month, 1, 0)) + 1
  
  return year,month,day
end

-- == General time functions ==

-- 13.7
function time.zone_from_longitude(l)
  return l / 360
end

-- 13.8
function time.universal_from_local(jd, loc)
  return jd - zone_from_longitude(loc['lon'])
end

-- 13.9
function time.local_from_universal(jd, loc)
  return jd + zone_from_longitude(loc['lon'])
end

-- 13.10
function time.standard_from_universal(jd, loc)
  return jd + loc['tz'] / 24
end

-- 13.11
function time.universal_from_standard(jd, loc)
  return jd - loc['tz'] / 24
end

-- 13.14
function astro.ephemeris_correction(jd)
  local year,month,day = time.jd_to_gregorian(fulld(jd))
  local c = (time.gregorian_to_jd(year, 7, 1, 0) - time.gregorian_to_jd(1900, 1, 1, 0)) / JulianCentury
  local c2 = c * c
  local c3 = c * c2
  local c4 = c * c3
  local c5 = c * c4
  local c6 = c * c5
  local c7 = c * c6
  local c8 = c * c7
  local c9 = c * c8
  local c10= c * c9
  local ec

      if (1988 <= year) and (year <= 2019) then
    ec = (year - 1933) / 86400
  elseif (1900 <= year) and (year <= 1987) then
    ec = -0.00002 + 0.000297 * c + 0.025184 * c2 -
          0.181133 * c3 + 0.553040 * c4 -
          0.861938 * c5 + 0.677066 * c6 -
          0.212591 * c7
  elseif (1800 <= year) and (year <= 1899) then
    ec = -0.00009 + 0.003844 * c + 0.083563 * c2 +
          0.865736 * c3 + 4.867575 * c4 +15.845535 * c5 +
         31.332267 * c6 +38.291999 * c7 +28.316289 * c8 +
         11.636204 * c9 + 2.043794 * c10
  elseif (1700 <= year) and (year <= 1799) then
    local y17 = year - 1700
    ec = (8.118780842 - 0.005092142 * y17 + 0.003336121 * y17 * y17 - 0.0000266484 * y17 * y17 * y17) / 86400
  elseif (1620 <= year) and (year <= 1699) then
    local y16 = year - 1600
    ec = (196.58333 - 4.0675 * y16 + 0.0219167 * y16 * y16) / 86400
  else
    local x = (time.gregorian_to_jd(year, 1, 1, 0) - time.gregorian_to_jd(1810, 1, 1, 0)) + 0.5
    ec = (x*x / 41048480 - 15) / 86400
  end

  return ec
end

-- 13.15
function time.dynamical_from_universal(t)
  return t + astro.ephemeris_correction(t)
end

-- 13.16
function time.universal_from_dynamical(t)
  return t - astro.ephemeris_correction(t)
end

-- == main astrological functions ==

-- get the sun's longitude at given Julian Date (JD)
function astro.solar_longitude(jd)
  local L0, M, C
  local T = (jd - J2000) / JulianCentury
  local T2 = T * T

  L0 = fixangle(280.46646 + (36000.76983 * T) + ( 0.0003032 * T2))
  M  = fixangle(357.52911 + (35999.05029 * T) + (-0.0001537 * T2))
  C  = ((1.914602 + (-0.004817 * T) + (-0.000014 * T2)) * dsin(M)) +
       ((0.019993 - ( 0.000101 * T)) * dsin(2 * M)) +
       (0.000289 * dsin(3 * M))
  return L0 + C
end

-- get the date after jd, at which sun is at given position l = [0,360)
-- own adaptation of 13.33
function astro.solar_longitude_after(l, jd)
  local dl, t0
  local rate = TropicalYear / 360

  t0 = fulld(jd)

  -- first approximation of time
  t0 = t0 + rate * mod((l - astro.solar_longitude(t0)), 360)

  -- determine the search range
  local a = math.max(jd, t0 - 5)
  local b = t0 + 5
  local lp = 0

  -- initial condition
  t0 = math.max(a, t0)

  -- repeat until we get 0.1 deg. accuracy (~2.4h), fall out of <a,b> range or get more than 10 loops
  repeat
    lp = lp + 1
    dl = mod(l - astro.solar_longitude(t0) + 180, 360) - 180
    t0 = t0 + rate * dl
  until ((math.abs(dl) < 0.1) or (t0 < a) or (t0 > b) or (lp > 10))

  return t0,a,b,lp
end

-- 13.41
function astro.estimate_prior_solar_longitude(l, jd)
  local rate = TropicalYear / 360
  local t, d

  t = jd - rate*mod(astro.solar_longitude(jd) - l, 360)
  d = mod(astro.solar_longitude(t) - l + 180, 360) - 180

  return math.min(jd, t - rate * d)
end

-- 13.44
function astro.nth_new_moon(n)
  local n0 = 24724
  local k = n - n0
  local c = k / 1236.85
  local c2 = c * c
  local c3 = c * c2
  local c4 = c * c3

  local approx        = J2000 + 5.09765 + MeanSynodicMonth * 1236.85 * c + 0.0001337 * c2 - 0.000000150 * c3 + 0.00000000073 * c4
  local E             = 1 - 0.002516 * c - 0.0000074 * c2
  local solar_anomaly = 2.5534   + 1236.85      * 29.10535669 * c - 0.0000218 * c2 - 0.00000011 * c3
  local lunar_anomaly = 201.5643 + 385.81693528 * 1236.85     * c + 0.0107438 * c2 + 0.00001239 * c3 - 0.000000058 * c4
  local moon_argument = 160.7108 + 390.67050274 * 1236.85     * c - 0.0016341 * c2 - 0.00000227 * c3 + 0.000000011 * c4
  local omega         = 124.7746 - 1.56375580   * 1236.85     * c + 0.0020691 * c2 + 0.00000215 * c3

  local correction = -0.00017 * dsin(omega)
  local e = {1, E, E*E}
  for i=1,24 do
    correction = correction +
        moon.T133v[i] * e[moon.T133w[i] + 1] * dsin(moon.T133x[i] * solar_anomaly + moon.T133y[i] * lunar_anomaly + moon.T133z[i] * moon_argument)
  end

  local extra = 0.000325 * dsin(299.77 + 132.8475848 * c - 0.009173 * c2)

  local additional = 0
  for i=1,13 do
    additional = additional + moon.T134l[i] * dsin(moon.T134i[i] + moon.T134j[i] * k)
  end

  return time.universal_from_dynamical(approx + correction + extra + additional)
end

-- 13.47, based on Meeus
function astro.lunar_longitude(jd)
  local T, T2, T3, T4, LP, D, M, MP, F
  local A1, A2, E, E2, Sl
  local Eterm
  local i
 
  T = (jd - J2000) / JulianCentury
  T2 = T * T
  T3 = T * T2
  T4 = T * T3
 
  LP = 218.3164477 + 481267.88123421 * T - 0.0015786 * T2 + T3 /   538841.0 - T4 /  65194000.0
  D  = 297.8501921 + 445267.1114034  * T - 0.0018819 * T2 + T3 /   545868.0 - T4 / 113065000.0
  M  = 357.5291092 +  35999.0502909  * T - 0.0001536 * T2 + T3 / 24490000.0
  MP = 134.9633964 + 477198.8675055  * T + 0.0087414 * T2 + T3 /    69699.0 - T4 /  14712000.0
  F  =  93.2720950 + 483202.0175233  * T - 0.0036539 * T2 - T3 /  3526000.0 + T4 / 863310000.0
 
  A1 = 119.75 +    131.849 * T
  A2 =  53.09 + 479264.290 * T
 
  E  = 1 - 0.002516 * T - 0.0000074 * T2
  E2 = E * E
 
  Sl = 0.0
  for i=1,60 do
    Eterm = 1
    if math.abs(moon.T45AM[i]) == 1 then
      Eterm = E
    elseif math.abs(moon.T45AM[i]) == 2 then
      Eterm = E2
    end
 
    Sl = Sl + moon.T45AL[i] * Eterm * dsin(fixangle(moon.T45AD[i] * D + moon.T45AM[i] * M + moon.T45AMP[i] * MP + moon.T45AF[i] * F))
  end
 
  Sl = Sl + 3958 * dsin(fixangle(A1))   + 1962 * dsin(fixangle(LP-F))  + 318 * dsin(fixangle(A2))
 
  return mod(LP + Sl / 1000000.0, 360)
end

-- 13.53
function astro.lunar_phase(jd)
  local f = mod(astro.lunar_longitude(jd) - astro.solar_longitude(jd), 360)
  local t0 = astro.nth_new_moon(0)
  local n = round((jd-t0)/MeanSynodicMonth)
  local fp = 360 * mod((jd - astro.nth_new_moon(n)) / MeanSynodicMonth, 1)
  
  return common.cv(math.abs(f - fp) > 180, fp, f)
end

function astro.lunar_phase_after(l, jd)
  local dl, t0
  local rate = MeanSynodicMonth / 360
  
  t0 = fulld(jd)
  
  t0 = t0 + rate * mod((l - astro.lunar_phase(t0)), 360)
  
  local a = math.max(jd, t0 - 5)
  local b = t0 + 5
  local lp = 0
  
  t0 = math.max(a, t0)
  repeat
    lp = lp + 1
    dl = mod(l - astro.lunar_phase(t0) + 180, 360) - 180
    t0 = t0 + rate * dl
  until ((math.abs(dl) < 0.1) or (t0 < a) or (t0 > b) or (lp > 10))
  
  return t0,a,b,lp
end

-- 13.45
function astro.new_moon_before(jd)
  local t0 = astro.nth_new_moon(0)
  local f  = astro.lunar_phase(jd)
  local n  = round((jd - t0) / MeanSynodicMonth - f / 360)

  local k = n-1
  while astro.nth_new_moon(k+1) < jd do
    k = k + 1
  end

  return astro.nth_new_moon(k)
end

-- 13.46
function astro.new_moon_at_or_after(jd)
  local t0 = astro.nth_new_moon(0)
  local f  = astro.lunar_phase(jd)
  local n = round((jd - t0) / MeanSynodicMonth - f / 360)

  local k = n
  while astro.nth_new_moon(k) < jd do
    k = k + 1
  end

  return astro.nth_new_moon(k)
end

-- 17.1
function astro.current_major_solar_term(jd)
  local s = astro.solar_longitude(time.universal_from_standard(jd, time.location(jd)))

  return amod(2 + math.floor(s/30), 12) 
end

-- 17.3
function astro.gensokyo_solar_longitude_on_or_after(l, jd)
  local t = astro.solar_longitude_after(l, universal_from_standard(jd, time.location(jd)))
  return time.standard_from_universal(t, time.location(t))
end

-- 17.7
function time.midnight_in_gensokyo(jd)
  return time.universal_from_standard(jd, time.location(jd))
end

-- 17.4
function astro.major_solar_term_on_or_after(jd)
  local s = astro.solar_longitude(time.midnight_in_gensokyo(jd))
  local l = mod(30 * math.ceil(s / 30), 360)
  return astro.gensokyo_solar_longitude_on_or_after(l, jd)
end

-- 17.5
function astro.current_minor_solar_term(jd)
  local s = astro.solar_longitude(time.universal_from_standard(jd, time.location(jd)))

  return amod(3 + math.floor((s - 15) / 30), 12)
end

-- 17.6
function astro.minor_solar_term_on_or_after(jd)
  s = astro.solar_longitude(jd, time.midnight_in_gensokyo(jd))
  l = mod(30 * math.floor((s - 15) / 30) + 15, 360)
  return astro.gensokyo_solar_longitude_or_or_after(l, jd)
end

-- 17.8, originally winter_solstice_on_or_before
function astro.gensokyo_spring_equinox_on_or_before(jd)
  local approx = astro.estimate_prior_solar_longitude(astro.spring, time.midnight_in_gensokyo(jd + 1))
  -- local approx = astro.estimate_prior_solar_longitude(astro.winter, time.midnight_in_gensokyo(jd + 1))

  -- modified for looking for spring
  local day = fulld(approx) - 1
  while mod(astro.solar_longitude(time.midnight_in_gensokyo(day+1)) + 180, 360) <= (astro.spring + 180) do
  -- while mod(astro.solar_longitude(time.midnight_in_gensokyo(day+1)), 360) <= astro.winter do
    day = day + 1
  end

  return day
end

-- 17.9
function astro.gensokyo_new_moon_on_or_after(jd)
  local t = astro.new_moon_at_or_after(time.midnight_in_gensokyo(jd))
  return fulld(time.standard_from_universal(t, time.location(t)))
end

-- 17.10
function astro.gensokyo_new_moon_before(jd)
  local t = astro.new_moon_before(time.midnight_in_gensokyo(jd))
  return fulld(time.standard_from_universal(t, time.location(t)))
end

-- 17.11
function astro.gensokyo_no_major_solar_term(jd)
  return astro.current_major_solar_term(jd) == astro.current_major_solar_term(astro.gensokyo_new_moon_on_or_after(jd + 1))
end

-- 17.12
function astro.gensokyo_prior_leap_month(mp, m)
  return (m >= mp) and (astro.gensokyo_no_major_solar_term(m) or astro.gensokyo_prior_leap_month(mp, astro.gensokyo_new_moon_before(m)))
end

-- 17.13
function astro.gensokyo_new_year_in_sui(jd)
  local s1 = astro.gensokyo_spring_equinox_on_or_before(jd)
  local s2 = astro.gensokyo_spring_equinox_on_or_before(s1 + 370)
  local m12 = astro.gensokyo_new_moon_on_or_after(s1 + 1)
  local m13 = astro.gensokyo_new_moon_on_or_after(m12 + 1)
  local nm11 = astro.gensokyo_new_moon_before(s2 + 1)

  if (round((nm11 - m12) / MeanSynodicMonth) == 12) and (astro.gensokyo_no_major_solar_term(m12) or astro.gensokyo_no_major_solar_term(m13)) then
    return astro.gensokyo_new_moon_on_or_after(m13 + 1)
  else
    return m13
  end
end

-- 17.14
function astro.gensokyo_new_year_on_or_before(jd)
  local new_year = astro.gensokyo_new_year_in_sui(jd)
  if jd >= new_year then
    return new_year
  else
    return astro.gensokyo_new_year_in_sui(jd - 180)
  end
end

-- 17.16
function time.jd_to_gensokyo(jd)
  local s1   = astro.gensokyo_spring_equinox_on_or_before(jd)
  local s2   = astro.gensokyo_spring_equinox_on_or_before(s1 + 370)
  local m12  = astro.gensokyo_new_moon_on_or_after(s1 + 1)
  local nm11 = astro.gensokyo_new_moon_before(s2 + 1)
  local m    = astro.gensokyo_new_moon_before(jd + 1)

  local leap_year = (round((nm11 - m12) / MeanSynodicMonth) == 12)
  local month = amod(round((m - m12) / MeanSynodicMonth) - common.cv(leap_year and astro.gensokyo_prior_leap_month(m12, m), 1, 0), 12)

  local leap_month = leap_year and astro.gensokyo_no_major_solar_term(m) and not astro.gensokyo_prior_leap_month(m12, astro.gensokyo_new_moon_before(m))

  local season = math.floor(1.5 - month / 12 + (jd - GENSOKYO_EPOCH) / TropicalYear) - 1

  local day = jd - m + 1

  return season,month,leap_month,day
end

-- 17.17
function time.gensokyo_to_jd(season, month, leap_month, day)
  local mid_year = GENSOKYO_EPOCH + math.floor((season + 0.5) * TropicalYear)
  local new_year = astro.gensokyo_new_year_on_or_before(mid_year)
  local p        = astro.gensokyo_new_moon_on_or_after(new_year + (month-1) * 29)
  local ds,dm,dl,dd = time.jd_to_gensokyo(p)

  local prior_new_moon = p
  if ((month ~= dm) or (leap_month ~= dl)) then
    prior_new_moon = astro.gensokyo_new_moon_on_or_after(p+1)
  end

  return prior_new_moon + day - 1
end

function time.gensokyo_leap(season)
  local s1       = time.gensokyo_to_jd(season, 1, false, 1)
  local s2       = time.gensokyo_to_jd(season + 1, 1, false, 1)

  return (round((s2 - s1) / MeanSynodicMonth) == 13)
end

return {astro = astro, time = time}
-- [[Category:Lua Libraries|{{PAGENAME}}]]