Skip to content

Instantly share code, notes, and snippets.

@ponapalt

ponapalt/qreki.py

Created October 11, 2024 07:17
Show Gist options
  • Save ponapalt/f4be2109ad7575e7e89e347e5f3e62b4 to your computer and use it in GitHub Desktop.
Save ponapalt/f4be2109ad7575e7e89e347e5f3e62b4 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
qreki.py
import math
import datetime
from typing import List
# 本スクリプトは高野 英明氏のQREKI.AWKをもとに、Claude 3.5 Sonnetを使用してpython用に変換したものです。
# GitHub Gistにはファイル単体しか置けませんので、配布規定通りオリジナルを同梱することができません。
# オリジナルは下記からダウンロードお願いします。
# https://www.vector.co.jp/soft/dl/dos/personal/se016093.html
# 以下は使用サンプルです。
# def main():
# if len(sys.argv) <= 3:
# date = datetime.date.today()
# else:
# date = datetime.date(int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3]))
#
# lunar_date = Qreki.from_gregorian(date)
# rokuyou = Qreki.get_rokuyou(date)
# is_leap_month = Qreki.get_leap_month(date)
#
# k_month = " 閏" if is_leap_month else ""
# k_month += "正月" if lunar_date[1] == 1 else f"{lunar_date[1]:2d}月"
#
# k_date = "朔日" if lunar_date[2] == 1 else f"{lunar_date[2]:2d}日"
#
# print(f"西暦{date.year:4d}年 {date.month:2d}月 {date.day:2d}日は、旧暦{lunar_date[0]:4d}年 {k_month} {k_date} {rokuyou}です。")
#
# if __name__ == "__main__":
# main()
class Qreki:
PI: float = 3.141592653589793238462
k: float = PI / 180.0
@staticmethod
def calc_kyureki(tm0: float) -> List[int]:
"""
旧暦を計算する
:param tm0: ユリウス日
:return: [年, 閏月フラグ, 月, 日]
"""
kyureki = [0] * 4
chu = [[0.0, 0.0] for _ in range(4)]
saku = [0.0] * 5
m = [[0, 0, 0] for _ in range(5)]
Qreki._before_nibun(tm0, chu)
for i in range(1, 4):
Qreki._calc_chu(chu[i-1][0] + 32.0, i, chu)
saku[0] = Qreki._calc_saku(chu[0][0])
for i in range(1, 5):
tm = saku[i-1] + 30.0
saku[i] = Qreki._calc_saku(tm)
if abs(int(saku[i-1]) - int(saku[i])) <= 26.0:
saku[i] = Qreki._calc_saku(saku[i-1] + 35.0)
if int(saku[1]) <= int(chu[0][0]):
saku = saku[1:] + [Qreki._calc_saku(saku[3] + 35.0)]
elif int(saku[0]) > int(chu[0][0]):
saku = [Qreki._calc_saku(saku[0] - 27.0)] + saku[:-1]
lap = 1 if int(saku[4]) <= int(chu[3][0]) else 0
m[0][0] = int(chu[0][1] / 30.0) + 2
if m[0][0] > 12:
m[0][0] -= 12
m[0][2] = int(saku[0])
m[0][1] = 0
for i in range(1, 5):
if lap == 1 and i != 1:
if int(chu[i-1][0]) <= int(saku[i-1]) or int(chu[i-1][0]) >= int(saku[i]):
m[i-1][0] = m[i-2][0]
m[i-1][1] = 1
m[i-1][2] = int(saku[i-1])
lap = 0
m[i][0] = m[i-1][0] + 1
if m[i][0] > 12:
m[i][0] -= 12
m[i][2] = int(saku[i])
m[i][1] = 0
state = 0
for i in range(5):
if int(tm0) < int(m[i][2]):
state = 1
break
elif int(tm0) == int(m[i][2]):
state = 2
break
if state == 0 or state == 1:
i -= 1
kyureki[1] = m[i][1]
kyureki[2] = m[i][0]
kyureki[3] = int(tm0) - int(m[i][2]) + 1
a = [0] * 6
Qreki._JD2YMDT(tm0, a)
kyureki[0] = a[0]
if kyureki[2] > 9 and kyureki[2] > a[1]:
kyureki[0] -= 1
return kyureki
@staticmethod
def _calc_chu(tm: float, i: int, chu: List[List[float]]) -> None:
"""中気の時刻を計算する"""
tm1 = int(tm)
tm2 = tm - tm1
tm2 -= 9.0 / 24.0
t = (tm2 + 0.5) / 36525.0
t += (tm1 - 2451545.0) / 36525.0
rm_sun = Qreki._LONGITUDE_SUN(t)
rm_sun0 = 30.0 * int(rm_sun / 30.0)
delta_t1, delta_t2 = 1.0, 1.0
while abs(delta_t1 + delta_t2) > (1.0 / 86400.0):
t = (tm2 + 0.5) / 36525.0
t += (tm1 - 2451545.0) / 36525.0
rm_sun = Qreki._LONGITUDE_SUN(t)
delta_rm = rm_sun - rm_sun0
if delta_rm > 180.0:
delta_rm -= 360.0
elif delta_rm < -180.0:
delta_rm += 360.0
delta_t1 = int(delta_rm * 365.2 / 360.0)
delta_t2 = delta_rm * 365.2 / 360.0 - delta_t1
tm1 = tm1 - delta_t1
tm2 = tm2 - delta_t2
if tm2 < 0:
tm2 += 1.0
tm1 -= 1.0
chu[i][0] = tm2 + 9.0 / 24.0 + tm1
chu[i][1] = rm_sun0
@staticmethod
def _before_nibun(tm: float, nibun: List[List[float]]) -> None:
"""直前の二分二至の時刻を求める"""
tm1 = int(tm)
tm2 = tm - tm1
tm2 -= 9.0 / 24.0
t = (tm2 + 0.5) / 36525.0
t += (tm1 - 2451545.0) / 36525.0
rm_sun = Qreki._LONGITUDE_SUN(t)
rm_sun0 = 90 * int(rm_sun / 90.0)
delta_t1, delta_t2 = 1.0, 1.0
while abs(delta_t1 + delta_t2) > (1.0 / 86400.0):
t = (tm2 + 0.5) / 36525.0
t += (tm1 - 2451545.0) / 36525.0
rm_sun = Qreki._LONGITUDE_SUN(t)
delta_rm = rm_sun - rm_sun0
if delta_rm > 180.0:
delta_rm -= 360.0
elif delta_rm < -180.0:
delta_rm += 360.0
delta_t1 = int(delta_rm * 365.2 / 360.0)
delta_t2 = delta_rm * 365.2 / 360.0 - delta_t1
tm1 = tm1 - delta_t1
tm2 = tm2 - delta_t2
if tm2 < 0:
tm2 += 1.0
tm1 -= 1.0
nibun[0][0] = tm2 + 9.0 / 24.0 + tm1
nibun[0][1] = rm_sun0
@staticmethod
def _calc_saku(tm: float) -> float:
"""朔の時刻を計算する"""
tm1 = int(tm)
tm2 = tm - tm1
tm2 -= 9.0 / 24.0
delta_t1, delta_t2 = 1.0, 1.0
lc = 1
while abs(delta_t1 + delta_t2) > (1.0 / 86400.0):
t = (tm2 + 0.5) / 36525.0
t += (tm1 - 2451545.0) / 36525.0
rm_sun = Qreki._LONGITUDE_SUN(t)
rm_moon = Qreki._LONGITUDE_MOON(t)
delta_rm = rm_moon - rm_sun
if lc == 1 and delta_rm < 0.0:
delta_rm = Qreki._NORMALIZATION_ANGLE(delta_rm)
elif rm_sun >= 0 and rm_sun <= 20 and rm_moon >= 300:
delta_rm = Qreki._NORMALIZATION_ANGLE(delta_rm)
delta_rm = 360.0 - delta_rm
elif abs(delta_rm) > 40.0:
delta_rm = Qreki._NORMALIZATION_ANGLE(delta_rm)
delta_t1 = int(delta_rm * 29.530589 / 360.0)
delta_t2 = delta_rm * 29.530589 / 360.0 - delta_t1
tm1 = tm1 - delta_t1
tm2 = tm2 - delta_t2
if tm2 < 0.0:
tm2 += 1.0
tm1 -= 1.0
lc += 1
if lc == 15 and abs(delta_t1 + delta_t2) > (1.0 / 86400.0):
tm1 = int(tm - 26)
tm2 = 0
elif lc > 30 and abs(delta_t1 + delta_t2) > (1.0 / 86400.0):
return tm
return tm2 + tm1 + 9.0 / 24.0
@staticmethod
def _LONGITUDE_SUN(t: float) -> float:
"""太陽の黄経を計算する"""
th = 0.0
angles = [
(31557.0, 161.0), (29930.0, 48.0), (2281.0, 221.0), (155.0, 118.0),
(33718.0, 316.0), (9038.0, 64.0), (3035.0, 110.0), (65929.0, 45.0),
(22519.0, 352.0), (45038.0, 254.0), (445267.0, 208.0), (19.0, 159.0),
(32964.0, 158.0), (71998.1, 265.1), (35999.05, 267.52)
]
for coef, const in angles:
ang = Qreki._NORMALIZATION_ANGLE(coef * t + const)
th += 0.0004 * math.cos(Qreki.k * ang)
th += 0.0200 * math.cos(Qreki.k * Qreki._NORMALIZATION_ANGLE(71998.1 * t + 265.1))
th -= 0.0048 * t * math.cos(Qreki.k * Qreki._NORMALIZATION_ANGLE(35999.05 * t + 267.52))
th += 1.9147 * math.cos(Qreki.k * Qreki._NORMALIZATION_ANGLE(35999.05 * t + 267.52))
ang = Qreki._NORMALIZATION_ANGLE(36000.7695 * t + 280.4659)
th = Qreki._NORMALIZATION_ANGLE(th + ang)
return th
@staticmethod
def _LONGITUDE_MOON(t: float) -> float:
"""月の黄経を計算する"""
th = 0.0
angles = [
(2322131.0, 191.0), (4067.0, 70.0), (549197.0, 220.0), (1808933.0, 58.0),
(349472.0, 337.0), (381404.0, 354.0), (958465.0, 340.0), (12006.0, 187.0),
(39871.0, 223.0), (509131.0, 242.0), (1745069.0, 24.0), (1908795.0, 90.0),
(2258267.0, 156.0), (111869.0, 38.0), (27864.0, 127.0), (485333.0, 186.0),
(405201.0, 50.0), (790672.0, 114.0), (1403732.0, 98.0), (858602.0, 129.0),
(1920802.0, 186.0), (1267871.0, 249.0), (1856938.0, 152.0), (401329.0, 274.0),
(341337.0, 16.0), (71998.0, 85.0), (990397.0, 357.0), (818536.0, 151.0),
(922466.0, 163.0), (99863.0, 122.0), (1379739.0, 17.0), (918399.0, 182.0),
(1934.0, 145.0), (541062.0, 259.0), (1781068.0, 21.0), (133.0, 29.0),
(1844932.0, 56.0), (1331734.0, 283.0), (481266.0, 205.0), (31932.0, 107.0),
(926533.0, 323.0), (449334.0, 188.0), (826671.0, 111.0), (1431597.0, 315.0),
(1303870.0, 246.0), (489205.0, 142.0), (1443603.0, 52.0), (75870.0, 41.0),
(513197.9, 222.5), (445267.1, 27.9), (441199.8, 47.4), (854535.2, 148.2),
(1367733.1, 280.7), (377336.3, 13.2), (63863.5, 124.2), (966404.0, 276.5),
(35999.05, 87.53), (954397.74, 179.93), (890534.22, 145.7), (413335.35, 10.74),
(477198.868, 44.963)
]
for coef, const in angles:
ang = Qreki._NORMALIZATION_ANGLE(coef * t + const)
th += 0.0006 * math.cos(Qreki.k * ang)
th += 6.2888 * math.cos(Qreki.k * Qreki._NORMALIZATION_ANGLE(477198.868 * t + 44.963))
ang = Qreki._NORMALIZATION_ANGLE(481267.8809 * t + 218.3162)
th = Qreki._NORMALIZATION_ANGLE(th + ang)
return th
@staticmethod
def _YMDT2JD(year: int, month: int, day: int, hour: int, min: int, sec: int) -> float:
"""年月日時分秒からユリウス日を計算する"""
if month < 3:
year -= 1
month += 12
jd = int(365.25 * year)
jd += int(year / 400.0)
jd -= int(year / 100.0)
jd += int(30.59 * (month - 2.0))
jd += 1721088
jd += day
t = sec / 3600.0 + min / 60.0 + hour
t = t / 24.0
jd += t
return jd
@staticmethod
def _JD2YMDT(JD: float, TIME: List[int]) -> None:
"""ユリウス日から年月日時分秒を計算する"""
x0 = int(JD + 68570.0)
x1 = int(x0 / 36524.25)
x2 = x0 - int(36524.25 * x1 + 0.75)
x3 = int((x2 + 1) / 365.2425)
x4 = x2 - int(365.25 * x3) + 31.0
x5 = int(int(x4) / 30.59)
x6 = int(int(x5) / 11.0)
TIME[2] = x4 - int(30.59 * x5)
TIME[1] = x5 - 12 * x6 + 2
TIME[0] = 100 * (x1 - 49) + x3 + x6
if TIME[1] == 2 and TIME[2] > 28:
if TIME[0] % 100 == 0 and TIME[0] % 400 == 0:
TIME[2] = 29
elif TIME[0] % 4 == 0:
TIME[2] = 29
else:
TIME[2] = 28
tm = 86400.0 * (JD - int(JD))
TIME[3] = int(tm / 3600.0)
TIME[4] = int((tm - 3600.0 * TIME[3]) / 60.0)
TIME[5] = int(tm - 3600.0 * TIME[3] - 60 * TIME[4])
@staticmethod
def _rokuyou(month: int, day: int) -> str:
"""六曜を計算する"""
rokuyou_tbl = ["先勝", "友引", "先負", "仏滅", "大安", "赤口"]
x = (month + day - 2) % 6
return rokuyou_tbl[x]
@staticmethod
def _NORMALIZATION_ANGLE(angle: float) -> float:
"""角度の正規化を行う"""
if angle < 0.0:
angle1 = -angle
angle2 = int(angle1 / 360.0)
angle1 -= 360.0 * angle2
angle1 = 360.0 - angle1
else:
angle1 = int(angle / 360.0)
angle1 = angle - 360.0 * angle1
return angle1
@staticmethod
def from_gregorian(date: datetime.date) -> List[int]:
"""
グレゴリオ暦の日付から旧暦の年月日を計算する
:param date: datetime.date オブジェクト
:return: [年, 月, 日] の整数リスト
"""
tm = Qreki._YMDT2JD(date.year, date.month, date.day, 0, 0, 0)
q_reki = Qreki.calc_kyureki(tm)
return [q_reki[0], q_reki[2], q_reki[3]]
@staticmethod
def get_rokuyou(date: datetime.date) -> str:
"""
グレゴリオ暦の日付から六曜を計算する
:param date: datetime.date オブジェクト
:return: 六曜の文字列
"""
lunar_date = Qreki.from_gregorian(date)
return Qreki._rokuyou(lunar_date[1], lunar_date[2])
@staticmethod
def get_leap_month(date: datetime.date) -> bool:
"""
指定された日付が旧暦の閏月かどうかを判定する
:param date: datetime.date オブジェクト
:return: 閏月の場合True、そうでない場合False
"""
tm = Qreki._YMDT2JD(date.year, date.month, date.day, 0, 0, 0)
q_reki = Qreki.calc_kyureki(tm)
return q_reki[1] == 1
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment