Python で月食計算

Python で月食計算 (Skyfield版)

Skyfield を用いると、月食が起こる時刻を計算することができます。
（なお、月食が起こる概略の日時は既知のものとします。）

参考にしたページ：
Eclipse Calculations using Python
GitHub - Various scripts to calculate eclipses

必要モジュールをインポートします。

from skyfield.api import load
from pytz import timezone
import numpy as np

Timescale オブジェクトを読み込み、計算する日時を設定します。日時は、世界時で指定します。
時刻は配列で指定することが可能です。以下、時刻を指定した計算値は、配列で出力されます。
月食の継続時間は最大でも4時間程度のため、ここでは15000秒分、計算します。
また、タイムゾーンを日本に設定します。

ts = load.timescale()
t = ts.utc(2021, 5, 26, 9, 0, range(0, 15000))
tz = timezone('Asia/Tokyo')

天体暦のファイルを読み込み、太陽、月、地球のオブジェクト（ここではsun, moon, earth）を作成します。

eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']

地球中心から見た太陽と月の位置を計算します。
時刻に対応して、結果は配列で出力されます。

sun_app = earth.at(t).observe(sun).apparent()
moon_app = earth.at(t).observe(moon).apparent()

太陽・月の半径、距離（ここではkm単位）から、見かけの大きさを計算します。
結果は配列で出力されます。

r_sun = 696000
sun_dist = sun_app.distance().km
sun_rad = np.arctan2(r_sun, sun_dist)

r_moon = 1737
moon_dist = moon_app.distance().km
moon_rad = np.arctan2(r_moon, moon_dist)

次に太陽と月の視差を計算します。
視差は (地球半径)/(太陽または月までの距離) で計算できます。
結果は配列で出力されます。

r_earth = 6378
parallax_sun = r_earth / sun_dist
parallax_moon = r_earth / moon_dist

視差から、地球の本影の視半径を求めます。
本影の視半径は (太陽の視差)+(月の視差)-(太陽の視半径) で計算できます。（下図）
地球の本影の視半径の大きさは、大気による屈折で約1/50大きくなることから、51/50 を掛けています。
結果は配列で出力されます。

umbra = (parallax_moon - sun_rad + parallax_sun) * 51/50

次に月と地球の本影の角距離を計算します。
本影は太陽の反対側になるので、太陽と月の角距離を求めて180度を引きます。
角距離は separation_from メソッドを使うことで、簡単に計算できます。
結果は配列で出力されます。

app_sep = abs(sun_app.separation_from(moon_app).radians - np.deg2rad(180))

食分を計算します。結果は配列で出力されます。
食分がプラス（(umbra + r_moon) > s ・・・地球の本影と月の半径の和よりも、本影中心と月の角距離が小さい）なら、月が欠けています。

percent_eclipse = (umbra + moon_rad - app_sep) / (moon_rad * 2)

食分が最大となる、配列のインデックスを探します。

max_i = np.argmax(percent_eclipse)

食分が最大となる時刻と食分値を出力します。

print('食の最大:', t[max_i].astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
print('最大食分: {0:.2f}'.format(percent_eclipse[max_i]))

欠け始めと食の終わりの時刻を探します。
最初に食分が0以上になった時刻が欠け始め（第一接触）、再び食分が0になった時刻が食の終わり（第四接触）です。

eclipse = False

for ti, pi in zip(t, percent_eclipse):
    if pi > 0 :
        if eclipse == False:
            print('欠け始め:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
            eclipse = True
    else :
        if eclipse == True:
            print('食の終わり:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
            eclipse = False

計算結果

食の最大: 2021-05-26 20:18:43 JST
最大食分: 1.02
欠け始め: 2021-05-26 18:44:28 JST
食の終わり: 2021-05-26 21:52:59 JST

全体は次のようになります。

from skyfield.api import load
from pytz import timezone
import numpy as np

# 初期時刻設定
ts = load.timescale()
t = ts.utc(2021, 5, 26, 9, 0, range(0, 15000))
tz = timezone('Asia/Tokyo')

# 太陽・月・地球
eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']

# 太陽・月の位置計算
sun_app = earth.at(t).observe(sun).apparent()
moon_app = earth.at(t).observe(moon).apparent()

# 太陽・月の見かけの大きさ計算
r_sun = 696000
sun_dist = sun_app.distance().km
sun_rad = np.arctan2(r_sun, sun_dist)

r_moon = 1737
moon_dist = moon_app.distance().km
moon_rad = np.arctan2(r_moon, moon_dist)

# 視差・本影の視半径計算
r_earth = 6378
parallax_sun = r_earth / sun_dist
parallax_moon = r_earth / moon_dist
umbra = (parallax_moon - sun_rad + parallax_sun) * 51/50

# 月・地球の本影の角距離の計算
app_sep = abs(sun_app.separation_from(moon_app).radians - np.deg2rad(180))

# 食分の計算
percent_eclipse = (umbra + moon_rad - app_sep) / (moon_rad * 2)
    
# 食の最大の検索	
max_i = np.argmax(percent_eclipse)
    
print('食の最大:', t[max_i].astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
print('最大食分: {0:.2f}'.format(percent_eclipse[max_i]))

# 欠け始めと食の終わりの検索
eclipse = False

for ti, pi in zip(t, percent_eclipse):
    if pi > 0 :
        if eclipse == False:
            print('欠け始め:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
            eclipse = True
    else :
        if eclipse == True:
            print('食の終わり:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
            eclipse = False

江越 航のホームページ

Python で月食計算 (Skyfield版)

　江越　航のホームページ