惑星と小天体の軌道を求める方法を説明します。

観測時刻を定義する

観測時刻 T_eph[JED] はユリウス日とする。(JED : Julian Ephemeris Date)

年 月 日 時 分 秒

世界時からの時差 ＋ － 時間

T_eph = [JED]

軌道要素を定義する

惑星 小天体

惑星小天体
軌道の形状 軌道長半径 a0[au]

軌道長半径 a_dot[au/century]

離心率 e0

離心率 e_dot[1/century]

軌道長半径 a[au]

離心率 e
軌道の存在する平面 軌道傾斜角 I0[deg]

軌道傾斜角 I_dot[deg/century]

昇交点経度 Omega0[deg]

昇交点経度 Omega_dot[deg/century]

軌道傾斜角(I)[deg]

昇交点経度(Omega)[deg]
軌道の向き 近日点経度(varpi0)[deg]

近日点経度(varpi_dot)[deg/century]

 近日点引数(omega)[deg]
軌道上の位置

平均経度 L0[deg]

平均経度 L_dot[deg/century]

西暦が1800年未満または2050年以降のとき

木星から冥王星までの平均近点角を求める係数
b * T^2
b =

木星から海王星までの平均近点角を求める係数
c * cos(f * T), s * sin(f * T)
f =
c =
s =

元期(Epoch)[JED]

平均近点角(M0)[deg]

惑星の軌道要素は以下のページを参照しました。
NASA : Keplerian Elements for Approximate Positions of the Major Planets

赤道座標を計算する

計算ボタンを押すと赤道座標が表示されます。計算手順の詳細は以下に記載しました。

赤経 h m s ( °)
赤緯 ° m s

赤経(J2000.0) h m s ( °)
赤緯(J2000.0) ° m s

計算結果を比較するために下記ページの参照をおすすめします。
CASIO : 惑星の黄道系、赤道系の位置計算

このページで計算した結果と比較して惑星は1分角程度の誤差が発生することは避けられません。
また、小天体は軌道要素が異なると位置が大幅にずれる可能性があります。
精度が高い位置を求めることは困難であることがこの結果からもわかると思います。

ガウス引力定数を求める

1天文単位 AU = 149597870691[m] とする。
太陽中心重力定数 GM_Sun[m^3*s^-2] は、
万有引力定数G[m^3*kg^-1*s^-2]と太陽の質量M_Sun[kg]の積である。
したがって、GM_Sun = 1.32712440018e+20
ガウス引力定数 K0 は、太陽を回る天体の軌道長半径[au]の1.5乗と平均運動n[rad/day]の積となる。
したがって、K0 = sqrt(GM_Sun / (AU * AU * AU)) * 86400 = 0.01720209895

近日点引数を求める (惑星のみ)

時刻 T_eph[JED] における近日点引数 omega[rad] は、
近日点引数 omega[rad] = 近日点経度 varpi[rad] - 昇交点経度 Omega[rad]

2000年1月1日正午を起点としたユリウス世紀
T[century] = (T_eph - 2451545.0) / 36525

近日点経度 varpi[rad]
= 近日点経度 varpi0[rad] + 近日点経度 varpi_dot[rad/century] * T[century]

昇交点経度 Omega[rad]
= 昇交点経度 Omega0[rad] + 昇交点経度 Omega_dot[rad/century] * T[century]

軌道面座標を黄道面座標に変換するための行列を求める

軌道面座標を黄道面座標に変換するための行列 Matrix は、
Matrix = (M11 M12), (M21 M22), (M31 M32)

cos_I = cos(I)
sin_I = sin(I)
cos_Omega = cos(Omega)
sin_Omega = sin(Omega)
cos_omega = cos(omega)
sin_omega = sin(omega)

M11 = cos_omega * cos_Omega - sin_omega * sin_Omega * cos_I
M12 = -sin_omega * cos_Omega - cos_omega * sin_Omega * cos_I
M21 = cos_omega * sin_Omega + sin_omega * cos_Omega * cos_I
M22 = -sin_omega * sin_Omega + cos_omega * cos_Omega * cos_I
M31 = sin_omega * sin_I
M32 = cos_omega * sin_I

平均運動を求める

小天体の平均運動 n[rad/day] は、

楕円軌道のとき (離心率e < 1)
n = K0 / (a * sqrt(a))

楕円軌道のときは軌道短半径b[au]も求める
b = a * sqrt(1 - e * e)

双曲線軌道のとき (離心率e > 1)
n = K0 / (-a * sqrt(-a))

平均近点角を求める

時刻 T_eph[JED] における平均近点角 M[rad] は、

惑星のとき

2000年1月1日正午を起点としたユリウス世紀
T[century] = (T_eph - 2451545.0) / 36525

平均経度 L[rad] = 平均経度 L0[rad] + 平均経度 L_dot[rad/century] * T[century]

M[rad] = 平均経度 L[rad] - 近日点経度 varpi[rad]

小天体のとき

M[rad] = M0 + n * (T_eph - Epoch)

離心近点角を求める

平均近点角が M[rad] のとき、離心近点角 E[rad] は、ケプラーの運動方程式を解いて、

(1) E0 = M
(2) delta_E = (M - E0 + e * sin(E0)) / (1 - e * cos(E0))
(3) E = E0 + delta_E
(4) E0 = E

delta_E->0(delta_Eの絶対値が0.000001未満)になるまで(2)から(4)を繰り返す。

日心軌道面座標を求める

近日点方向をx_prime軸とする日心軌道面座標 r_prime[au,au,0] は、

x_prime = a * (cos(E) - e)
y_prime = b * sin(E)
r_prime = (x_prime y_prime 0)

日心黄道座標を求める

春分点方向をx軸とする黄道面上の日心黄道座標 r_ecl[au,au,au] は、
r_ecl = Matrix * r_prime

x = M11 * x_prime + M12 * y_prime
y = M21 * x_prime + M22 * y_prime
z = M31 * x_prime + M32 * y_prime
r_ecl = (x y z)

黄道傾斜角を求める

時刻 T_eph[JED] における黄道傾斜角 epsilon[rad] は、

2000年1月1日正午を起点としたユリウス世紀
T[century] = (T_eph - 2451545.0) / 36525

epsilon[deg] = (84381.406 - 46.836769 * T - 0.00059 * T^2 + 0.001813 * T^3) / 3600
epsilon[rad] = epsilon[deg] * PI / 180
ただし、歳差運動を考慮しないJ2000.0のときはT=0となるので、 84381.406 / 3600 [deg] とする。

赤道座標を求める

赤道座標を求める天体の日心黄道座標を r_ecl = (x y z)、
地球の日心黄道座標を r_ecl_earth = (x0 y0 z0) とすると、
その天体の赤道座標である赤経 RA[rad] と 赤緯 Dec[rad] は、

x1 = x - x0
y1 = (y - y0) * cos(epsilon) - (z - z0) * sin(epsilon)
z1 = (y - y0) * sin(epsilon) + (z - z0) * cos(epsilon)
RA = atan(y1 / x1)
Dec = atan(z1 / sqrt(x1^2 + y1^2))

x1 > 0 かつ y1 < 0 のとき
RA = RA + PI * 2

x1 < 0 のとき
RA = RA + PI

以上でJ2000.0の赤道座標を求めることができます。
歳差運動を考慮する場合は、この結果を補正する必要があります。