暦に興味を持ったので、旧暦計算プログラム QREKI.AWK を Python に移植してみた
前回、前々回と祝日について触れた際、暦に興味を持ったので Wikipedia の暦法関連を読み漁っていた。天文の話だけでなく、歴史や文化の話も混じっていておもしろい。
よくわかってなかった旧暦についても知ることができた。狭義の旧暦は今の日本の暦（グレゴリオ暦、ただしうるう年の計算基準が皇紀）の前の暦の意味で、明治5年12月2日まで使われていた天保暦。俗に言う旧暦は、天保暦の置閏法を借りてきて、今に無理やり適用したもの、と。
天保暦っぽいのを算出
そして、この旧暦を算出する QREKI.AWK 、およびその移植、派生品を発見。 .js, .pl, .php etc. 。でも qreki.py は検索しても見つからなかったので作ってみた。
>>> from qreki import Kyureki >>> k = Kyureki.from_ymd(2009, 11, 27) >>> print repr(k) Kyureki(2009, 10, 0, 11) >>> print k 2009年10月11日 >>> print k.rokuyou() 友引 >>> k.year, k.month, k.leap_month, k.day (2009, 10, 0, 11) >>> >>> from datetime import date >>> d = date(2009, 11, 27) >>> Kyureki.from_date(d) Kyureki(2009, 10, 0, 11)
できるのは旧暦を求めることと、その日の六曜 (大安, 赤口, 先勝, 友引, 先負, 仏滅)を求めること。
変換対象は datetime.date に収まる範囲。つまり proleptic Gregorian 暦*1の 1/1/1 から 9999/12/31 まで。この範囲内では QREKI.AWK と同じ結果が得られることを確認した。ただし、動作が同じということと、算出結果に誤りがないかどうかは別問題。詳しくはオリジナルの QREKI.DOC に。
9999年分の算出にかかった時間は、作者環境では、オリジナルの AWK 版で２日、ベタ移植した最初期の Python 版で 4 時間。
ベタ移植した時点では Python 的には無駄と思われる部分がいくつか見られたので少しずつ手直しして、動作が変わらないことを確認、の繰り返しで今に至る。 ver 0.4.2 は 9999 年算出に同環境で 1.5 時間。それなりに速くなったと思う。
コード
- qreki.py ver 0.4.6 ダウンロード (2009/12/1)
- qreki.py ver 0.4.5 ダウンロード (2009/11/30)
- qreki.py ver 0.4.4 ダウンロード (2009/11/28)
更新履歴
- ver 0.4.6 (2009/12/1)
- 内部関数 kyureki_from_jd, chuki_from_jd, before_nibun_from jd, saku_from_jd, jd2yearmonth も C 拡張化して高速化した。 C 拡張なしでも今までどおり動作する。 9999年算出速度は作者環境で 18 分。約 5 倍速。
- ver 0.4.5 (2009/11/30)
- コマンドラインからの使用もできるように。
C:\>python -m qreki -h Usage: qreki.py [year, [month, [day]]] Options: --version show program's version number and exit -h, --help show this help message and exit C:\>python -m qreki 2009年11月30日 2009年10月14日 C:\>python -m qreki 2009 12 1 2009年12月1日 2009年10月15日
- ver 0.4.4 (2009/11/30)
- docstring に少々追加。テストを行いやすくするため関数名の調整を行い pure Python 版と C 拡張のコードを共存させた (qreki._longitude_of_sun と qreki.longitude_of_sun の関係など ) 。
- ver 0.4.3 (2009/11/28)
- スピードアップ用の qreki.c および setup.py を追加。C 拡張なしでの動作もいままでどおり可能。 setup.py build して C 拡張を作っておくと動作が速くなる（作者環境で 4 倍速）。 C 拡張として作ったのは処理時間の 60 % 以上を占めていた太陽、月の黄経を計算する関数のみだが、十分な効果が得られたと思う。
- ver 0.4.2 (2009/11/28)
- 公開開始。
qreki.py (ver 0.4.6)
# coding: utf-8 from __future__ import division import math import datetime u"""新暦、旧暦変換スクリプト Arrenged by fgshun 2009 (http://d.hatena.ne.jp/fgshun/) 使用例 >>> from qreki import Kyureki >>> k = Kyureki.from_ymd(2009, 11, 27) >>> print repr(k) Kyureki(2009, 10, 0, 11) >>> print k 2009年10月11日 >>> print k.rokuyou() 友引 >>> k.year, k.month, k.leap_month, k.day (2009, 10, 0, 11) >>> >>> from datetime import date >>> d = date(2009, 11, 27) >>> Kyureki.from_date(d) Kyureki(2009, 10, 0, 11) このスクリプトで求められる『旧暦』とは、 天保暦法を参考にしつつも、中気や朔の刻の計算方法を 現代流の略算式に置き換えたりするなどの変更が加わっています。 詳しくはオリジナルに含まれる QREKI.DOC などを参考にしてください。 このスクリプトの内部で使われているローカル補正込みのユリウス通日とは 本来の Astronomical Julian Day (AJD) から、 - 0.5 + TZ(+9.0/24.0=0.375) ずらしたものであり、 JST における 0 時にて整数となるようになっています。 なお、 + 0.5 することで UTC における 0 時を整数とした Chronological Julian Day (CJD) とは補正の正負が異なります。 オリジナル 旧暦計算サンプルプログラム $Revision: 1.1 $ Coded by H.Takano 1993,1994 オリジナルは高野 英明氏の AWK 版です。下記より入手できます。 http://www.vector.co.jp/soft/dos/personal/se016093.html なお、これを移植したり、参考にして作成されたりしたものが いくつか存在するようです。 私が見つけたものをいくつか上げておきます （動作を確認していないものも含んでいます）。 JavaScript http://park1.wakwak.com/~y-nagano/Programs/koyomi/ Perl http://homepage2.nifty.com/sophia0/sub04.html Ruby http://www.310f.com/exocet/hiki/?jqreki.rb PHP http://www.2chan.net/script/ http://wppluginsj.sourceforge.jp/wp-koyomi/ JAVA http://homepage1.nifty.com/ave/kaihatsu/qreki.htm """ VERSION_INFO = (0, 4, 6) VERSION = u'.'.join(map(unicode, VERSION_INFO)) ORGINAL_VERSION_INFO = (1, 1) ORGINAL_VERSION = u'.'.join(map(unicode, ORGINAL_VERSION_INFO)) DEG_TO_RAD = math.pi / 180.0 # （角度の）度からラジアンに変換する係数 TZ = 0.375 # +9.0/24.0 (JST) __all__ = [ 'Kyureki', 'kyureki_from_ymd', 'kyureki_from_date', 'rokuyou_from_ymd', 'rokuyou_from_date'] class Kyureki(tuple): u"""旧暦を表すクラス""" __slots__ = () ROKUYOU = (u'大安', u'赤口', u'先勝', u'友引', u'先負', u'仏滅') def __new__(cls, year, month, leap_month, day): return super(Kyureki, cls).__new__( cls, (year, month, leap_month, day)) @classmethod def from_ymd(cls, year, month, day, tz=TZ): u"""新暦を表す 3 整数より旧暦を得る""" date = datetime.date(year, month, day) return cls.from_date(date, tz) @classmethod def from_date(cls, date, tz=TZ): u"""datetime.date より旧暦を得る""" jd = date2jd(date) kyureki = kyureki_from_jd(jd, tz, (date.year, date.month)) return cls(*kyureki) @property def year(self): u"""旧暦の年""" return self[0] @property def month(self): u"""旧暦の月""" return self[1] @property def leap_month(self): u"""閏月フラグ 0 ならば通常の月、 1 ならば閏月を表す。""" return self[2] @property def day(self): u"""旧暦の日""" return self[3] def rokuyou(self): u"""六曜を得る 戻り値は六曜 (大安, 赤口, 先勝, 友引, 先負, 仏滅) の文字列。""" return self.ROKUYOU[(self.month + self.day) % 6] def __repr__(self): return '%s(%r, %r, %r, %r)' % ( self.__class__.__name__, self.year, self.month, self.leap_month, self.day) def __unicode__(self): return u'%d年%s%d月%d日' % ( self.year, u'閏' if self.leap_month else u'', self.month, self.day) def __str__(self): return str(self.__unicode__()) kyureki_from_ymd = Kyureki.from_ymd kyureki_from_date = Kyureki.from_date # old name calc_kyureki = Kyureki.from_ymd def _kyureki_from_jd(tm, tz, shinreki_ym=None): u"""新暦に対応する、旧暦を求める 引数: tm0: ローカル補正込みのユリウス通日 shinreki_ym: 新暦年月（速度向上用オプション） tz: タイムゾーン 戻り値: 旧暦を表すタプル 0: 旧暦年 1: 旧暦月 2: 閏月フラグ (平月: 0, 閏月: 1） 3: 旧暦日""" tm0 = int(tm) # 計算対象の直前にあたる二分二至の時刻を求める # chu[0][0]:二分二至の時刻 chu[0][1]:その時の太陽黄経 chu = [] chu.append(before_nibun_from_jd(tm, tz)) # 中気の時刻を計算 for i in xrange(1, 4): chu.append(chuki_from_jd(chu[i - 1][0] + 32.0, tz)) # 計算対象の直前にあたる二分二至の直前の朔の時刻を求める saku = [] saku.append(saku_from_jd(chu[0][0], tz)) # 朔の時刻を求める for i in xrange(1, 5): saku.append(saku_from_jd(saku[i - 1] + 30.0, tz)) if abs(int(saku[i - 1]) - int(saku[i])) <= 26: # 前と同じ時刻を計算した場合（両者の差が26日以内）には、 # 初期値を +33日にして再実行させる saku[i] = saku_from_jd(saku[i - 1] + 35.0, tz) if int(saku[1]) <= int(chu[0][0]): # saku[1]が二分二至の時刻以前になってしまった場合には、 # 朔をさかのぼり過ぎたと考えて、朔の時刻を繰り下げて修正する。 # その際、計算もれ（saku[4]）になっている部分を補うため、 # 朔の時刻を計算する。 # （近日点通過の近辺で朔があると起こる事があるようだ...？） for i in xrange(0, 4): saku[i] = saku[i + 1] else: saku[4] = saku_from_jd(saku[3] + 35.0, tz) elif int(saku[0]) > int(chu[0][0]): # saku[0]が二分二至の時刻以後になってしまった場合には、 # 朔をさかのぼり足りないと見て、朔の時刻を繰り上げて修正する。 # その際、計算もれ（saku[0]）になっている部分を補うため、 # 朔の時刻を計算する。 # （春分点の近辺で朔があると起こる事があるようだ...？） for i in xrange(4, 0, -1): saku[i] = saku[i - 1] else: saku[0] = saku_from_jd(saku[0] - 27.0, tz) # 閏月検索Ｆｌａｇセット # （節月で４ヶ月の間に朔が５回あると、閏月がある可能性がある。） # leap=0:平月 leap=1:閏月 leap = 1 if int(saku[4]) <= int(chu[3][0]) else 0 # 朔日行列の作成 # m[i][0] ... 月名（1:正月 2:２月 3:３月 ....） # m[i][1] ... 閏フラグ（0:平月 1:閏月） # m[i][2] ... 朔日のjd m = [] for i in xrange(5): m.append([None, None, None]) m[0][0] = int(chu[0][1] / 30.0) + 2 m[0][1] = 0 m[0][2] = int(saku[0]) for i in xrange(1, 5): if leap == 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]) leap = 0 m[i][0] = m[i - 1][0] + 1 if m[i][0] > 12: m[i][0] -= 12 m[i][1] = 0 m[i][2] = int(saku[i]) # 朔日行列から旧暦を求める。 state = 0 for i in range(0, 5): if tm0 < m[i][2]: state = 1 break elif tm0 == m[i][2]: state = 2 break else: # break しなかったときに期待する i の値は 4 ではなく 5 i += 1 if state == 0 or state == 1: i -= 1 kyureki_month = m[i][0] kyureki_leap = m[i][1] kyureki_day = tm0 - m[i][2] + 1 # 旧暦年の計算 # （旧暦月が10以上でかつ新暦月より大きい場合には、 # まだ年を越していないはず...） if shinreki_ym: shinreki_year, shinreki_month = shinreki_ym else: # 引数に新暦の年月が与えられなかったので # ユリウス通日から逆算する shinreki_year, shinreki_month = jd2yearmonth(tm) kyureki_year = shinreki_year if kyureki_month > 9 and kyureki_month > shinreki_month: kyureki_year -= 1 return kyureki_year, kyureki_month, kyureki_leap, kyureki_day def _chuki_from_jd(tm, tz): u"""中気の時刻を求める 引数: tm: 計算対象となる時刻（ローカル補正込みのユリウス通日） tz: タイムゾーン 戻り値: 中気の時刻（ローカル補正込みのユリウス通日）と その時の黄経のタプル""" # 時刻引数を分解する tm2, tm1 = math.modf(tm) # JST ==> DT （補正時刻=0.0sec と仮定して計算） tm2 -= tz # 中気の黄経 λsun0 を求める t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0 rm_sun = longitude_of_sun(t) rm_sun0 = rm_sun - rm_sun % 30.0 # 繰り返し計算によって中気の時刻を計算する # （誤差が±1.0 sec以内になったら打ち切る。） delta_t1 = 0.0 delta_t2 = 1.0 while abs(delta_t1 + delta_t2) > 1.0 / 86400.0: # λsun を計算 t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0 rm_sun = longitude_of_sun(t) # 黄経差 Δλ＝λsun −λsun0 delta_rm = rm_sun - rm_sun0 # Δλの引き込み範囲（±180°）を逸脱した場合には、補正を行う if delta_rm > 180.0: delta_rm -= 360.0 elif delta_rm < -180.0: delta_rm += 360.0 # 時刻引数の補正値 Δt # delta_t = delta_rm * 365.2 / 360.0; delta_t2, delta_t1 = math.modf(delta_rm * 365.2 / 360.0) # 時刻引数の補正 # tm -= delta_t; tm1 -= delta_t1 tm2 -= delta_t2 if tm2 < 0.0: tm2 += 1.0 tm1 -= 1.0 return (tm1 + tm2 + tz, rm_sun0) def _before_nibun_from_jd(tm, tz): u"""直前の二分二至の時刻を求める 引数: tm: 計算対象となる時刻（ローカル補正込みのユリウス通日） tz: タイムゾーン 戻り値: 二分二至の時刻（ローカル補正込みのユリウス通日）と その時の黄経のタプル""" # 時刻引数を分解する tm2, tm1 = math.modf(tm) # JST ==> DT （補正時刻=0.0sec と仮定して計算） tm2 -= tz # 直前の二分二至の黄経 λsun0 を求める t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0 rm_sun = longitude_of_sun(t) rm_sun0 = rm_sun - rm_sun % 90.0 # 繰り返し計算によって直前の二分二至の時刻を計算する # （誤差が±1.0 sec以内になったら打ち切る。） delta_t1 = 0.0 delta_t2 = 1.0 while abs(delta_t1 + delta_t2) > (1.0 / 86400.0): # λsun を計算 t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0 rm_sun = longitude_of_sun(t) # 黄経差 Δλ＝λsun −λsun0 delta_rm = rm_sun - rm_sun0 # Δλの引き込み範囲（±180°）を逸脱した場合には、補正を行う if delta_rm > 180.0: delta_rm -= 360.0 elif delta_rm < -180.0: delta_rm += 360.0 # 時刻引数の補正値 Δt # delta_t = delta_rm * 365.2 / 360.0; delta_t2, delta_t1 = math.modf(delta_rm * 365.2 / 360.0) # 時刻引数の補正 # tm -= delta_t; tm1 -= delta_t1 tm2 -= delta_t2 if tm2 < 0.0: tm2 += 1.0 tm1 -= 1.0 # [0] 時刻引数を合成するのと、DT ==> JST 変換を行い、戻り値とする # （補正時刻=0.0sec と仮定して計算） # [1] 黄経 return (tm1 + tm2 + tz, rm_sun0) def _saku_from_jd(tm, tz): u"""与えられた時刻の直近の朔の時刻（JST）を求める 引数: tm: 計算対象となる時刻（ローカル補正込みのユリウス通日） 戻り値: 朔の時刻（ローカル補正込みのユリウス通日）""" # 時刻引数を分解する tm2, tm1 = math.modf(tm) # JST ==> DT （補正時刻=0.0sec と仮定して計算） tm2 -= tz # 繰り返し計算によって朔の時刻を計算する # （誤差が±1.0 sec以内になったら打ち切る。） delta_t1 = 0.0 delta_t2 = 1.0 for lc in xrange(1, 30): # 太陽の黄経λsun ,月の黄経λmoon を計算 # t = (tm - 2451548.0 + 0.5)/36525.0; t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0 rm_sun = longitude_of_sun(t) rm_moon = longitude_of_moon(t) # 月と太陽の黄経差Δλ # Δλ＝λmoon−λsun delta_rm = rm_moon - rm_sun if lc == 1 and delta_rm < 0.0: # ループの１回目（lc=1）で delta_rm < 0.0 の場合には # 引き込み範囲に入るように補正する delta_rm = delta_rm % 360.0 elif rm_sun >= 0.0 and rm_sun <= 20.0 and rm_moon >= 300.0: # 春分の近くで朔がある場合（0 ≦λsun≦ 20）で、 # 月の黄経λmoon≧300 の場合には、Δλ＝ 360.0 − Δλ と # 計算して補正する delta_rm = delta_rm % 360.0 delta_rm = 360.0 - delta_rm elif abs(delta_rm) > 40.0: # Δλの引き込み範囲（±40°）を逸脱した場合には、補正を行う delta_rm = delta_rm % 360.0 # 時刻引数の補正値 Δt # delta_t = delta_rm * 29.530589 / 360.0; delta_t2, delta_t1 = math.modf(delta_rm * 29.530589 / 360.0) # 時刻引数の補正 # tm -= delta_t; tm1 -= delta_t1 tm2 -= delta_t2 if tm2 < 0.0: tm2 += 1.0 tm1 -= 1.0 if abs(delta_t1 + delta_t2) > 1.0 / 86400.0: if lc == 15: # ループ回数が15回になったら、初期値 tm を tm-26 とする。 tm1 = tm - 26.0 tm2 = 0.0 else: # 朔の時刻の算出が終了した break else: # 初期値を補正したにも関わらず、 # 振動を続ける場合には異常終了させる。 raise ValueError(u'朔の計算が収束せず') # 時刻引数を合成するのと、DT ==> JST 変換を行い、戻り値とする # （補正時刻=0.0sec と仮定して計算） return tm1 + tm2 + tz def _longitude_of_sun(t): u"""太陽の黄経 λsun を計算する 引数: tm: 計算対象となる時刻 (AJD) / 36525.0 戻り値: 太陽の黄経 λsun""" k = DEG_TO_RAD cos = math.cos # 摂動項の計算 ang = (31557.0 * t + 161.0) % 360.0 th = .0004 * cos(k * ang) ang = (29930.0 * t + 48.0) % 360.0 th += .0004 * cos(k * ang) ang = (2281.0 * t + 221.0) % 360.0 th += .0005 * cos(k * ang) ang = (155.0 * t + 118.0) % 360.0 th += .0005 * cos(k * ang) ang = (33718.0 * t + 316.0) % 360.0 th += .0006 * cos(k * ang) ang = (9038.0 * t + 64.0) % 360.0 th += .0007 * cos(k * ang) ang = (3035.0 * t + 110.0) % 360.0 th += .0007 * cos(k * ang) ang = (65929.0 * t + 45.0) % 360.0 th += .0007 * cos(k * ang) ang = (22519.0 * t + 352.0) % 360.0 th += .0013 * cos(k * ang) ang = (45038.0 * t + 254.0) % 360.0 th += .0015 * cos(k * ang) ang = (445267.0 * t + 208.0) % 360.0 th += .0018 * cos(k * ang) ang = (19.0 * t + 159.0) % 360.0 th += .0018 * cos(k * ang) ang = (32964.0 * t + 158.0) % 360.0 th += .0020 * cos(k * ang) ang = (71998.1 * t + 265.1) % 360.0 th += .0200 * cos(k * ang) ang = (35999.05 * t + 267.52) % 360.0 th -= 0.0048 * t * cos(k * ang) th += 1.9147 * cos(k * ang) # 比例項の計算 ang = (36000.7695 * t) % 360.0 ang = (ang + 280.4659) % 360.0 th = (th + ang) % 360.0 return th def _longitude_of_moon(t): u"""月の黄経 λmoon を計算する 引数: tm: 計算対象となる時刻 (AJD) / 36525.0 戻り値: 月の黄経 λmoon""" k = DEG_TO_RAD cos = math.cos # 摂動項の計算 ang = (2322131.0 * t + 191.0) % 360.0 th = .0003 * cos(k * ang) ang = (4067.0 * t + 70.0) % 360.0 th += .0003 * cos(k * ang) ang = (549197.0 * t + 220.0) % 360.0 th += .0003 * cos(k * ang) ang = (1808933.0 * t + 58.0) % 360.0 th += .0003 * cos(k * ang) ang = (349472.0 * t + 337.0) % 360.0 th += .0003 * cos(k * ang) ang = (381404.0 * t + 354.0) % 360.0 th += .0003 * cos(k * ang) ang = (958465.0 * t + 340.0) % 360.0 th += .0003 * cos(k * ang) ang = (12006.0 * t + 187.0) % 360.0 th += .0004 * cos(k * ang) ang = (39871.0 * t + 223.0) % 360.0 th += .0004 * cos(k * ang) ang = (509131.0 * t + 242.0) % 360.0 th += .0005 * cos(k * ang) ang = (1745069.0 * t + 24.0) % 360.0 th += .0005 * cos(k * ang) ang = (1908795.0 * t + 90.0) % 360.0 th += .0005 * cos(k * ang) ang = (2258267.0 * t + 156.0) % 360.0 th += .0006 * cos(k * ang) ang = (111869.0 * t + 38.0) % 360.0 th += .0006 * cos(k * ang) ang = (27864.0 * t + 127.0) % 360.0 th += .0007 * cos(k * ang) ang = (485333.0 * t + 186.0) % 360.0 th += .0007 * cos(k * ang) ang = (405201.0 * t + 50.0) % 360.0 th += .0007 * cos(k * ang) ang = (790672.0 * t + 114.0) % 360.0 th += .0007 * cos(k * ang) ang = (1403732.0 * t + 98.0) % 360.0 th += .0008 * cos(k * ang) ang = (858602.0 * t + 129.0) % 360.0 th += .0009 * cos(k * ang) ang = (1920802.0 * t + 186.0) % 360.0 th += .0011 * cos(k * ang) ang = (1267871.0 * t + 249.0) % 360.0 th += .0012 * cos(k * ang) ang = (1856938.0 * t + 152.0) % 360.0 th += .0016 * cos(k * ang) ang = (401329.0 * t + 274.0) % 360.0 th += .0018 * cos(k * ang) ang = (341337.0 * t + 16.0) % 360.0 th += .0021 * cos(k * ang) ang = (71998.0 * t + 85.0) % 360.0 th += .0021 * cos(k * ang) ang = (990397.0 * t + 357.0) % 360.0 th += .0021 * cos(k * ang) ang = (818536.0 * t + 151.0) % 360.0 th += .0022 * cos(k * ang) ang = (922466.0 * t + 163.0) % 360.0 th += .0023 * cos(k * ang) ang = (99863.0 * t + 122.0) % 360.0 th += .0024 * cos(k * ang) ang = (1379739.0 * t + 17.0) % 360.0 th += .0026 * cos(k * ang) ang = (918399.0 * t + 182.0) % 360.0 th += .0027 * cos(k * ang) ang = (1934.0 * t + 145.0) % 360.0 th += .0028 * cos(k * ang) ang = (541062.0 * t + 259.0) % 360.0 th += .0037 * cos(k * ang) ang = (1781068.0 * t + 21.0) % 360.0 th += .0038 * cos(k * ang) ang = (133.0 * t + 29.0) % 360.0 th += .0040 * cos(k * ang) ang = (1844932.0 * t + 56.0) % 360.0 th += .0040 * cos(k * ang) ang = (1331734.0 * t + 283.0) % 360.0 th += .0040 * cos(k * ang) ang = (481266.0 * t + 205.0) % 360.0 th += .0050 * cos(k * ang) ang = (31932.0 * t + 107.0) % 360.0 th += .0052 * cos(k * ang) ang = (926533.0 * t + 323.0) % 360.0 th += .0068 * cos(k * ang) ang = (449334.0 * t + 188.0) % 360.0 th += .0079 * cos(k * ang) ang = (826671.0 * t + 111.0) % 360.0 th += .0085 * cos(k * ang) ang = (1431597.0 * t + 315.0) % 360.0 th += .0100 * cos(k * ang) ang = (1303870.0 * t + 246.0) % 360.0 th += .0107 * cos(k * ang) ang = (489205.0 * t + 142.0) % 360.0 th += .0110 * cos(k * ang) ang = (1443603.0 * t + 52.0) % 360.0 th += .0125 * cos(k * ang) ang = (75870.0 * t + 41.0) % 360.0 th += .0154 * cos(k * ang) ang = (513197.9 * t + 222.5) % 360.0 th += .0304 * cos(k * ang) ang = (445267.1 * t + 27.9) % 360.0 th += .0347 * cos(k * ang) ang = (441199.8 * t + 47.4) % 360.0 th += .0409 * cos(k * ang) ang = (854535.2 * t + 148.2) % 360.0 th += .0458 * cos(k * ang) ang = (1367733.1 * t + 280.7) % 360.0 th += .0533 * cos(k * ang) ang = (377336.3 * t + 13.2) % 360.0 th += .0571 * cos(k * ang) ang = (63863.5 * t + 124.2) % 360.0 th += .0588 * cos(k * ang) ang = (966404.0 * t + 276.5) % 360.0 th += .1144 * cos(k * ang) ang = (35999.05 * t + 87.53) % 360.0 th += .1851 * cos(k * ang) ang = (954397.74 * t + 179.93) % 360.0 th += .2136 * cos(k * ang) ang = (890534.22 * t + 145.7) % 360.0 th += .6583 * cos(k * ang) ang = (413335.35 * t + 10.74) % 360.0 th += 1.2740 * cos(k * ang) ang = (477198.868 * t + 44.963) % 360.0 th += 6.2888 * cos(k * ang) # 比例項の計算 ang = (481267.8809 * t) % 360.0 ang = (ang + 218.3162) % 360.0 th = (th + ang) % 360.0 return th def date2jd(date): u"""datetime.date からローカル補正込みのユリウス通日を得る""" return float(date.toordinal() + 1721424) def _jd2yearmonth(jd): u"""ローカル補正込みのユリウス通日から年月を得る""" f0 = math.floor(jd + 68570.0) f1 = math.floor(f0 / 36524.25) f2 = f0 - math.floor(36524.25 * f1 + 0.75) f3 = math.floor((f2 + 1.0) / 365.2425) f4 = f2 - math.floor(365.25 * f3) + 31.0 f5 = math.floor(f4 / 30.59) f6 = math.floor(f5 / 11.0) i1 = int(f1) i3 = int(f3) i5 = int(f5) i6 = int(f6) year = 100 * (i1 - 49) + i3 + i6 month = i5 - 12 * i6 + 2 return year, month # スピードアップ用の C 言語版が存在すればそちらを使う try: import _qreki kyureki_from_jd = _qreki._kyureki_from_jd chuki_from_jd = _qreki._chuki_from_jd before_nibun_from_jd = _qreki._before_nibun_from_jd saku_from_jd = _qreki._saku_from_jd longitude_of_sun = _qreki._longitude_of_sun longitude_of_moon = _qreki._longitude_of_moon jd2yearmonth = _qreki._jd2yearmonth except ImportError: kyureki_from_jd = _kyureki_from_jd chuki_from_jd = _chuki_from_jd before_nibun_from_jd = _before_nibun_from_jd saku_from_jd = _saku_from_jd longitude_of_sun = _longitude_of_sun longitude_of_moon = _longitude_of_moon jd2yearmonth = _jd2yearmonth def rokuyou_from_ymd(year, month, day): u"""六曜算出ショートカット 引数: 新暦年月日 戻り値: 六曜 (大安, 赤口, 先勝, 友引, 先負, 仏滅) の文字列 六曜を求めるにあたって、旧暦を算出している。 旧暦も必要とする場合、 Kyureki.from_ymd(year, month, day) で 旧暦オブジェクトをつくり、 これの rokuyou() メソッドを呼ぶほうが効率がよい。""" date = datetime.date(year, month, day) return rokuyou_from_date(date) # old name get_rokuyou = rokuyou_from_ymd def rokuyou_from_date(date): u"""六曜算出ショートカット 引数: datetime.date （新暦） 戻り値: 六曜 (大安, 赤口, 先勝, 友引, 先負, 仏滅) の文字列 六曜を求めるにあたって、旧暦を算出している。 旧暦も必要とする場合、 Kyureki.from_date(date) で 旧暦オブジェクトをつくり、 これの rokuyou() メソッドを呼ぶほうが効率がよい。""" kyureki = Kyureki.from_date(date) return kyureki.rokuyou() def main(): import optparse usage = u'%prog [year, [month, [day]]]' version = u'%%prog %s' % VERSION parser = optparse.OptionParser(usage=usage, version=version) option, args = parser.parse_args() try: def print_date(shinreki, kyureki): print u'%d年%d月%d日 %s' % ( shinreki.year, shinreki.month, shinreki.day, kyureki) len_args = len(args) if len_args == 0: d = datetime.date.today() k = Kyureki.from_date(d) print_date(d, k) elif len_args == 1: year = int(args[0]) d = datetime.date(year, 1, 1) d1 = datetime.timedelta(1) while d.year == year: k = Kyureki.from_date(d) print_date(d, k) d += d1 elif len_args == 2: year, month = map(int, args) d = datetime.date(year, month, 1) d1 = datetime.timedelta(1) while d.month == month: k = Kyureki.from_date(d) print_date(d, k) d += d1 elif len_args == 3: d = datetime.date(*map(int, args)) k = Kyureki.from_date(d) print_date(d, k) else: parser.error(u'引数が多すぎます') except StandardError, e: parser.error(unicode(e)) if __name__ == '__main__': main()
qreki.c
#include "Python.h" #include <stdlib.h> #include <math.h> static double normalize_angle(double angle); static PyObject * kyureki_from_jd(PyObject *self, PyObject *args, PyObject *kwargs); static PyObject * chuki_from_jd(PyObject *self, PyObject *args, PyObject *kwargs); static PyObject * before_nibun_from_jd(PyObject *self, PyObject *args, PyObject *kwargs); static PyObject * saku_from_jd(PyObject *self, PyObject *args, PyObject *kwargs); static PyObject * longitude_of_sun(PyObject *self, PyObject *args, PyObject *kwargs); static PyObject * longitude_of_moon(PyObject *self, PyObject *args, PyObject *kwargs); static PyObject * jd2yearmonth(PyObject *self, PyObject *args, PyObject *kwargs); static void chuki_from_jd2( double tm, double tz, double *chuki, double *longitude); static void before_nibun_from_jd2( double tm, double tz, double *nibun, double *longitude); static int saku_from_jd2(double tm, double tz, double *saku); static double longitude_of_sun2(double t); static double longitude_of_moon2(double t); static void jd2yearmonth2(double jd, int *year, int *month); static const double degToRad = Py_MATH_PI / 180.0; static double normalize_angle(double angle) { angle = fmod(angle, 360.0); if (angle < 0.0) { angle += 360.0; } return angle; } static PyObject * kyureki_from_jd(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"tm", "shinreki_ym", "tz", NULL}; double tm, tz; int shinreki_year=0, shinreki_month=0; int kyureki_year, kyureki_month, kyureki_leap, kyureki_day; int tm0; double chu[4][2]; double saku[5]; int m[4][3]; int leap; int i, state; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "dd|(ii)", kwlist, &tm, &tz, &shinreki_year, &shinreki_month)) return NULL; tm0 = (int)tm; before_nibun_from_jd2(tm, tz, &chu[0][0], &chu[0][1]); for (i=1; i < 4; i++) { chuki_from_jd2(chu[i-1][0] + 32.0, tz, &chu[i][0], &chu[i][1]); } if (saku_from_jd2(chu[0][0], tz, &saku[0]) == -1) return NULL; for (i=1; i < 5; i++) { if (saku_from_jd2(saku[i-1] + 30.0, tz, &saku[i]) == -1) return NULL; if (abs((int)saku[i - 1] - (int)saku[i]) <= 26) { if (saku_from_jd2(saku[i-1] + 35.0, tz, &saku[i]) == -1) return NULL; } } if ((int)(saku[1]) <= (int)(chu[0][0])) { for (i=0; i < 4; i++) saku[i] = saku[i+1]; if (saku_from_jd2(saku[3] + 35.0, tz, &saku[i]) == -1) return NULL; } else if((int)saku[0] > (int)chu[0][0]) { for (i=4; i > 0; i--) saku[i] = saku[i-1]; if (saku_from_jd2(saku[0] - 27.0, tz, &saku[i]) == -1) return NULL; } leap = ((int)saku[4] <= (int)chu[3][0]) ? 1 : 0; m[0][0] = (int)(chu[0][1] / 30.0) + 2; m[0][1] = 0; m[0][2] = (int)saku[0]; for (i=1; i < 5; i++) { if (leap == 1 && i != 1) { if ((int)chu[i-1][0] <= (int)saku[i-1] || (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]; leap = 0; } } m[i][0] = m[i-1][0] + 1; if (m[i][0] > 12) m[i][0] -= 12; m[i][1] = 0; m[i][2] = (int)saku[i]; } state = 0; for (i=0; i < 5; i++) { if (tm0 < m[i][2]) { state = 1; break; } else if (tm0 == m[i][2]) { state = 2; break; } } if (state == 0 || state == 1) i -= 1; kyureki_month = m[i][0]; kyureki_leap = m[i][1]; kyureki_day = tm0 - m[i][2] + 1; if (shinreki_year == 0) jd2yearmonth2(tm, &shinreki_year, &shinreki_month); kyureki_year = shinreki_year; if (kyureki_month > 9 && kyureki_month > shinreki_month) kyureki_year -= 1; return Py_BuildValue( "iiii", kyureki_year, kyureki_month, kyureki_leap, kyureki_day); } static PyObject * chuki_from_jd(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"tm", "tz", NULL}; double tm, tz; double chuki, longitude; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "dd", kwlist, &tm, &tz)) return NULL; chuki_from_jd2(tm, tz, &chuki, &longitude); return Py_BuildValue("dd", chuki, longitude); } static void chuki_from_jd2( double tm, double tz, double *chuki, double *longitude) { double tm1, tm2, t; double rm_sun, rm_sun0; double delta_rm, delta_t1, delta_t2; tm2 = modf(tm, &tm1); tm2 -= tz; t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0; rm_sun = longitude_of_sun2(t); rm_sun0 = rm_sun - fmod(rm_sun, 30.0); delta_t1 = 0.0; delta_t2 = 1.0; while (fabs(delta_t1 + delta_t2) > 1.0 / 86400.0) { t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0)/ 36525.0; rm_sun = longitude_of_sun2(t); delta_rm = rm_sun - rm_sun0; if (delta_rm > 180.0) { delta_rm -= 360.0; } else if (delta_rm < -180.0) { delta_rm += 360.0; } delta_t2 = modf(delta_rm * 365.2 / 360.0, &delta_t1); tm1 -= delta_t1; tm2 -= delta_t2; if (tm2 < 0.0) { tm2 += 1.0; tm1 -= 1.0; } } *chuki = tm1 + tm2 + tz; *longitude = rm_sun0; return; } static PyObject * before_nibun_from_jd(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"tm", "tz", NULL}; double tm, tz; double nibun, longitude; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "dd", kwlist, &tm, &tz)) return NULL; before_nibun_from_jd2(tm, tz, &nibun, &longitude); return Py_BuildValue("dd", nibun, longitude); } static void before_nibun_from_jd2( double tm, double tz, double *nibun, double *longitude) { double tm1, tm2, t; double rm_sun, rm_sun0; double delta_rm, delta_t1, delta_t2; tm2 = modf(tm, &tm1); tm2 -= tz; t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0; rm_sun = longitude_of_sun2(t); rm_sun0 = rm_sun - fmod(rm_sun, 90.0); delta_t1 = 0.0; delta_t2 = 1.0; while (fabs(delta_t1 + delta_t2) > 1.0 / 86400.0) { t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0)/ 36525.0; rm_sun = longitude_of_sun2(t); delta_rm = rm_sun - rm_sun0; if (delta_rm > 180.0) { delta_rm -= 360.0; } else if (delta_rm < -180.0) { delta_rm += 360.0; } delta_t2 = modf(delta_rm * 365.2 / 360.0, &delta_t1); tm1 -= delta_t1; tm2 -= delta_t2; if (tm2 < 0.0) { tm2 += 1.0; tm1 -= 1.0; } } *nibun = tm1 + tm2 + tz; *longitude = rm_sun0; return; } static PyObject * saku_from_jd(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"tm", "tz", NULL}; double tm, tz; double saku; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "dd", kwlist, &tm, &tz)) return NULL; if (saku_from_jd2(tm, tz, &saku) == -1) return NULL; return PyFloat_FromDouble(saku); } static int saku_from_jd2(double tm, double tz, double *saku) { double tm1, tm2, t; double rm_sun, rm_moon; double delta_rm, delta_t1, delta_t2; int lc; tm2 = modf(tm, &tm1); tm2 -= tz; delta_t1 = 0.0; delta_t2 = 1.0; for (lc = 1; lc < 30; lc++) { t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0)/ 36525.0; rm_sun = longitude_of_sun2(t); rm_moon = longitude_of_moon2(t); delta_rm = rm_moon - rm_sun; if (lc == 1 && delta_rm < 0.0) { delta_rm = normalize_angle(delta_rm); } else if (rm_sun >= 0.0 && rm_sun <= 20.0 && rm_moon >= 300.0) { delta_rm = normalize_angle(delta_rm); delta_rm = 360.0 - delta_rm; } else if (fabs(delta_rm) > 40.0) { delta_rm = normalize_angle(delta_rm); } delta_t2 = modf(delta_rm * 29.530589 / 360.0, &delta_t1); tm1 -= delta_t1; tm2 -= delta_t2; if (tm2 < 0.0) { tm2 += 1.0; tm1 -= 1.0; } if (fabs(delta_t1 + delta_t2) > 1.0 / 86400.0) { if (lc == 15) { tm1 = tm - 26.0; tm2 = 0.0; } } else { break; } } if (lc >= 30) { PyErr_SetString(PyExc_ValueError, ""); return -1; } *saku = tm1 + tm2 + tz; return 0; } static PyObject * longitude_of_sun(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"t", NULL}; double t; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "d", kwlist, &t)) return NULL; return PyFloat_FromDouble(longitude_of_sun2(t)); } static double longitude_of_sun2(double t) { double ang, th; ang = normalize_angle(31557.0 * t + 161.0); th = .0004 * cos(degToRad * ang); ang = normalize_angle(29930.0 * t + 48.0); th += .0004 * cos(degToRad * ang); ang = normalize_angle(2281.0 * t + 221.0); th += .0005 * cos(degToRad * ang); ang = normalize_angle(155.0 * t + 118.0); th += .0005 * cos(degToRad * ang); ang = normalize_angle(33718.0 * t + 316.0); th += .0006 * cos(degToRad * ang); ang = normalize_angle(9038.0 * t + 64.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(3035.0 * t + 110.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(65929.0 * t + 45.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(22519.0 * t + 352.0); th += .0013 * cos(degToRad * ang); ang = normalize_angle(45038.0 * t + 254.0); th += .0015 * cos(degToRad * ang); ang = normalize_angle(445267.0 * t + 208.0); th += .0018 * cos(degToRad * ang); ang = normalize_angle(19.0 * t + 159.0); th += .0018 * cos(degToRad * ang); ang = normalize_angle(32964.0 * t + 158.0); th += .0020 * cos(degToRad * ang); ang = normalize_angle(71998.1 * t + 265.1); th += .0200 * cos(degToRad * ang); ang = normalize_angle(35999.05 * t + 267.52); th -= 0.0048 * t * cos(degToRad * ang); th += 1.9147 * cos(degToRad * ang); ang = normalize_angle(36000.7695 * t); ang = normalize_angle(ang + 280.4659); th = normalize_angle(th + ang); return th; } static PyObject * longitude_of_moon(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"t", NULL}; double t; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "d", kwlist, &t)) return NULL; return PyFloat_FromDouble(longitude_of_moon2(t)); } static double longitude_of_moon2(double t) { double ang, th; ang = normalize_angle(2322131.0 * t + 191.0); th = .0003 * cos(degToRad * ang); ang = normalize_angle(4067.0 * t + 70.0); th += .0003 * cos(degToRad * ang); ang = normalize_angle(549197.0 * t + 220.0); th += .0003 * cos(degToRad * ang); ang = normalize_angle(1808933.0 * t + 58.0); th += .0003 * cos(degToRad * ang); ang = normalize_angle(349472.0 * t + 337.0); th += .0003 * cos(degToRad * ang); ang = normalize_angle(381404.0 * t + 354.0); th += .0003 * cos(degToRad * ang); ang = normalize_angle(958465.0 * t + 340.0); th += .0003 * cos(degToRad * ang); ang = normalize_angle(12006.0 * t + 187.0); th += .0004 * cos(degToRad * ang); ang = normalize_angle(39871.0 * t + 223.0); th += .0004 * cos(degToRad * ang); ang = normalize_angle(509131.0 * t + 242.0); th += .0005 * cos(degToRad * ang); ang = normalize_angle(1745069.0 * t + 24.0); th += .0005 * cos(degToRad * ang); ang = normalize_angle(1908795.0 * t + 90.0); th += .0005 * cos(degToRad * ang); ang = normalize_angle(2258267.0 * t + 156.0); th += .0006 * cos(degToRad * ang); ang = normalize_angle(111869.0 * t + 38.0); th += .0006 * cos(degToRad * ang); ang = normalize_angle(27864.0 * t + 127.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(485333.0 * t + 186.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(405201.0 * t + 50.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(790672.0 * t + 114.0); th += .0007 * cos(degToRad * ang); ang = normalize_angle(1403732.0 * t + 98.0); th += .0008 * cos(degToRad * ang); ang = normalize_angle(858602.0 * t + 129.0); th += .0009 * cos(degToRad * ang); ang = normalize_angle(1920802.0 * t + 186.0); th += .0011 * cos(degToRad * ang); ang = normalize_angle(1267871.0 * t + 249.0); th += .0012 * cos(degToRad * ang); ang = normalize_angle(1856938.0 * t + 152.0); th += .0016 * cos(degToRad * ang); ang = normalize_angle(401329.0 * t + 274.0); th += .0018 * cos(degToRad * ang); ang = normalize_angle(341337.0 * t + 16.0); th += .0021 * cos(degToRad * ang); ang = normalize_angle(71998.0 * t + 85.0); th += .0021 * cos(degToRad * ang); ang = normalize_angle(990397.0 * t + 357.0); th += .0021 * cos(degToRad * ang); ang = normalize_angle(818536.0 * t + 151.0); th += .0022 * cos(degToRad * ang); ang = normalize_angle(922466.0 * t + 163.0); th += .0023 * cos(degToRad * ang); ang = normalize_angle(99863.0 * t + 122.0); th += .0024 * cos(degToRad * ang); ang = normalize_angle(1379739.0 * t + 17.0); th += .0026 * cos(degToRad * ang); ang = normalize_angle(918399.0 * t + 182.0); th += .0027 * cos(degToRad * ang); ang = normalize_angle(1934.0 * t + 145.0); th += .0028 * cos(degToRad * ang); ang = normalize_angle(541062.0 * t + 259.0); th += .0037 * cos(degToRad * ang); ang = normalize_angle(1781068.0 * t + 21.0); th += .0038 * cos(degToRad * ang); ang = normalize_angle(133.0 * t + 29.0); th += .0040 * cos(degToRad * ang); ang = normalize_angle(1844932.0 * t + 56.0); th += .0040 * cos(degToRad * ang); ang = normalize_angle(1331734.0 * t + 283.0); th += .0040 * cos(degToRad * ang); ang = normalize_angle(481266.0 * t + 205.0); th += .0050 * cos(degToRad * ang); ang = normalize_angle(31932.0 * t + 107.0); th += .0052 * cos(degToRad * ang); ang = normalize_angle(926533.0 * t + 323.0); th += .0068 * cos(degToRad * ang); ang = normalize_angle(449334.0 * t + 188.0); th += .0079 * cos(degToRad * ang); ang = normalize_angle(826671.0 * t + 111.0); th += .0085 * cos(degToRad * ang); ang = normalize_angle(1431597.0 * t + 315.0); th += .0100 * cos(degToRad * ang); ang = normalize_angle(1303870.0 * t + 246.0); th += .0107 * cos(degToRad * ang); ang = normalize_angle(489205.0 * t + 142.0); th += .0110 * cos(degToRad * ang); ang = normalize_angle(1443603.0 * t + 52.0); th += .0125 * cos(degToRad * ang); ang = normalize_angle(75870.0 * t + 41.0); th += .0154 * cos(degToRad * ang); ang = normalize_angle(513197.9 * t + 222.5); th += .0304 * cos(degToRad * ang); ang = normalize_angle(445267.1 * t + 27.9); th += .0347 * cos(degToRad * ang); ang = normalize_angle(441199.8 * t + 47.4); th += .0409 * cos(degToRad * ang); ang = normalize_angle(854535.2 * t + 148.2); th += .0458 * cos(degToRad * ang); ang = normalize_angle(1367733.1 * t + 280.7); th += .0533 * cos(degToRad * ang); ang = normalize_angle(377336.3 * t + 13.2); th += .0571 * cos(degToRad * ang); ang = normalize_angle(63863.5 * t + 124.2); th += .0588 * cos(degToRad * ang); ang = normalize_angle(966404.0 * t + 276.5); th += .1144 * cos(degToRad * ang); ang = normalize_angle(35999.05 * t + 87.53); th += .1851 * cos(degToRad * ang); ang = normalize_angle(954397.74 * t + 179.93); th += .2136 * cos(degToRad * ang); ang = normalize_angle(890534.22 * t + 145.7); th += .6583 * cos(degToRad * ang); ang = normalize_angle(413335.35 * t + 10.74); th += 1.2740 * cos(degToRad * ang); ang = normalize_angle(477198.868 * t + 44.963); th += 6.2888 * cos(degToRad * ang); ang = normalize_angle(481267.8809 * t); ang = normalize_angle(ang + 218.3162); th = normalize_angle(th + ang); return th; } static PyObject * jd2yearmonth(PyObject *self, PyObject *args, PyObject *kwargs) { static char *kwlist[] = {"jd", NULL}; double jd; int year, month; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "d", kwlist, &jd)) return NULL; jd2yearmonth2(jd, &year, &month); return Py_BuildValue("ii", year, month); } static void jd2yearmonth2(double jd, int *year, int *month) { double f0, f1, f2, f3, f4, f5, f6; int i1, i3, i5, i6; f0 = floor(jd + 68570.0); f1 = floor(f0 / 36524.25); f2 = f0 - floor(36524.25 * f1 + 0.75); f3 = floor((f2 + 1.0) / 365.2425); f4 = f2 - floor(365.25 * f3) + 31.0; f5 = floor(f4 / 30.59); f6 = floor(f5 / 11.0); i1 = (int)f1; i3 = (int)f3; i5 = (int)f5; i6 = (int)f6; *year = 100 * (i1 - 49) + i3 + i6; *month = i5 - 12 * i6 + 2; return; } static PyMethodDef module_methods[] = { {"_kyureki_from_jd", (PyCFunction)kyureki_from_jd, METH_VARARGS | METH_KEYWORDS, "calculate kyureki."}, {"_chuki_from_jd", (PyCFunction)chuki_from_jd, METH_VARARGS | METH_KEYWORDS, "calculate chuki."}, {"_before_nibun_from_jd", (PyCFunction)before_nibun_from_jd, METH_VARARGS | METH_KEYWORDS, "calculate before nibun."}, {"_saku_from_jd", (PyCFunction)saku_from_jd, METH_VARARGS | METH_KEYWORDS, "calculate saku."}, {"_longitude_of_sun", (PyCFunction)longitude_of_sun, METH_VARARGS | METH_KEYWORDS, "calculate longitude of sun."}, {"_longitude_of_moon", (PyCFunction)longitude_of_moon, METH_VARARGS | METH_KEYWORDS, "calculate longitude of moon."}, {"_jd2yearmonth", (PyCFunction)jd2yearmonth, METH_VARARGS | METH_KEYWORDS, "calculate shinreki year and shinreki month from jd."}, {NULL, NULL, 0, NULL} /* Sentinel */ }; PyMODINIT_FUNC init_qreki(void) { PyObject* m; m = Py_InitModule3("_qreki", module_methods, "calculate longitude of sun, moon."); if (m == NULL) return; }
setup.py
from distutils.core import setup, Extension py_modules = ['qreki'] ext_modules = [ Extension('_qreki', sources = ['qreki.c'])] setup( name='qreki', version='0.4.6', description='Qreki', py_modules=py_modules, ext_modules=ext_modules)
自分用メモ、ベタ移植 ver 0.0 からの修正の概要
処理時間の多くを占めていたのは太陽、月の黄経の算出の箇所。ここを中心に速度アップをめざした。
まずは、角度を 360 度内に正規化する
def NORMALIZATION_ANGLE(angle): u"""角度の正規化を行う。すなわち引数の範囲を ０≦θ＜３６０ にする。""" if angle < 0.0: angle1 = -angle angle2 = int(angle1 / 360.0) angle1 -= 360 * angle2 angle1 = 360.0 - angle1 else: angle1 = int(angle / 360.0) angle1 = angle - 360.0 * angle1 return angle1
は % 360.0 で事足りるので置き換えた。
そして、グローバル変数に何十回とアクセスする場合、ローカル変数に代入しておくようにした。
係数をタプル化してジェネレータと sum 関数を使うようにしようとしたが止めた。似たような式の集まりだがベタ書きしたほうが動作は速い。
その他の箇所では、整数部、小数部に分けるのに
tm1 = int(tm)
tm2 = tm - tm1
しているので math.modf に変更したり、小数を整数に戻したのにすぐに浮動小数点と掛けたりするようなとこを math.floor に置き換えたり、小数のままにしてみたり、などなど。
自分用メモ、より速くするために
目的を「ある日の旧暦を求める」ではなく「ある期間の旧暦の一覧を作る」に変更するならば、一度求めた中気の時刻、朔の時刻を保持しておくことで効率を上げることは可能。
太陽、月の黄経の算出部分を C 言語で書く。実もふたもない対策だが効果は高そう（追記、 C 言語で書き終えた。約 4 倍の速度が得られた）。
アンダーフロー対策の整数部、小数部の分解を止めれば、精度は犠牲になるが少し速くなるはず。
2009/11/29 追記。角度ラジアン間の変換。 Python では math.radians とかどうだろう？ やっていることは同じではあるが…（30日追記、 deg_to_rad = math.pi / 180.0 しておいて deg_to_rad * deg のほうが速い）。 C のほうはこのままでよいが、関数マクロ化すると見た目だけはよくなるかも。
