Python で日食計算

Python で日食計算 (PyEphem版)

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

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

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

import ephem
import datetime
import numpy as np
from operator import itemgetter

計算を開始する日時を設定します。日時は、世界時で指定します。

initial_date = datetime.datetime(2019, 12, 26, 5, 0, 0, 0, tzinfo=datetime.timezone.utc)

観測者のインスタンス（ここでは osaka）を作成します。
作成した観測者 osaka に、計算を始める日時を設定します。

osaka = ephem.Observer()
osaka.lat = '34.6914'
osaka.lon = '135.4917'
osaka.date = initial_date

太陽と月のインスタンス（ここではsun, moon）を作成します。

sun = ephem.Sun()
moon = ephem.Moon()

結果を入れる配列を初期化します。

results = []

これより初期時刻から1秒ずつ時間を進めながら、太陽と月の位置、および半径の計算を繰り返します。
日食の継続時間は概ね3時間程度のため、ここでは11000秒分、計算します。
太陽と月の大きさ（sun.size, moon.size）の単位は、角度の秒です（直径およそ1800秒=0.5度）。

for x in range(0, 11000):
    osaka.date = initial_date + datetime.timedelta(seconds = x)
    sun.compute(osaka)
    moon.compute(osaka)
    r_sun = sun.size/2
    r_moon = moon.size/2

次に太陽と月の角距離を計算します。角距離は、太陽・月それぞれの方位・高度の値をもとに、separation メソッドを使うことで、簡単に計算できます。
separation メソッドで計算される角度の単位は度なので、秒の単位にするために3600を掛けています。

なお赤経・赤緯の値を用いないのは、月や太陽の赤経・赤緯の値は地球中心から見た値になっていることから、地平視差が生じるためです。

    s = np.rad2deg(ephem.separation((sun.az, sun.alt), (moon.az, moon.alt))) * 60 * 60

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

    percent_eclipse = (r_sun + r_moon - s) / (r_sun * 2)

計算結果（時刻、食分）を配列 results に追加します。

    results.append([ephem.localtime(osaka.date), percent_eclipse])

食分の最大値と、その時刻を探します。
配列 results の1番目の要素の最大値が、食分の最大値です。

max_eclipse = max(results, key = itemgetter(1))

print('食の最大:', max_eclipse[0])
print('最大食分: {0:.2f}'.format(max_eclipse[1]))

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

eclipse = False

for x in results:
    if x[1] > 0 :
        if eclipse == False:
            print('欠け始め:', x[0])
            eclipse = True
    else :
        if eclipse == True:
            print('食の終わり:', x[0])
            eclipse = False

計算結果

食の最大: 2019-12-26 15:31:37
最大食分: 0.37
欠け始め: 2019-12-26 14:22:19
食の終わり: 2019-12-26 16:33:08

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

import ephem
import datetime
import numpy as np
from operator import itemgetter

# 初期時刻設定
initial_date = datetime.datetime(2019, 12, 26, 5, 0, 0, 0, tzinfo=datetime.timezone.utc)

# 観測地設定
osaka = ephem.Observer()
osaka.lat = '34.6914'
osaka.lon = '135.4917'
osaka.date = initial_date

# 太陽・月
sun = ephem.Sun()
moon = ephem.Moon()

# 結果のリストを初期化
results = []

# 1秒ずつ3時間分　計算繰り返し
for x in range(0, 11000):
    # 時刻を1秒進める
    osaka.date = initial_date + datetime.timedelta(seconds = x)

    # 太陽・月の位置・半径計算
    sun.compute(osaka)
    moon.compute(osaka)
    r_sun = sun.size/2
    r_moon = moon.size/2

    # 太陽・月の角距離の計算
    s = np.rad2deg(ephem.separation((sun.az, sun.alt), (moon.az, moon.alt))) * 60 * 60

    # 食分の計算
    percent_eclipse = (r_sun + r_moon - s) / (r_sun * 2)
    
    # 計算結果を追加（時刻、食分）
    results.append([ephem.localtime(osaka.date), percent_eclipse])

# 食の最大の検索	
max_eclipse = max(results, key = itemgetter(1))
    
print('食の最大:', max_eclipse[0])
print('最大食分: {0:.2f}'.format(max_eclipse[1]))

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

for x in results:
    if x[1] > 0 :
        if eclipse == False:
            print('欠け始め:', x[0])
            eclipse = True
    else :
        if eclipse == True:
            print('食の終わり:', x[0])
            eclipse = False

江越 航のホームページ

Python で日食計算 (PyEphem版)

　江越　航のホームページ