Excelで解くケプラーの方程式

「J2000.0から視位置への変換（太陽編）」で年周光行差をちゃんと計算しなくちゃ、という話になったのですが、年周光行差をちゃんと計算しようと思うと「近日点黄経」が必要になります。
　　太陽の真近点角 = 太陽黄経 - 近日点黄経
という関係があります。
「太陽黄経」は「太陽の黄経・黄緯を求め月相を計算する」に書いたように赤経・赤緯からすぐに（とは言ってもそのためには真黄道傾斜角を求める必要がありますから、“すぐに”ではあっても“簡単に”ではありませんが）求めることができます。
つまり真近点角を求めれば近日点黄経はすぐに求められます。真近点角は離心近点角から求まります。そして離心近点角は平均近点角から求めることができます。そして平均近点角は章動の計算をするときに求めています。
なんだか遠回りしているような気もしますが、これで近日点黄経は求められそうです。

ところで平均近点角から離心近点角を求めるところで「ケプラーの方程式」
　　u - e sin u = l
を解く必要があります（u: 離心近点角、l: 平均近点角）
そこで今回はExcelで「ケプラーの方程式」を解く方法です。

実例（Excelのシート）は数日中に書く予定の「J2000.0から視位置への変換」の続編にありますので今回は省略します。

なお離心近点角から真近点角の求め方はExcel風に書くと次のようになります。
　　真近点角
　　　　=ATAN2((COS(離心近点角)-離心率)/(1-離心率*COS(離心近点角)),
　　　　　　SQRT(1-離心率^2)*SIN(離心近点角)/(1-離心率*COS(離心近点角)))

atan2()を使っていますので象限がどうのこうのというのは考える必要はありません。

--------

ケプラーの方程式は解析的には解けないので数値的に解を求める必要があります。

離心近点角uの初期値を適当に決めて上式の残差

　　R= l - u + e sin u

を求め、残差Rが0に（近く）なるように計算を繰り返すことにします。

残差 R を u で微分すると

　　dR/du = -1 + e cos u

となりますから、まず

　　R1 = l - u1 + e sin u1

を求め

　　Δu1 = - R1 / ( dR/du ) = - R1 / ( -1 + e cos u1 )
　　．
　　．
　　Δui = - Ri / ( dR/du ) = - Ri / ( -1 + e cos ui )

から

　　u2 = u1 + Δu1
　　．
　　．
　　ui = ui-1 + Δui-1

として新たな R を求めるということを繰り返しRが十分小さくなるu2,..,uiを求めればいいはずです。

u1をどう決めるか問題ですが

　　平均近点角≒離心近点角

と考え

　　u1 = l + e sin l

とします。

Excelでこれを計算してみます。

各セルに次のように計算式をセットします。
離心率eは0.01672とします。これは地球のじっさいの離心率です。
平均近点角l=1（57.29...度）の場合を計算します。

じっさいに計算するとあっけなく収束してしまい離心近点角 u=1.01420 が求まります。



-----

ところでExcelにはソルバーという機能があります。これを使うともっと安直です。

「ソルバー」はアドインなので最初にインストールする必要があるのですが、一度インストールしてしまえばExcelの一つの機能としてふつうに使えます。

次のようなシートを作ります。



データ=>分析=>ソルバーからソルバーを起動します。

uつまりA5を変化させてRつまりB5セルを0にすることが目標ですから
　　目的セル： ＄B$5
　　目標値： 値： 0
　　変化させるセル： ＄A$5
と設定します。



そして実行ボタンをクリックすると次の結果が一瞬で得られます。



　　u=1.01420

となりもちろん上の繰り返し計算を行ったときと同じ結果が得られます。

関連
　　「続・J2000.0から視位置への変換（太陽編）」
　　実際にケプラーの方程式を利用しています。

参考
　　『長沢工「天体の位置計算」地人書館』