ユリウス日とグレゴリウス暦日付との間の換算

こんにちは。
グレゴリウス暦日付 Gregorian calendar date $(year, month, day) $ と、修正ユリウス日（MJD）や Julian day number (JDN) との間の換算を計算してみました。下記のような変換計算の流れです。計算式については Wikipedia のユリウス通日#西暦からの換算 (Julian day#Calculation) そのままです^{1 2 3}。

(year, month, day) \leftrightarrow (y, m, d) \leftrightarrow n \leftrightarrow mjd  \leftrightarrow jdn

変換計算の中核部は、計算の都合上、導入した変数の間で行なっています。すなわち、$(y, m, d) \leftrightarrow n$ （関数 ymd2n、n2ymd）。

$ ./date.py --test
True: ymd2n(0,0,0) == 0
True: date2ymd(0,3,1) == (0,0,0)
True: date2mjd(1858, 11, 17) == 0

date.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division
from __future__ import print_function

# n: the continuous count of days since March first in 1 B.C. proleptic Gregorian
# date = (year, month, day) in the proleptic Gregorian calendar
# ymd  = (y, m, d); modified date in the proleptic Gregorian calendar (see date2ymd() below)

DELTA_N_MJD = {"Julian": 678883, "Gregorian": 678881}  # calendars

def mjd2jdn(mjd):
    return mjd + 2400001

def jdn2jd(jdn):
    return jdn - 0.5  # 00:00:00

def n2mjd(n):
    return n - DELTA_N_MJD["Gregorian"]

def mjd2n(mjd):
    return mjd + DELTA_N_MJD["Gregorian"]

def date2mjd(year, month, day):
    return n2mjd(date2n(year, month, day))

def mjd2date(mjd):
    return ymd2date(*n2ymd(mjd2n(mjd)))

def diffdate(date1, date2):
    return date2n(*date1) - date2n(*date2)

def adddate(date, d):
    n = date2n(*date) + d
    return ymd2date(n2ymd(n))

def date2n(year, month, day):
    return ymd2n(*date2ymd(year, month, day))

def date2ymd(year, month, day):
    y, m, d = year, month-3, day-1
    if m < 0:
        y, m = y-1, m%12
    return y, m, d

def ymd2date(y,m,d):
    year, month, day = y, m+3, d+1
    if month > 12:
        year, month = year+1, month%12
    return year, month, day

def ymd2n(y,m,d):
    n = d + (153*m+2)//5 + 365*y + y//4  # d + (306*m+4)//10 + (1461*y)//4
    n += -(y//100) + y//400  # (-3*(y//100))//4  [Gregorian]
    return n

def n2ymd(n):
    a = 4*n + 3
    a += 4*((3*((4*(n+1))//146097+1))//4)  # [Gregorian]
    y, b = a//1461, 5*((a%1461)//4) + 2
    m, d = b//153, (b%153)//5
    return y, m, d

def n2dayofweek(n):
   dow = ['sun', 'mon', 'tue', 'wed', 'thu', 'fri', 'sat']
   return dow[(n+3)%7]  # [Gregorian]

#===== test =====
def test_eval(str):
    print(eval(str), str)
    return 0
    
def test_epoch():
    test_eval('ymd2n(0,0,0) == 0')
    test_eval('date2ymd(0,3,1) == (0,0,0)')
    test_eval('date2mjd(1858, 11, 17) == 0')
    return 0

def test_ymd2n():
    for n in range(-146097,146097):
        y,m,d = n2ymd(n)
        if n != ymd2n(y,m,d):
            year, month, day = ymd2date(y,m,d)
            print(n,'', year, month, day, 'error')
            break
    return 0

def test_n2ymd():
    m, d = 0, 0
    for y in range(-401,400):
        n = ymd2n(y,m,d)
        if (y,m,d) != n2ymd(n):
            year, month, day = ymd2date(y,m,d)
            print(n,'', year, month, day, 'error')
            break
    return 0

def main():
    """
    {f}: Conversion between modified Julian day number and Gregorian calendar date.

    usage: {f} [-h] [--test]
    
    options:
        -h, --help    show this help message and exit
        --test        test
    """
    import docopt, textwrap
    args = docopt.docopt(textwrap.dedent(main.__doc__.format(f=__file__)))
    if args.get("--test", 0):
        test_epoch()
        test_ymd2n()
        test_n2ymd()
    return 0

if __name__ == '__main__':
    main()

1. またこの DELTA_N_MJD の Gregorian and Julian calendars の差 2 は、両暦の一致期間を西暦200年3月1日から100年間とすることから来ていて（西暦元年ではなく）、これは春分の日を3月21日に定めるとした西暦325年のニケーア宗教会議から来ています。 ↩

2. この値は1858年11月17日正子UT（jd=2400000.5）をMJD元期としていることからです（1858.680356 * 365.25 = 678883.000）。 ↩

3. なお、$(y, m, d) \leftrightarrow n$ の計算部分に表れていますが、一月平均日数が 30.6 (=153/5) 日、4年分の日数が 1461 日、グレゴリウス暦400年分の日数が 146097 日です。 ↩