自然リリス計算法メモ

自然リリス

NotebookLMの記事より覚書

SkyfieldとSwiss Ephemerisの主な計算精度の違い（不一致の原因）は、主に以下の5つの要素に起因します。両者は基礎となる設計思想が異なるため、同じ計算結果を得るためには設定を厳密に合わせる必要があります。

1. エフェメリス（天体暦）カーネルと圧縮モデルの違い
Skyfieldは、NASA JPLが提供する生の高精度なSPKファイル（DE421やDE440など）を直接読み込んで計算を行います。一方、Swiss Ephemerisは独自の内蔵モデル（Moshierモデル）や、データをビット圧縮したファイル（DE431ベースなど）を使用することが多いです。この基礎データの違いや圧縮処理により、月などの位置で数秒角のベースライン誤差が生じることがあります。

2. 座標系の基準と変換（ICRSと日付の黄道）
Skyfieldの中核エンジンは、恒星を基準とした固定の赤道座標系である「ICRS（国際天球基準座標系）」で動作します。対してSwiss Ephemerisは、占星術で標準的な「日付の地心黄道座標（Geocentric Ecliptic of Date）」を優先します。
両者の結果を一致させるには、Skyfield側で歳差（Precession）と章動（Nutation）の両方を適用した「真の黄道と春分点（True Ecliptic and Equinox of Date）」の座標フレームを明示的に指定する必要があります。章動を無視した「平均（Mean）」のフレームを使ってしまうと、数秒角の誤差が発生します。

3. 時間基準と$\Delta T$（地球の自転の遅れ）の補正
天文学における精度は、UTC（協定世界時）からTT（地球時）への変換精度に依存します。地球の自転の遅れを補正する「$\Delta T$」について、Skyfieldは定期的に更新される内部テーブルを使用しますが、Swiss Ephemerisが異なる多項式モデルや古いテーブルを使っている場合、同じUTC時刻を指定しても内部的な計算時刻（TT）に数秒の差が生じます。月は1秒間に約0.5秒角移動するため、このわずかな時間のズレがノード（交点）やリリス（遠地点）の明確な位置誤差につながります。

4. 視位置（Apparent）と幾何学的位置（Geometric）の扱いの違い
月の真のノード（True Node）や真のリリス（True Lilith）は、軌道の瞬間的な数学的特性です。Swiss Ephemerisはこれらを計算する際、光が地球に届くまでの時間（光行差）を無視した「幾何学的位置（Geometric position）」をデフォルトで出力します。
もしSkyfieldで光行差を含む「視位置（Apparent position）」を計算してしまうと、光が月から地球に届くまでの約1.3秒分、要素が過去にシフトしてしまい、結果がずれてしまいます。

5. 観測地の視差（地心か局所的か）
Swiss Ephemerisは、ノードやリリスの計算において観測地の座標（Topocentricフラグ）を意図的に無視し、常に地球の中心（Geocentric）からの位置を出力します。Skyfieldで特定の緯度・経度からの見え方を計算してしまうと、視差が加わってしまいSwiss Ephemerisの標準出力と一致しなくなります。

まとめ
Swiss Ephemerisと同等の精度や結果をSkyfieldで再現するには、「同じエフェメリスを使用する」「時間をTT（地球時）で揃える」「章動を含んだ日付の黄道座標系を指定する」「光行差を含まない地球中心（幾何学的）のベクトルを用いる」という厳密なパラメータの同期が求められます。

SkyfieldでSwiss Ephemerisと同じ章動（Nutation）モデルを適用し、章動と歳差の両方を含んだ座標系を再現するには、Skyfieldの framelib モジュールに用意されている ecliptic_frame を使用します。

Swiss Ephemerisのデフォルト出力は「日付の真の地心黄道座標（True Geocentric Ecliptic of Date）」であり、これをSkyfieldで正確に一致させる手順は以下の通りです。

1. ecliptic_frame の指定
skyfield.framelib の ecliptic_frame は、章動と歳差の両方を適用した「日付の真の黄道と春分点」の座標フレームです。章動を無視した「平均（Mean）」のフレームや、固定のJ2000黄道フレームを使ってしまうと、数秒角の誤差が発生してしまいます。

2. 軌道要素の計算での適用
真のノードやリリス（月の遠地点）などの接触軌道要素（Osculating elements）を計算する場合は、osculating_elements_of 関数の reference_frame 引数にこの ecliptic_frame を渡します。

from skyfield.framelib import ecliptic_frame
from skyfield.elementslib import osculating_elements_of

# 軌道要素を計算する際にecliptic_frameを指定
elements = osculating_elements_of(geometric_position, reference_frame=ecliptic_frame)

3. 黄道座標（黄経・黄緯）の取得
単純に天体の位置を黄道座標（黄経・黄緯）で取得したい場合は、計算した位置に対して frame_latlon を使用します。

from skyfield.framelib import ecliptic_frame

# 位置ベクトルを黄道座標の緯度・経度に変換
lat, lon, distance = geometric_position.frame_latlon(ecliptic_frame)

注意点（計算の前提）
Swiss Ephemerisの出力と完全に一致させるためには、座標フレームを合わせるだけでなく、以下の点も満たす必要があります。

• 位置ベクトルを求める際、光行差を含まない「幾何学的位置（Geometric position）」を使用する（.apparent() などの視位置変換は行わない）。
• 天体のベクトルは、太陽系重心（Barycentric）ではなく、地球の中心から対象天体へのベクトル（Geocentric）を自身で引き算して求める。
• 地球の自転の遅れを補正するため、計算する時刻をTT（地球時）で厳密に揃える。

Skyfieldで光行差（光が地球に届くまでの遅れ）を無視した「幾何学的位置（Geometric position）」を計算するには、対象となる天体の位置ベクトルと地球の位置ベクトルを、同一時刻において直接引き算します。

Skyfieldで用意されている .observe() メソッドを使用すると、光が天体から地球に届くまでの時間（光行差）が自動的に計算され、過去の位置（Astrometric position）が返されてしまいます。また、そこからさらに .apparent() を呼び出すと光の湾曲なども含まれた視位置（Apparent position）になってしまいます。

そのため、光行差を完全に無視して瞬間的な数学的座標（幾何学的位置）を取得するには、ベクトル同士の演算のみを用います。以下がそのコード例です。

from skyfield.api import load

# 1. エフェメリス（天体暦）と対象天体をロード
# ※Swiss Ephemerisと一致させる場合はDE431やDE441などを利用します
eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']

# 2. タイムスケールの作成と時刻の指定
# ※地球の自転の遅れ（ΔT）によるズレをなくすため TT（地球時）を指定します
ts = load.timescale()
t = ts.tt(2023, 1, 1, 0, 0)

# 3. 幾何学的位置（Geometric position）の計算
# observe() を使わずに、対象天体から地球を引き算した関数に対して .at(t) を実行します
geometric_position = (moon - earth).at(t)

# ---------------------------------------------------------
# 【参考】以下のように各天体の位置を計算してから引き算しても同じ結果になります
# earth_pos = earth.at(t)
# moon_pos = moon.at(t)
# geometric_position = moon_pos - earth_pos
# ---------------------------------------------------------

# この geometric_position に対して、前回解説した ecliptic_frame などを渡して
# 軌道要素や黄道座標を計算します。（ .apparent() は絶対に呼び出さないでください ）

このコードのポイント:

• (moon - earth) のように天体（ベクトル関数）同士を引き算することで、地球中心（Geocentric）から対象天体へのベクトルを定義できます。
• これを特定の時間 t で評価する（.at(t)）ことで、同一時刻における数学的な相対位置が得られます。
• 得られた geometric_position をそのまま osculating_elements_of() などに渡すことで、Swiss Ephemerisと同様に光行差を含まない真のノードやリリス（月の遠地点）などを算出することができます。

Skyfieldで計算した光行差を含まない「幾何学的位置」から、Swiss Ephemerisと同等の精度で月の真のノード（True Node）や真のリリス（True Lilith / 接触遠地点）を計算するには、osculating_elements_of 関数を使用します。

計算を正確に行うためには、「地球と月の重力パラメータ（GM）の和」を指定することと、リリスの位置を近点（Periapsis）の反対側として「180度シフト」させる計算が必要です。

以下が具体的なPythonコード例です。

from skyfield.api import load
from skyfield.framelib import ecliptic_frame
from skyfield.elementslib import osculating_elements_of
from skyfield.data.gravitational_parameters import GM_dict

# 1. エフェメリスのロード
eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']

# 2. 地球と月の重力パラメータ（GM）の和を算出
# ※Swiss Ephemerisと一致させるために必須の処理です
mu_earth = GM_dict  # 地球のGM
mu_moon = GM_dict   # 月のGM
mu_sys = mu_earth + mu_moon

# 3. 日時の設定（TT: 地球時を使用）
ts = load.timescale()
t = ts.tt(2023, 1, 1, 0, 0)

# 4. 幾何学的位置（Geometric position）の取得
geometric_position = (moon - earth).at(t)

# 5. 接触軌道要素（Osculating elements）の計算
# 座標系に ecliptic_frame（日付の真の黄道）を指定し、重力パラメータに mu_sys を渡す
elements = osculating_elements_of(
    geometric_position,
    reference_frame=ecliptic_frame,
    gm_km3_s2=mu_sys
)

# ---------------------------------------------------------
# 【真のノード（True Node）の計算】
# 昇交点黄経（Longitude of Ascending Node）をそのまま取得します
true_node = elements.longitude_of_ascending_node.degrees
print(f"真のノード (True Node): {true_node:.4f}度")

# ---------------------------------------------------------
# 【真のリリス（True Lilith / 遠地点）の計算】
# Skyfieldは近点引数（Periapsis）を出力するため、遠地点（Apogee）を求めるには
# 「昇交点黄経 ＋ 近点引数 ＋ 180度」を計算し、360度で丸めます
omega = elements.longitude_of_ascending_node.degrees
omega_peri = elements.argument_of_periapsis.degrees

true_lilith = (omega + omega_peri + 180.0) % 360.0
print(f"真のリリス (True Lilith): {true_lilith:.4f}度")

重要なポイント

1. 重力パラメータ（GM）の加算 (mu_sys)
   高精度な接触軌道要素を求める場合、中心天体の重力だけでなく、対象天体の重力も考慮する必要があります。地球・月系においてSwiss Ephemerisの結果を再現するには、gm_km3_s2 引数に地球と月のGMの合計値を明示的に渡す必要があります。Skyfieldのデフォルトでは中心天体（地球）の質量しか考慮されない場合があるため、この設定が不可欠です。
2. 真のリリスは「遠地点（Apogee）」である
   elementslib モジュールは、軌道上の最も近い点である「近点（Periapsis）」の角度（近点引数）を出力します。真のリリスは軌道上で最も遠い「遠地点」を指すため、昇交点黄経と近点引数を足し合わせた近点の黄経に対し、さらに180度を加算して反対側に向ける必要があります。これを忘れると、結果がリリスではなくプリアポス（Priapus / 近地点）になってしまい、ちょうど180度のズレが生じます。
3. 視位置変換（.apparent()）の絶対回避
   軌道要素の計算元となるベクトルには、光行差を含んだ視位置ではなく必ず相対引き算のみで求めた幾何学的位置（geometric_position）を使用してください。光行差が含まれていると、光が到達するまでの約1.3秒分だけ過去の軌道要素が計算されてしまい、Swiss Ephemerisの標準結果と一致しなくなります。

ノードが「月の軌道と黄道面が交差する点」であり、日蝕や月蝕という視覚的な天文現象によって位置を確認できるのに対し、リリス（月の遠地点）は空間に存在する数学的な点であるため、光を放つこともなく視覚的に直接観測することはできません。

そのため、天文観測的にリリス（遠地点）の正しい位置を判断する基準は、「地球の中心から月の中心までの物理的な距離が最大（極大）になる瞬間とその位置」を測定することにあります。

具体的には、以下のようなアプローチでその正しさが検証・判断されています。

1. 物理的な最大距離の特定（Actual / Measured Apogee）
天文学（NASAのJPL Horizonsなどの高精度なデータ）を用いて、月が地球から最も遠ざかる（距離が極大値から減少に転じる）正確な時刻を秒単位で割り出します。この「距離が最大になるというイベント（瞬間）」において月が実際に位置している黄経が、観測的に正しい真の遠地点（Actual Apogee または Measured Apogee）とみなされます。

2. 各種計算モデルと実際の最大距離地点との照合
占星術ソフトウェアやSwiss Ephemerisなどの天文計算エンジンでは、リリスの計算に主に3つのモデルが用いられますが、これらが「実際に月が最も遠ざかるタイミングと位置」にどれだけ一致するかが、正確さの指標となります。

• 平均リリス（Mean Lilith）: 軌道の変動を平均化したものであり、完全な円運動を前提としているため、月が実際に最も遠ざかるタイミングや位置とは一致しません。
• 真のリリス（True / Osculating Lilith）: 瞬間的な軌道要素（接触軌道）から算出されるこの点は、月が実際に遠地点に到達するタイミングと位置に非常に正確に（誤差1秒や数秒角以内で）一致することが確認されています。
• 自然/補間リリス（Natural / Interpolated Lilith）: 月が実際に最大距離に達するポイントそのものを捉えるために軌道を補間して導き出されたモデルです。実際の最大距離のポイントをOsculatingモデル以上に正確に（実際の極大地点から数分角以内の誤差で）トレースできることが検証されています。

3. 「視位置」の概念がない（光行差を含まない）
リリスは実体を持たない空間上の点であり、そこから光が発せられるわけではありません。したがって、光が地球に届くまでの遅れ（光行差）を含んだ「視位置（Apparent position）」で観測的な正しさを定義することはできず、光行差を無視した「幾何学的な真の位置（True position）」同士で月の位置とリリスの位置を比較しなければならないという厳密なルールがあります。

このように、リリスの天文学的な正しさは「目に見える現象」ではなく、「地球と月の物理的距離の極大値」という定量的な観測データを基準にした照合によって判断されています。

Skyfieldには、Swiss Ephemerisのように「自然リリス（Interpolated Lilith）」を直接算出するための組み込み関数や機能（Swiss EphemerisにおけるSE_INTP_APOGのようなもの）は用意されていません。

自然リリス（Interpolated）は、Swiss Ephemerisの開発者であるDieter Kochが独自のヒューリスティクス（発見的手法）を用いて考案・プログラムした特殊なモデルです。Swiss Ephemerisはこうした複雑な計算を単一の呼び出しで実行できる専門的なライブラリですが、Skyfieldは汎用的な天文計算のモジュール式ツールキットであるため、このような占星術特有のアルゴリズムは標準実装されていません。

Skyfieldで同等の概念（実際の遠地点）を算出する代替アプローチ
自然リリスは、純粋な楕円軌道ではなく「現在の月のサイクルにおいて、実際に月が地球から最も遠ざかる（最大距離に達する）ポイント」を捉えようとするモデルです。Skyfieldでこの「実際の遠地点（Actual/Measured Apogee）」を求めるには、ユーザー自身で距離の最大値を検索するコードを記述する必要があります。

具体的には以下のアプローチで算出が可能です：

1. 距離の極大値（イベント）の検索: Skyfieldには、時間の関数が返す値の局所的な最大値を検索する機能があります。これを用いて、地球と月の距離が極大（増加から減少に転じる瞬間）となる正確な時刻を特定します。
2. その瞬間の位置の取得: 算出された「距離が最大になる瞬間」における月の幾何学的位置（光行差を含まない位置）を取得します。

結論として、関数を1つ呼び出すだけで自然リリスを取得することはできませんが、Skyfieldの関数を使って「地球と月の物理的な距離が最大になる瞬間」をプログラムで検索・特定することで、自然リリスが意図している観測的な遠地点を自力で算出することは可能です。

Skyfieldを使用して、地球と月の距離が物理的に最大（極大）になる瞬間を検索するには、skyfield.searchlib モジュールに用意されている find_maxima 関数を利用します。

これは、光行差を含まない「幾何学的位置」から地球と月との距離を計算する独自の関数を定義し、その関数が返す値の極大値（局所的な最大値）となる時間を探索するというアプローチです。この手法で特定された時刻と位置は、Swiss Ephemerisの議論でも「実際の遠地点（Measured Apogee / Actual Apogee）」や自然リリスの観測的な根拠として重視されているものです。

以下が、その瞬間を検索するための具体的なコード例です。

from skyfield.api import load
from skyfield.searchlib import find_maxima
from skyfield.framelib import ecliptic_frame

# 1. エフェメリスのロード
eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']

# 2. 地球と月の距離を時間の関数として定義する
def earth_moon_distance(t):
    # 光行差を無視した幾何学的位置（Geometric position）を計算
    geometric_position = (moon - earth).at(t)
    # 距離（キロメートル）を返す
    return geometric_position.distance().km

# 時間に依存して変化する関数としてSkyfieldに認識させるための設定
earth_moon_distance.step_days = 1.0  # 1日刻みで大まかに探索し、その後詳細に特定する

# 3. 検索期間の指定（例として2023年1月の1ヶ月間）
ts = load.timescale()
t0 = ts.tt(2023, 1, 1, 0, 0)
t1 = ts.tt(2023, 1, 31, 0, 0)

# 4. find_maxima で距離が極大となる瞬間を検索
t_max, max_distances = find_maxima(t0, t1, earth_moon_distance)

# 5. 検索結果の表示と、その瞬間の位置（実際の遠地点）の取得
for t, dist in zip(t_max, max_distances):
    print(f"距離が最大になる瞬間 (TT): {t.tt_strftime('%Y-%m-%d %H:%M:%S')}")
    print(f"最大距離: {dist:.3f} km")

    # その瞬間の幾何学的位置ベクトルを取得
    geometric_position_at_max = (moon - earth).at(t)

    # 日付の真の黄道（章動と歳差を含む）における座標を取得
    lat, lon, _ = geometric_position_at_max.frame_latlon(ecliptic_frame)
    print(f"実際の遠地点の黄経: {lon.degrees:.4f}度\n")

検索処理の重要なポイント

1. find_maxima の使用
   Skyfieldの find_maxima は、時間の関数が返す値の局所的な最大値を高精度で特定するための関数です。これを使用することで、数式上の計算だけでなく、実際の軌道データに基づいて月が地球から最も遠ざかる（距離の変化がプラスからマイナスに転じる）瞬間を秒単位の精度で割り出すことができます。
2. step_days の指定
   関数に step_days を設定することで、探索のステップ幅をSkyfieldに指示します。月は約27.3日で軌道を一周するため、ステップ幅を 1.0（1日）程度に設定しておけば、遠地点を通り過ぎて見逃すことなく確実に関数の極大値をキャッチできます。
3. 距離計算メソッド distance().km
   引き算で求めた幾何学的位置（(moon - earth).at(t)）に対して .distance().km を呼び出すことで、スカラー値である距離（キロメートル）を取得します。これを関数の戻り値とすることで、find_maxima が最大距離を評価できるようになります。
4. 真の黄道フレームの適用
   特定された最大距離の瞬間 t を用いて再度幾何学的位置を計算し、それを ecliptic_frame に渡すことで、Swiss Ephemeris等の出力系（日付の真の黄道）と一致する「実際の遠地点」の黄経を正確に算出することができます。

Skyfieldで算出した「実際の遠地点（物理的な最大距離）」と、Swiss Ephemerisの「自然リリス（Interpolated Apogee: SE_INTP_APOG）」を比較するには、特定した最大距離の時刻を用いて、両者の出力する黄経（Longitude）を照合します。

比較を行う際の手順と、結果に生じる避けられないズレ（誤差）についての重要な背景は以下の通りです。

1. 比較の具体的な手順

PythonなどでSwiss Ephemerisのラッパー（pyswissephなど）が使える環境であれば、前回Skyfieldで特定した「距離が最大になる時刻」と同じ時間を使って SE_INTP_APOG を計算します。

Swiss Ephemeris側の定数指定:

• 天体IDとして 21 (SE_INTP_APOG) を指定します。
• フラグには、光行差や視差を除外するために、真位置（幾何学的位置）を示す設定を指定します。

コードのイメージ:

import swisseph as swe

# Skyfieldで算出した「距離が最大になる瞬間」のUTC時刻（例: 2023年1月1日 12:00）
year, month, day, hour = 2023, 1, 1, 12.0

# ユリウス日の取得
jd = swe.julday(year, month, day, hour)

# SE_INTP_APOG (21) の位置を計算（フラグには真位置・黄道座標などを指定）
# ※ swissephのデフォルトは幾何学的地心座標です
result, flag = swe.calc_ut(jd, 21)

# Swiss Ephemerisの自然リリスの黄経
se_intp_apog_lon = result

print(f"Swiss Ephemeris SE_INTP_APOG: {se_intp_apog_lon:.4f}度")
# これと、Skyfieldで求めた実際の月の黄経（Actual Apogee）を比較する

2. 両者の結果は「完全には一致しない」という前提

ここで最も重要なのは、Skyfieldで算出した「実際の月が最も遠ざかる黄経」と、Swiss Ephemerisの「SE_INTP_APOGの黄経」を比較した場合、完全には一致せず、数分角のズレが生じるという事実です。

Swiss Ephemerisのグループでの議論や開発者の検証において、以下のことが明らかになっています。

• SE_INTP_APOG の計算モデルの正体: Swiss Ephemerisの開発者であるDieter Koch氏によって考案された自然リリス（Interpolated apogee）は、「白魔術や発見的手法（heuristics）」とも表現される独自の補間アルゴリズムによって導き出されています。
• 誤差の存在: SE_INTP_APOGは、真のリリス（Osculating Apogee: 2体問題を前提とした瞬間的な軌道要素）よりも、実際の月が最大距離に到達する位置をよく近似しますが、完全にトレースできるわけではありません。
• 具体的なズレの例: Swiss Ephemerisの開発者Alois Treindl氏が行った検証（2009年11月22日の遠地点）によると、月が物理的な最大距離に到達した瞬間、実際の月の黄経は 308.345度 でした。その同日同時刻の SE_INTP_APOG は 308.3236度 であり、真のリリス（Osculating: 308.3139度）よりは実際の月に近いものの、約0.02度（1.2分角程度）のズレが残ることが確認されています。

まとめ

比較を行う際は、Skyfieldで算出した結果（実際の最大距離における月の位置）こそが物理的・観測的な「真の遠地点」であり、Swiss Ephemerisの SE_INTP_APOG は、その物理的地点を数学的に近似しようと試みた「独自の補間モデル」であるという違いを認識しておく必要があります。

両者を比較することで、Swiss Ephemerisの補間アルゴリズムが、実際の天体の物理的運動（Skyfieldの出力）に対してどれくらいの誤差（数分角程度）を持っているかを測る検証を行うことができます。

Skyfieldのfind_maxima関数で計算の精度をさらに高めるには、関数の引数としてepsilon（およびnum）を設定する方法があります。

Skyfieldの検索ルーチンであるfind_maxima（およびfind_minima）には、探索の精度やイテレーション（反復）に関するパラメータとしてepsilonとnum引数が用意されています。

以下のようにepsilon引数を指定することで、デフォルトの設定よりも細かな時間分解能で極大値（最大距離となる瞬間）を特定することが可能です。

# 1ミリ秒（1/86400000 日）の精度を指定
# （※epsilonに日数を指定する仕様や、以下のミリ秒の算出方法はソース外の一般的なPython/Skyfieldの知識に基づいています）
epsilon_days = 1.0 / (24 * 60 * 60 * 1000)

# epsilon引数を追加して高精度に検索を実行
t_max, max_distances = find_maxima(
    t0,
    t1,
    earth_moon_distance,
    epsilon=epsilon_days
)

精度を高めるためのポイント

1. epsilonの指定:
   検索において「どれくらいの時間誤差を許容するか」を決定します。この値を小さくする（例：秒単位からミリ秒単位にする）ことで、Swiss Ephemerisのような他エンジンと秒以下のレベルで結果を照合する際の、時間的な丸め誤差を最小限に抑えることができます。
2. num引数の活用:
   find_maxima関数はnum引数も受け取ります。これは内部の探索処理における反復回数や分割数に関わるパラメータであり、必要に応じて調整できます。
3. step_daysとの違い:
   前回設定した earth_moon_distance.step_days = 1.0 は、大まかに谷や山を「見逃さない」ための最初のステップ幅の指定です。これに対し、epsilonは見つけた山の頂上（極大値）を「どれだけピンポイントに絞り込むか」という最終的な到達精度を決定します。

Swiss Ephemerisの高精度な出力結果（SE_INTP_APOGなど）と比較検証を行う場合は、このepsilonを小さく設定してSkyfield側の時間特定のブレをなくしておくことが重要です。

Skyfieldの find_maxima 関数で指定する epsilon 引数は、「日数（days）」を単位として指定する仕様になっています。1日は86,400秒であるため、ミリ秒以下の精度（例えば1マイクロ秒 = 0.001ミリ秒）を設定するには、1日をマイクロ秒単位に換算して割り算を行います。

具体的なコードは以下のようになります。

from skyfield.api import load
from skyfield.searchlib import find_maxima

# --- 精度の計算 ---
# 1日 = 24時間 * 60分 * 60秒 = 86400秒
# 1ミリ秒 = 1.0 / (86400 * 1000) 日
# 1マイクロ秒（0.001ミリ秒） = 1.0 / (86400 * 1000 * 1000) 日

# epsilonを1マイクロ秒（ミリ秒以下）に設定
epsilon_microsec = 1.0 / (86400 * 1000 * 1000)

# --- 検索の実行 ---
# t0 から t1 の間で、earth_moon_distance が極大になる瞬間を高精度に探索
t_max, max_distances = find_maxima(
    t0,
    t1,
    earth_moon_distance,
    epsilon=epsilon_microsec
)

# 結果の確認
for t, dist in zip(t_max, max_distances):
    # マイクロ秒まで表示して精度を確認
    print(f"最大距離の瞬間 (TT): {t.tt_strftime('%Y-%m-%d %H:%M:%S.%f')}")
    print(f"最大距離: {dist:.6f} km")

ポイント:

• .%f によるマイクロ秒の表示: t.tt_strftime やその他の時間フォーマット関数において %f を使用することで、マイクロ秒の精度まで結果を出力・確認することができます。
• 計算負荷とのトレードオフ: epsilon を極端に小さくすると（例：ナノ秒レベル）、内部アルゴリズムがその精度に到達するまで反復（イテレーション）を繰り返すため、計算時間が長くなる可能性があります。Swiss Ephemeris（SE_INTP_APOG など）との比較においては、1マイクロ秒〜100マイクロ秒程度の設定で十分な精度が確保できます。

申し訳ありませんが、こちらのチャット環境ではPythonコードを直接実行することができません。

しかし、ご提供いただいたソース資料の中には、天文学の研究者やSwiss Ephemerisの開発者たちが、同等の高精度アルゴリズム（NASAのJPL Horizonsデータなど）を用いて計算した「実際の遠地点（最大距離の瞬間）」の実行結果・検証データが含まれています。

ご自身の環境でコードを実行した際に、正しく高精度な計算ができているかを確認するための「正解データ（ベンチマーク）」として、以下のデータをご活用いただけます。

過去の検証データから見る実際の遠地点の計算結果

事例1: 2021年5月11日の遠地点

• 最大距離到達の瞬間 (UTC): 2021年5月11日 22:09:34
• 最大距離: 406,513.4585317285 km
  ※ NASAのJPL Horizonsを用いて小数点以下8桁（ミリメートル以下の精度）まで計算された物理的な最大距離です。

事例2: 2021年6月8日の遠地点

• 最大距離到達の瞬間 (UTC): 2021年6月8日 02:39:21.8
• 最大距離: 406,211.04826419445 km
  ※ サブ秒（100ミリ秒単位）の時刻精度と高い距離精度で特定されています。

事例3: 2009年11月22日の遠地点

• 最大距離到達の瞬間 (UT): 2009年11月22日 20:10:00 付近
• 最大距離: 0.00270572328246 AU (約404,770 km)
• 実際の月の黄経: 308.3457度
  ※ Swiss Ephemeris開発者のAlois Treindl氏が、反復計算から実際の最大距離の瞬間を特定した際のデータです。

お手元でコードを実行し検証する際のアドバイス

前回ご提示した「1マイクロ秒（epsilon=1.0 / (86400 * 1000 * 1000)）」の精度を指定したSkyfieldコードをお手元のPC等で実行される場合、検索範囲（t0, t1）を上記の事例1（2021年5月11日付近）などに設定してみてください。出力された時刻や距離が上記のデータと一致すれば、コードがマイクロ秒単位の精度で完璧に動作していることが証明できます。

比較検証を行う際は、以下の点にご注意ください。

1. 時間スケールの確認: 上記の事例1・2の時刻は「UTC（協定世界時）」で記録されています。Skyfieldで照合する際は、計算用の時刻はTT（地球時）を維持しつつ、結果の表示を .utc_strftime('%Y-%m-%d %H:%M:%S.%f') にして出力すると比較が容易になります。
2. エフェメリスファイルの選択: このレベルの極めて高い精度（数十センチ、ミリメートル単位）でNASAの出力結果と一致させるには、de421.bsp よりも、より新しい月モデル（重力・マントルの摩擦減衰なども考慮されている）を含む de430.bsp や de440.bsp などを読み込んで使用することが推奨されます。

直接コードを実行してお見せできず恐縮ですが、お手元の環境での高精度な計算テストにぜひこのベンチマークデータをお役立てください。

Skyfieldで最新のJPLエフェメリス（DE440やDE441など）を使用する方法は非常にシンプルで、標準の load() 関数に該当のファイル名を指定するだけです。

これまでのDE421等の代わりに、最新の高精度モデルを読み込ませるための具体的な手順と注意点を解説します。

1. コードの書き方

Skyfieldの load() 関数に目的のエフェメリスファイル名を指定します。Skyfieldは指定されたファイルがローカルにない場合、NASA/JPLのサーバーから自動的にダウンロードして読み込みます。

from skyfield.api import load

# DE440s (DE440の期間短縮・軽量版：1849年〜2150年カバー)
eph = load('de440s.bsp')

# もしくは、フル期間のDE440を使用する場合
# eph = load('de440.bsp')

earth = eph['earth']
moon = eph['moon']

# 以降の計算は今までと同じです

2. 最新エフェメリス（DE440 / DE441）の特徴と選び方

用途に応じて適切なファイルを選択してください。

• DE440 (de440.bsp / de440s.bsp):
  1550年〜2650年をターゲットにした現代の最高精度モデルです。月の液体コアと固体マントルの摩擦による減衰（Damping）が計算に組み込まれており、現代の天文観測データ（Lunar Laser Rangingなど）に最もフィットするように作られています。
  ※通常の占星術計算や現代の現象を追うだけであれば、ファイルサイズが小さく1849年〜2150年をカバーする de440s.bsp（約32MB）が最も使いやすく推奨されます。
• DE441 (de441.bsp):
  紀元前13,200年〜紀元17,191年という超長期間をカバーするモデルです。数千年単位の過去・未来へ遡る際に計算が破綻しないよう、月の減衰モデルがあえて省略されています。歴史的な占星術研究などにはこちらを使用します。

3. Skyfieldのバージョンに関する重要な注意

de441.bsp のように1つの天体に対して複数のデータセグメントを内包している巨大なエフェメリスファイルを正しく読み込むためには、Skyfieldのバージョンが 1.51 以降である必要があります。それより古いバージョンを使用すると、ファイルの半分の日付範囲しか認識されないバグが存在したため、アップデートしておくことをお勧めします。

4. Swiss Ephemerisとの精度比較におけるメリット

前回の回答でも少し触れましたが、Swiss Ephemerisの高精度出力とSkyfieldの出力をミリ秒・ミリメートル単位で正確に一致させたい場合、古い de421.bsp では基礎となる月モデルが古いため、どうしても誤差が残ります。Swiss Ephemeris自体も近年はDE431やDE441といった最新のモデルをベースに最適化されているため、最新の de440.bsp や de441.bsp に切り替えることで、両者の計算結果のズレをさらに極小化することが可能です。