计算给定日期前后的节气具体时间

从小到大，每到过年前使用农历的频次会显著增多，比如腊八、小年、除夕、初一，还有腊月的年俗。注意我们使用的农历不是阴历也不是阳历，而是阴阳合历。农历是根据月相的变化周期，每一次月相朔望变化为一个月，参考太阳回归年为一年的长度，并加入二十四节气与设置闰月以使平均历年与回归年相适应。

采用阴阳合历主要是要让农业生产配合对的农时，保证最大化利用自然天气变化，保证一年的农业收成。那么，如何比较准确的计算农历呢？这里面不得不说二十四节气。一首耳熟能详的歌谣，“春雨惊春清谷天，夏满芒夏暑相连。秋处露秋寒霜降，冬雪雪冬小大寒。”二十四节气对应的都是重要的农时。之前也就是当作一首歌谣来记得，但是要想知道每个节气具体到来的时间，除了查黄历，要准确计算该怎么办呢？

廿四节气原依据北斗七星的斗转星移制定，即所谓“斗柄指东，天下皆春；斗柄指南，天下皆夏；斗柄指西，天下皆秋；斗柄指北，天下皆冬”的星象规律。

现行的“二十四节气”来自于三百多年前依据太阳黄经度数划分的方法，自1645年起沿用至今。“二十四节气”是根据太阳在回归黄道上的位置（太阳黄经度数）来定的，即在一个为360度圆周的“黄道”（一年当中太阳在天球上的视路径）上，划分为24等份，每15°一等份，以春分点为0度起点，按黄经度数编排。即，视太阳从春分点，也就是黄经零度出发，此刻太阳垂直照射赤道，每前进15度为一个节气；运行一周又回到春分点，为一回归年，24个节气正好360度，太阳在黄道上每运行15度为一个“节气”。

因此，除了古人使用的北斗七星斗柄指向、“土圭测影”等方法，以及民间一些特殊现象的经验性总结，使用太阳黄经划分是最精确的方法。

接下来，我们该如何获得太阳黄经呢？传统方法是根据天文记录和官方给出的星图，比如NASA，南京紫金山天文台等专业机构。

Ephem是一个用于计算天文学相关信息的Python库，它提供了丰富的功能来处理天文数据，包括恒星位置、日月位置、行星位置、星座边界、日出日落时间等。这个库最初由B.C. Mills开发，后来被Paul Schlyter扩展和更新。

以计算给定日期（公历本地时间）的前一个和后一个节气为例，我们的计算过程大概可以分为几步，首先将根据日期计算得到对应的太阳黄经度数，然后以15度为一个节气，计算出前一个和后一个节气的目标太阳黄经度数，接下来以当前时间为起点，搜索达到目标经度时的本地时间。

由于Ephem库使用的是UT时间，注意先将本地时间换算到UT时间，如北京时间-8小时。

由于一年中节对的是农历的12个月，因此我们以找节为例，找气类似。第一步和第二步比较简单，第三步稍微复杂一些。下面再深入讲讲思路，以找当前日期的后一个节时间为例。

确定目标经度后，从当前日期/经度开始，求目标经度的时间呢？已知是当前日期，当前经度，目标经度，求目标节时间。最简单的方法就是模拟天体运动，转到对应的角度。那就相当于从当前日期开始，以一个步长推进时间，然后计算模拟推进日期的经度，如果达到目标经度，那此时的时间就是节时间。

按照这个思路实现后，发现解的节时间精度不够，对比紫金山天文台给的数据最多的时候能差出来2天。那怎么办？查了查资料，提高精度的方法也很简单，在上述方法的基础上，在初步节时间基础上，再在一个小范围内搜索。第二阶段的搜索，由于搜索范围缩小了，为了提高效率，可以不再模拟小步长推进，而是使用二分查找。划一个搜索范围，看中点时间的经度，如果比目标经度小，那就再中点和上界中搜索，如果比目标经度大，那就再下界和中点中搜索。迭代多次。

经过这样一个折腾，得到的节气时间能够精确的分钟。这个精度对于自用，或者基于此做万年历给节气时间足够了，毕竟不会出现天为单位偏差。但是，根据Ephem库的说明，由于使用的历表数据限制，同时使用的是快速算法，其给出的黄经精度会随着时间偏移增大而增加。

最后，放上代码供参考。

import math
import ephem
from datetime import datetime, timedelta
from typing import Union

# 节气
JIEQI = [
    "立春", "雨水", "惊蛰", "春分", "清明", "谷雨", 
    "立夏", "小满", "芒种", "夏至", "小暑", "大暑", 
    "立秋", "处暑", "白露", "秋分", "寒露", "霜降", 
    "立冬", "小雪", "大雪", "冬至", "小寒", "大寒"]


# 中文和英文混杂的对齐有问题，需要特殊的填充
def get_chinese_char_number(char):
    count = 0
    for item in char:
       # chinese char and chinese punctuation mark
        if 0x4E00 <= ord(item) <= 0x9FA5 or 0xFF00 <= ord(item) <=0xFFEF or 0x3000 <= ord(item) <= 0x303F:
            count += 1
    return count


def format_with_padding(char, align, length):
    if align=='l' or align=='left' :
        anchor = '<'
    elif align=='c' or align=='center':
        anchor = '^'
    elif align=='r' or align=='right':
        anchor = '>'
    else:
        ut.print_error('Not support align type. Current support is l(left), c(center), r(right)')
        return -1
    add_len = get_chinese_char_number(char)
    if add_len >= length:
        p_len = 0
    else:
        p_len = length - add_len
    return f'{char:{anchor}{p_len}}'


# 使用Ephem计算节气, 准确. 偏差在几个小时
# 计算太阳黄经hlon
def get_sun_longitude_hlon(date, observer):
    observer.date = date
    sun = ephem.Sun(observer)
    sun_lon = sun.hlon * 180 / ephem.pi  # 转换为度数
    return sun_lon % 360  # 确保在 0-360范围内


# 计算上一个节气
def find_previous_jieqi_hlon(date, observer):
    current_date = date
    cur_lon = get_sun_longitude_hlon(current_date, observer)
    # print(current_date, cur_lon)

    jieqi_index = int(cur_lon) // 15 - 1 if int(cur_lon) % 15 >= 0 else 0
    jieqi_index = jieqi_index % 24
    if jieqi_index%2 == 0: # 修正到节, 立秋315度为节, 对应21, 为奇数. 因此, 遇到偶数是气, 需要找到下一个节.
        jieqi_index -= 1

    # print(jieqi_index)
    target_lon = (jieqi_index) * 15 # 目标黄经(0起每15一个节气)
    # print(target_lon)

    # 第一阶段：找到黄经首次超过目标值的日期
    date0 = ephem.Date(f"{date.year}/{date.month}/{date.day}")
    # print(ephem_to_datetime(date0))
    observer0 = ephem.Observer()
    observer0.lat = observer.lat
    observer0.lon = observer.lon

    while True:
        current_lon = get_sun_longitude_hlon(date0, observer0) # 计算视黄经
        # print(ephem_to_datetime(date0), current_lon)

        if abs(current_lon - target_lon) < 0.1:
            break
        date0 -= 1./24  # 逐小时逼近

    # print('First phase', ephem_to_datetime(date0))

    # 第二阶段：二分法精确到秒
    lower = date0 - 2  # 前一天
    upper = date0 + 2  # 后两天
    # print(ephem_to_datetime(lower), ephem_to_datetime(upper))
    for _ in range(40):  # 20次二分可达约0.1秒精度
        mid = (lower + upper) * 0.5
        current_lon = get_sun_longitude_hlon(mid, observer0) # 计算视黄经
        # print(ephem_to_datetime(mid), current_lon)

        if current_lon < target_lon:
            lower = mid
        else:
            upper = mid

    # print('Second phase', ephem_to_datetime(lower), ephem_to_datetime(mid), ephem_to_datetime(upper))
    jieqi_index = (jieqi_index-9)%24
    jieqi = JIEQI[jieqi_index] # offset
    jieqi_time = ephem_to_datetime(upper)
    return jieqi, jieqi_time, jieqi_index


# 计算下一个节气
def find_next_jieqi_hlon(date, observer):
    current_date = date
    cur_lon = get_sun_longitude_hlon(current_date, observer)
    # print(current_date, cur_lon)

    jieqi_index = int(cur_lon) // 15 + 1 if int(cur_lon) % 15 >= 0 else 0
    jieqi_index = jieqi_index % 24
    if jieqi_index%2 == 0: # 修正到节, 立秋315度为节, 对应21, 为奇数. 因此, 遇到偶数是气, 需要找到下一个节.
        jieqi_index += 1

    # print(jieqi_index)
    target_lon = (jieqi_index) * 15 # 目标黄经(0起每15一个节气)
    # print(target_lon)

    # 第一阶段：找到黄经首次超过目标值的日期
    date0 = ephem.Date(f"{date.year}/{date.month}/{date.day}")
    # print(ephem_to_datetime(date0))
    observer0 = ephem.Observer()
    observer0.lat = observer.lat
    observer0.lon = observer.lon

    while True:
        current_lon = get_sun_longitude_hlon(date0, observer0) # 计算视黄经
        # print(ephem_to_datetime(date0), current_lon)

        if abs(current_lon - target_lon) < 0.1:
            break
        date0 += 1./24  # 逐小时逼近

    # print('First phase', ephem_to_datetime(date0))

    # 第二阶段：二分法精确到秒
    lower = date0 - 2  # 前一天
    upper = date0 + 2  # 后两天
    # print(ephem_to_datetime(lower), ephem_to_datetime(upper))
    for _ in range(40):  # 20次二分可达约0.1秒精度
        mid = (lower + upper) * 0.5
        current_lon = get_sun_longitude_hlon(mid, observer0) # 计算视黄经
        # print(ephem_to_datetime(mid), current_lon)

        if current_lon < target_lon:
            lower = mid
        else:
            upper = mid

    # print('Second phase', ephem_to_datetime(lower), ephem_to_datetime(mid), ephem_to_datetime(upper))
    jieqi_index = (jieqi_index-9)%24 # offset
    jieqi = JIEQI[jieqi_index]
    jieqi_time = ephem_to_datetime(upper)
    return jieqi, jieqi_time, jieqi_index


"""计算指定年份的二十四节气时间"""
def jq_new_hlon(date, observer, order):
    # 获取上一个和下一个节气
    jieqi = None
    jieqi_time = None
    jieqi_index = 0

    if order == -1:
        jieqi, jieqi_time, jieqi_index = find_previous_jieqi_hlon(date, observer)
        # print('prev', prev_jieqi, prev_jieqi_time)
    else:
        jieqi, jieqi_time, jieqi_index = find_next_jieqi_hlon(date, observer)
        # print('next', next_jieqi, next_jieqi_time)

    return jieqi, jieqi_time, jieqi_index


# 时间转换
# 获取当前环境的时区偏移值
def cur_time_zone_offset() -> float:
    current_time = datetime.now()
    # 获取UTC时间
    utc_time = datetime.utcnow()
    # 计算小时差
    time_zone_offset = round((current_time - utc_time).seconds / 3600.0, 2)
    # 以0.5为基准进行舍入操作:
    if time_zone_offset >= 0:
        time_zone_offset = math.floor(time_zone_offset / 0.5 + 0.5) * 0.5
    else: 
        time_zone_offset = math.ceil(time_zone_offset / 0.5 - 0.5) * 0.5
    
    return time_zone_offset


# 计算指定本地时间的真太阳时
def calculate_true_solar_time(
    date_local: datetime, 
    observer_latitude: Union[str, float]=0, 
    observer_longitude: Union[str, float]=0) -> datetime:
    """
    根据观测者的经纬度和时间, 计算真太阳时.
    参数:
    date_local (datetime): 观察时间
    observer_latitude (str): 观测者的纬度, 以度为单位.
    observer_longitude (str): 观测者的经度, 以度为单位.
    返回:
    datetime: 调整后的真太阳时, 已考虑观测者的经纬度和时区偏移量.
    """

    # 根据观测者的经纬度和时间, 计算真太阳时:
    # 计算当前时区偏移time_zone_offset
    time_zone_offset = cur_time_zone_offset()    
    # 将观测者的纬度和经度转换为ephem格式
    lat, lon = ephem.degrees(str(observer_latitude)), ephem.degrees(str(observer_longitude))
 
    # 创建一个新的观测者对象
    observer = ephem.Observer()
    observer.lat, observer.lon = lat, lon
 
    # 将输入的当地时间转换为ephem的Date格式
    date_0 = date_local - timedelta(hours=time_zone_offset)# 将本地时间转换为UTC+0时间
    utc_date = ephem.Date(date_0)
 
    # 设置观测者的日期和时间
    observer.date = utc_date
 
    # 计算太阳的位置
    sun = ephem.Sun(observer)
    # 计算本地太阳正午时间（在给定日期太阳通过当地子午线的时间）
    next_noon = observer.next_transit(sun).datetime() + timedelta(hours=time_zone_offset)
    # 计算钟表正午时间
    clock_noon = next_noon.replace(hour=12, minute=0, second=0, microsecond=0)
    # 计算太阳正午和钟表正午之间的差值
    time_diff = (clock_noon - next_noon).total_seconds() / 3600.0
    # 调整本地时间以获得真太阳时(ephem通常已经考虑光行差)
    true_solar_time = date_local + timedelta(hours=time_diff)
    return true_solar_time


# 将ephem时间转成datetime时间, 注意都是UT时间, ephem.localtime函数会出问题
def ephem_to_datetime(date_ephem):
    # ephem的基准日期（1899-12-31 12:00:00 UTC）
    EPOCH = datetime(1899, 12, 31, 12, 0, 0)

    delta = timedelta(days=date_ephem)
    date_ut = EPOCH + delta

    return date_ut


date_str = '2025-2-1 0:35:30'
observer_latitude, observer_longitude = '39.9', '110.8333'

date_test = datetime.strptime(date_str, '%Y-%m-%d %H:%M:%S')

observer = ephem.Observer()
observer.lat = observer_latitude
observer.lon = observer_longitude

true_solar_time = calculate_true_solar_time(date_test, observer_latitude, observer_longitude)
date_ut = true_solar_time - timedelta(hours=8)

jq_prev, jqdate_prev, jq_index = jq_new_hlon(date_ut, observer, -1) # 需要UT时间, 不是本地时间
jq_next, jqdate_next, jq_index = jq_new_hlon(date_ut, observer, 1) # 需要UT时间, 不是本地时间

jqdate_prev += timedelta(hours=8)
jqdate_next += timedelta(hours=8)

print(date_test.strftime('%Y-%m-%d %H:%M:%S'), '前一个节是')
print(format_with_padding(jq_prev, 'l', 10)+format_with_padding(jqdate_prev.strftime('%Y-%m-%d %H:%M:%S'), 'l', 10))
print(date_test.strftime('%Y-%m-%d %H:%M:%S'), '后一个节是')
print(format_with_padding(jq_next, 'l', 10)+format_with_padding(jqdate_next.strftime('%Y-%m-%d %H:%M:%S'), 'l', 10))