Python で日食計算

Python で日食計算 (Skyfield版)

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

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

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

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

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

ts = load.timescale()
t = ts.utc(2020, 6, 21, 7, 0, range(0, 11000))
tz = timezone('Asia/Tokyo')

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

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

観測者の位置を指定する Topos オブジェクトを設定します。

osaka = earth + Topos('34.6914 N', '135.4917 E')

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

sun_app = osaka.at(t).observe(sun).apparent()
moon_app = osaka.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)

次に太陽と月の角距離を計算します。角距離は、太陽・月それぞれの見かけの位置をもとに、separation_from メソッドを使うことで、簡単に計算できます。
結果は配列で出力されます。

app_sep = sun_app.separation_from(moon_app).radians

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

percent_eclipse = (sun_rad + moon_rad - app_sep) / (sun_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

計算結果

食の最大: 2020-06-21 17:10:19 JST
最大食分: 0.54
欠け始め: 2020-06-21 16:06:24 JST
食の終わり: 2020-06-21 18:07:39 JST

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

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

# 初期時刻設定
ts = load.timescale()
t = ts.utc(2020, 6, 21, 7, 0, range(0, 11000))
tz = timezone('Asia/Tokyo')

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

# 観測地設定
osaka = earth + Topos('34.6914 N', '135.4917 E')

# 太陽・月の位置計算
sun_app = osaka.at(t).observe(sun).apparent()
moon_app = osaka.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)

# 太陽・月の角距離の計算
app_sep = sun_app.separation_from(moon_app).radians

# 食分の計算
percent_eclipse = (sun_rad + moon_rad - app_sep) / (sun_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版)

　江越　航のホームページ