大阪市立科学館 学芸員

 江越 航のホームページ

Python で日食計算

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.6883'
osaka.lon = '135.4933'
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_moon + r_sun) > s ・・・太陽と月の半径の和よりも、太陽と月の角距離が小さい)なら、太陽が欠けています。

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

計算結果(時刻、角距離、食分)を配列 results に追加します。

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

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

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

print(u"食の最大: %s" % max_eclipse[0])
print(u"最大食分: %.2f" % max_eclipse[2])

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

eclipse = False

for x in results:
    if x[2] > 0 :
        if eclipse == False:
            print(u"欠け始め: %s" % x[0])
            eclipse = True
    else :
        if eclipse == True:
            print(u"食の終わり: %s" % x[0])
            eclipse = False

計算結果

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

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

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.6883'
osaka.lon = '135.4933'
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), s, percent_eclipse])

# 食の最大の検索	
max_eclipse = max(results, key = itemgetter(2))
    
print(u"食の最大: %s" % max_eclipse[0])
print(u"最大食分: %.2f" % max_eclipse[2])

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

for x in results:
    if x[2] > 0 :
        if eclipse == False:
            print(u"欠け始め: %s" % x[0])
            eclipse = True
    else :
        if eclipse == True:
            print(u"食の終わり: %s" % x[0])
            eclipse = False