使用Python计算Algol视星等：代码无食变现象的问题排查

我尝试用Python计算恒星Algol的视星等，维基百科显示Algol的星等通常稳定在2.1左右，但每2.86天会因约10小时的偏食降至3.4。但我的代码运行后结果始终在2.2左右，完全没出现降至3.4的食变现象，请问代码哪里错了？

我的代码 #

from skyfield.api import Star, load
from skyfield.data import hipparcos
from datetime import timedelta
import math

# Load the JPL ephemeris DE421 (covers 1900-2050).
planets = load('de421.bsp')
earth = planets['earth']

# load Hipparcos ephemeris
with load.open(hipparcos.URL) as f:
    df = hipparcos.load_dataframe(f)

# Create a timescale
ts = load.timescale()
now = ts.now()

# load Algol from Hipparcos
algol = Star.from_dataframe(df.loc[14576])


absolute_magnitude = -0.07
        
for x in range(100):
        time = now + timedelta(hours=x)
        astrometric = earth.at(time).observe(algol)

        ra, dec, distance = astrometric.radec(epoch=time)

        apparent_magnitude = absolute_magnitude + (5 * math.log10((distance.au/206264.80749673))) - 5
        print(apparent_magnitude)

输出结果 #

2.20099154236185
2.200991542408758
2.2009915424554958
2.2009915425020647
etc.

错误原因分析 #

• 核心逻辑错误：你用的是基于地球到恒星距离的视星等计算公式，但Algol是食双星系统，它的亮度变化来自两颗恒星互相遮挡，和地球到它的距离无关。短时间内地球到Algol的距离变化极小，所以你的代码只会算出它的平均视星等，完全无法体现食变带来的亮度衰减。
• 公式误用：你写的视星等公式是距离模数公式，用来计算单颗恒星因距离变化产生的视星等差异，但Algol的距离在观测周期内几乎固定，所以结果只会在2.2左右小幅波动，和食变现象无关。

解决方案 #

要计算Algol的食变视星等，需要结合它的光变曲线参数（食周期、食持续时间、亮度变化幅度等）来判断观测时间是否处于食期，进而调整视星等数值。以下是一个简单的模拟示例：

from skyfield.api import load
from datetime import timedelta

ts = load.timescale()
now = ts.now()

# Algol的光变核心参数
average_mag = 2.1
minimum_mag = 3.4
period_days = 2.867329  # 精确食周期
eclipse_hours = 10       # 偏食持续时长

# 注意：实际需要从天文数据库获取精确的最近一次食始时间，这里用示例值
last_eclipse_start = ts.utc(2024, 5, 20, 18, 30, 0)

for x in range(100):
    current_time = now + timedelta(hours=x)
    # 计算当前时间距离上次食始的天数
    delta_days = (current_time - last_eclipse_start).total_seconds() / 86400
    # 取模得到当前在食周期内的位置
    phase_in_period = delta_days % period_days
    
    current_mag = average_mag
    # 判断是否处于食期（包括食前过渡、食中、食后过渡）
    half_eclipse_days = eclipse_hours / 24
    if phase_in_period <= half_eclipse_days:
        # 食始到食甚的亮度衰减
        phase = phase_in_period / half_eclipse_days
        current_mag = average_mag + (minimum_mag - average_mag) * phase
    elif phase_in_period >= (period_days - half_eclipse_days):
        # 食甚到食终的亮度恢复
        phase = (period_days - phase_in_period) / half_eclipse_days
        current_mag = average_mag + (minimum_mag - average_mag) * phase
    
    print(current_mag)

说明 #

这个示例是简化的线性模拟，实际Algol的光变曲线是不对称的，想要更精确的结果，需要从专业天文数据库获取它的光变曲线参数，包括食分、不同阶段的亮度变化率等，再优化计算逻辑。

内容的提问来源于stack exchange，提问作者velkyvont