Python で MJD
〜 修正ユリウス日を計算 〜
2023-07-07 作成 福島
TOP > asob > python-mjd
[ TIPS | TOYS | OTAKU | LINK | MOVIE | CGI | AvTitle | ConfuTerm | HIST | AnSt | Asob ]

本稿は他言語に対するプロトタイピングとして記述しています。

日付の加減算を装備していない他言語への実装や、学習のために本稿を参考にしてください。

本稿の内容は十進 BASIC 版と同じですが Python は付属モジュールが優れているので、紀元前を扱わない限り新規に MJD を作成する必要がありません。
(Modified Julian Day: 修正ユリウス日)

Python で日付の加減算をするには以下のように簡単に記述できます。
Python 3.6.8 (default, Sep 21 2021, 20:17:36)
[GCC 8.4.1 20200928 (Red Hat 8.4.1-1)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from datetime import datetime, timedelta          # 日付演算モジュールをインポートする。
>>> date = datetime(2023,7,7) - timedelta(days=10950) # 日数を減算する。
>>> date.year, date.month, date.day                   # 結果を表示する。
(1993, 7, 14)

扱う範囲は西暦 1 ～ 9999 年末。
Python 3.6.8 (default, Sep 21 2021, 20:17:36)
[GCC 8.4.1 20200928 (Red Hat 8.4.1-1)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from datetime import datetime                     # 日付演算モジュールをインポートする。
>>> datetime.min, datetime.max                        # 日付の最小値、最大値を表示する。
(datetime.datetime(1, 1, 1, 0, 0), datetime.datetime(9999, 12, 31, 23, 59, 59, 999999))

日付の加減算に関する問題

去年 (1 年前) の日付は「今年の西暦 - 1」を計算すれば求められます。
来年 (1 年後) の日付は「今年の西暦 + 1」ですね。
例: 2023-7-7 の 1 年前は 2022-7-7。

同様に、昨日の日付は「今日 - 1」で求められます。
例: 2023-7-7 の 1 日前は 2023-7-6。

では、10,950 日前の日付はどうでしょうか。(10,950 = 365 × 30)
本稿記述時の日付が 2023-7-7 なので、30 年前の日付は 1993-7-7。
···というわけには行きません。

1993 年 ～ 2023 年の間には「うるう年」が 7 回もあります。(1996, 2000, 2004, 2008, 2012, 2016, 2020)
このため、10,950 日前の日付は 1993-7-14 となります。

うるう年には、以下のような複雑なルールがあります。
  1. 4 年に一度はうるう年にする。
  2. 100 年に一度はうるう年にしない。
  3. 400 年に一度はうるう年にする。
西暦 2000 年は、上記 2 でうるう年にしないはずが、同 3 によってうるう年になります。
年月日の加減算に関する問題
時刻の計算をしてみます。

2 時間 35 分の映画が 18 時 30 分から上映されたら帰宅できるでしょうか。(例: シン・エヴァンゲリオン劇場版:||)
(途中で退席するとかはナシで)
答えは 21:05 なので、終電には余裕がありそうということが分かります。
未成年者の場合は、保護者の帯同が必要になるでしょう。

時計では 18 時 30 分 + 2 時間 35 分 = 21 時 5 分という計算をします。

同様に日付を計算して 2023-7-7 の 10 年 8 カ月 30 日後の日付はどうなるでしょう。
答えは 2033-3-4 ですが、上記うるう年のルールがあるため、即座に計算できる人はかなり少ないと思います。
(加算する 10 年 8 カ月 30 日間の中にもうるう年があることに注意)
それでも日付を計算したい
日付が暦(れき)のままだと計算が複雑になります。
一旦、基準日をかなり昔に設定して、そこから何日経過したか日数を揃えることにより、計算を単純にします。
その便利な基準日の定義が「修正ユリウス日」(MJD: Modified Julian Day) です。
「ユリウス通日」がオリジナルですが、そのままだと使いにくいので、これを使いやすくしたのが「修正ユリウス日」です。

すべての日付を「修正ユリウス日」という単純な整数に変換し、加減算が済んだのちに通常日付 (グレゴリオ暦) に戻す、という使い方をします。

Python で実装

アルゴリズムは Wikipedia からのパクリです。
Python の場合、床関数(ゆかかんすう)を演算子「//」(除算の商) で代用できますが、ここでは可読性を優先して math.floor() で記述しています。
int() を使うと負の値を除算したときの結果が floor() や // と異なるが、通常は負数を扱わない。

mjd.py
#!/usr/bin/python
# mjd.py version 1.0 written by fuku@rouge.gr.jp
#
# ユリウス日 / 修正ユリウス日 / グレゴリオ暦 の相互変換を行う。
#
# https://ja.wikipedia.org/wiki/ユリウス通日#修正ユリウス日（MJD）
# より。

""" Modified Julian Day (MJD)

Mutual conversion:
    Julian Date / Modified Julian Day / Gregorian Calendar
"""


def greg2mjd(y,m,d):
    """ Convert Gregorian-calendar to MJD.
    """
    from math import floor
    if m <= 2:
        m += 12
        y -= 1
    days = floor(365.25*y) + floor(y/400) - floor(y/100) \
         + floor(30.59*(m-2)) + d - 678912
    return days



def jd2mjd(y,m,d):
    """ Convert Julian-date to MJD.
    """
    from math import floor
    if m <= 2:
        m += 12
        y -= 1
    days = floor(365.25*y) + floor(30.59*(m-2)) + d - 678914
    return days



def mjd2jd(days):
    """ Convert MJD to Julian-date.
    """
    from math import floor
    n = days + 678883
    a = 4 * n + 3
    b = 5 * floor((a % 1461) / 4) + 2

    y = floor(a/1461)
    m = floor(b/153) + 3
    d = floor((b%153)/5) + 1

    if m > 12:
        m -= 12
        y += 1

    return y,m,d



def mjd2greg(days):
    """ Convert MJD to Gregorian-date.
    """
    from math import floor
    n = days + 678881
    a = 4*n + 3 + 4 * floor(3 / 4 * floor(4*(n+1)/146097+1))
    b = 5*floor((a % 1461)/4) + 2

    y = floor(a / 1461)
    m = floor(b / 153) + 3
    d = floor((b % 153)/5) + 1

    if m > 12:
        m -= 12
        y += 1

    return y,m,d



def mjd2week(days):
    """ Get day of week from MJD.
    """
    return days % 7     # 0:Wed, 1:Thu, .., 6:Tue



def mjd2zodiac(days):
    """ Get Oriental-Zodiac from MJD.
    """
    return days % 12    # 0:Tora, 1:U, 2:Tatsu, .., 11:Ushi



def mjd2stems(days):
    """ Get Heavenly-Stems from MJD.
    """
    return days % 10    # 0:Kinoe, 1:Kinoto, .., 9:Mizunoto



if __name__ == '__main__':

    y,m,d = 2012,1,1
    print('Gregorian y,m,d %d,%d,%d' % (y,m,d))

    mjd = greg2mjd(y,m,d)
    w   = mjd2week(mjd)
    print('mjd %d %s' % (mjd,['Wed','Thu','Fri','Sat','Sun','Mon','Tue'][w]))

    y,m,d = mjd2greg(mjd)
    print('Gregorian y,m,d %d,%d,%d' % (y,m,d))

    print('-' * 3)

    y,m,d = 1582,2,1
    print('Julian y,m,d %d,%d,%d' % (y,m,d))

    mjd = jd2mjd(y,m,d)
    print('mjd ' + str(mjd))

    y,m,d = mjd2jd(mjd)
    print('Julian y,m,d %d,%d,%d' % (y,m,d))