Julian Day Number conversion, Gregorian & Julian

julian_day.py

""" Julian day number to date conversion
    Proleptic Gregorian and Julian, with Astronomical years
    i.e., 1 AD = year 1, 1 BC = year 0, 2 BC = year -1, etc

    Derived from RG Tantzen (1963), ACM.
    Algorithm 199: conversions between calendar date and Julian day number.

    https://en.wikipedia.org/wiki/Julian_day
    Julian day number 0 assigned to the day starting at noon on
    January 1, 4713 BC, proleptic Julian calendar
    November 24, 4714 BC, proleptic Gregorian calendar

    Written by PM 2Ring 2022.01.28
"""

# Offset for Julian day numbers.
# With offset=0, day 0 is 0000-Feb-29
greg_offset = 1721119
julian_offset = 1721117

def greg_ymd_to_jdn(y, m, d, divmod=divmod):
    yy, m = divmod(m + 9, 12)
    y += yy - 1
    c, y = divmod(y, 100)
    d += greg_offset + c * 146097 // 4 + y * 1461 // 4
    return d + (m * 153 + 2) // 5

def jdn_to_greg_ymd(d, divmod=divmod):
    d = (d - greg_offset) * 4 - 1
    c, d = divmod(d, 146097)
    d = d // 4 * 4 + 3
    y, d = divmod(d, 1461)
    d = d // 4 * 5 + 2
    m, d = divmod(d, 153)
    d = d // 5 + 1
    yy, m = divmod(m + 2, 12)
    return 100 * c + y + yy, m + 1, d

def julian_ymd_to_jdn(y, m, d, divmod=divmod):
    yy, m = divmod(m + 9, 12)
    y += yy - 1
    d += julian_offset + y * 1461 // 4
    return d + (m * 153 + 2) // 5

def jdn_to_julian_ymd(d, divmod=divmod):
    d = (d - julian_offset) * 4 - 1
    d = d // 4 * 4 + 3
    y, d = divmod(d, 1461)
    d = d // 4 * 5 + 2
    m, d = divmod(d, 153)
    d = d // 5 + 1
    yy, m = divmod(m + 2, 12)
    return y + yy, m + 1, d

# Tests

if __name__ == '__main__':
    print("JDN 0")
    ymd = jdn_to_julian_ymd(0)
    print("Julian", ymd, ymd == (-4712, 1, 1))
    ymd = jdn_to_greg_ymd(0)
    print("Gregorian", ymd, ymd == (-4713, 11, 24))
    print()

    print("J2000.0 epoch")
    j2000 = (2000, 1, 1)
    jdn = greg_ymd_to_jdn(*j2000)
    ymd = jdn_to_greg_ymd(jdn)
    print("JDN", jdn, jdn == 2451545)
    print(ymd, ymd == j2000)
    print()

    def test_jdn(jdn_to_ymd, ymd_to_jdn, rng, verbose=False):
        print("Testing", jdn_to_ymd.__name__, ymd_to_jdn.__name__, rng, sep="\n")
        for jdn in rng:
            ymd = jdn_to_ymd(jdn)
            jd = ymd_to_jdn(*ymd)
            if verbose:
                print(jdn, ymd, jd)
            assert jdn == jd, (jdn, ymd, jd)
        print("ok\n")

    lo, hi, step = 0, 10_000_000, 101
    rng = range(lo, hi, step)
    verbose = False
    test_jdn(jdn_to_greg_ymd, greg_ymd_to_jdn, rng, verbose=verbose)
    test_jdn(jdn_to_julian_ymd, julian_ymd_to_jdn, rng, verbose=verbose)

    # Test against stdlib date
    # date.fromordinal(1) is 1 jan 1 AD = JDN 1721426
    # minimum lo = 1721426, maximum hi = 5373485
    def test_std(lo, hi, step=1):
        from datetime import date
        std_offset = 1721425
        rng = range(lo, hi, step)
        print("Testing against stdlib")
        print(rng)
        for jdn in rng:
            std = date.fromordinal(jdn - std_offset)
            std_ymd = std.year, std.month, std.day
            greg_jdn = greg_ymd_to_jdn(*std_ymd)
            assert jdn == greg_jdn, (jdn, std, greg_jdn)
        print("ok\n")

    # JDN of 1 jan 1 AD, Gregorian
    lo = 1721426
    # JDN of 1 jan 10,000 AD, Gregorian
    hi = 5373485
    test_std(lo, hi, 11)

Live version on SageMathCell

Live Converter interface, using Sage.
With weekday