实历冬至合朔表（BC722-BC104）

利用python的pyephem模块（需先安装）可以很容易获得合朔及冬至时间。

关于朔旦冬至的判断：冬至月判断可能存在的几个误差。

1. 现代测算误差：现代的天文计算，是近代天文学发展起来以后利用精密仪器测量天体的实际运行状况，由此构建的理论模型。但天体的运动状况非常复杂，当今的测量结果能难适用于久远前或久远后的宇宙，距离时间越远，误差越大。对于先秦时期的历表，应注意定朔或冬至时间在0-1时或23-24时段，可能几分钟的误差在进位时表现为1天。
2. 古代测算误差：战国前后是以冬至所在的月份作为正月，但月朔是无法观测的，推朔会有误差，而冬至的测量也存在误差。以《左传》两次“日南至”的记载为例，皆与今测算差2日。

历表的排布：每个冬至合朔点有3项数据（干支、日期、时间），其中时间部分主要用于参考可能的误差问题，分三行排列。每年从冬至月起排，如果到下一个冬至月出现了13次合朔，该年则排出13个月份，但并不表示该年存在闰月，月名仅供参考。古代的置闰会因为不同的历元和置闰方法有变化。如果连续两个月朔日干支的天干部分相同，说明前一个月为大月，否则为小月。如果连续三个月朔日干支的天干部分相同，说明前两个月为连大月。古历的朔策较今值为大，大小月的情况也与实历不符。

import ephem

nlrq = ["初一", "初二", "初三", "初四", "初五", "初六", "初七", "初八", "初九", "初十", "十一", "十二", "十三", "十四", "十五", "十六", "十七", "十八",
        "十九", "二十", "二一", "二二", "二三", "二四", "二五", "二六", "二七", "二八", "二九", "三十"]


def jd2ce(JD):
    JD = JD + 0.5  # 以-4812年1月1日0时为历元
    Z = int(JD)
    F = JD - Z
    if Z < 2299161:  # 儒略历
        A = Z
    else:  # 格里历
        a = int((Z - 1867216.25) / 36524.25)
        A = Z + 1 + a - int(a / 4)
    B = A + 1524
    C = int((B - 122.1) / 365.25)
    D = int(365.25 * C)
    E = int((B - D) / 30.6001)
    day = B - D - int(30.6001 * E) + F
    if E < 14:
        month = E - 1
    elif E < 16:
        month = E - 13
    if month > 2:
        year = C - 4716
    elif month in [1, 2]:
        year = C - 4715
    d = day - int(day)  # 取出一日的小数部分，转换为时分秒
    h = int(d * 24)
    m = int((d * 24 - h) * 60)
    s = d * 86400 - h * 3600 - m * 60
    if year <= 0:  # 将天文年表示为公历年
        year -= 1
    return year, month, int(day), h, m, round(s)


def ut8(t):  # 由UT0h转为UT8h
    jd0 = ephem.julian_date(t)
    jd = jd0 + 8 / 24
    ce = jd2ce(jd)
    return ce


def jd2gz(t):
    JD = ephem.julian_date(t)
    tiangan = ["甲", "乙", "丙", "丁", "戊", "己", "庚", "辛", "壬", "癸"]
    dizhi = ["子", "丑", "寅", "卯", "辰", "巳", "午", "未", "申", "酉", "戌", "亥"]
    gz = [''] * 60  # 六十甲子表
    for i in range(60):        gz[i] = tiangan[i % 10] + dizhi[i % 12]
    return gz[int(JD + 0.5 + 8 / 24 + 49) % 60]


def time_format(t, type):
    if type == 0:  # 日期
        return "BC{}-{:02d}-{:02d}".format(-ut8(t)[0], ut8(t)[1], ut8(t)[2])
    elif type == 1:  # 时间
        return "{:02d}:{:02d}:{:02d}".format(ut8(t)[3], ut8(t)[4], ut8(t)[5])


header = ["参考年", "冬至", '冬至日期', "冬至月朔", "二月（丑）", "三月（寅）", "四月", "五月", "六月", "七月", "八月", "九月", "十月", "十一月", "十二月",
          "十三月"]
dzhsb = open("冬至合朔表.xls", 'w')
for i in range(len(header)):
    dzhsb.write(header[i] + '\t')

for i in range(620):  # BC722-BC104
    year = -723 + i
    time = str(year) + '/12'
    dz = ephem.next_solstice(time)  # 冬至时间
    dz_jd = ephem.julian_date(dz)  # 冬至儒略日
    dzs0 = jd2ce(dz_jd - 29)  # 冬至前一月
    dzs = ephem.next_new_moon(dzs0)  # 预判的冬至朔
    dzs_jd = ephem.julian_date(dzs)  # 在UT8h下比较冬至朔是否为实际的冬至所在月朔日
    if int(dz_jd + 0.5 + 8 / 24) < int(dzs_jd + 0.5 + 8 / 24):
        dzs = ephem.next_new_moon(jd2ce(dz_jd - 30))
        dzs_jd = ephem.julian_date(dzs)
    # 写入上一年历表
    if i != 0:
        if str(heshuo) == str(dzs):
            hsb.pop()  # 上年最后一月朔为是年冬至朔，则删除
        dzhsb.write('\n')
        for m in range(3):
            dzhsb.write(hsb[0][0] + '\t')  # 参考年
            for n in range(len(hsb)):  # 月表
                if n == 0: continue
                if n == 2:  # 在冬至后插入冬至的日期
                    if m == 1: dzhsb.write(dzrq)
                    dzhsb.write('\t')
                dzhsb.write(hsb[n][m + 0] + '\t')
            dzhsb.write('\n')
        hsb.clear()
    # 生成历表
    hsb = [[str(year + 1)]]
    hsb.extend([[jd2gz(dz), time_format(dz, 0), time_format(dz, 1)]])  # 干支，年月日，时分秒
    hsb.extend([[jd2gz(dzs), time_format(dzs, 0), time_format(dzs, 1)]])
    dzrq = nlrq[int(dz_jd + 0.5 + 8 / 24) - int(dzs_jd + 0.5 + 8 / 24)]
    for j in range(12):  # 冬至以后的每月合朔表
        hs = jd2ce(dzs_jd + 29 * (j + 1))  # 估测的合朔日期
        heshuo = ephem.next_new_moon(hs)  # 实际合朔日期
        hsb.extend([[jd2gz(heshuo), time_format(heshuo, 0), time_format(heshuo, 1)]])