« ICRS(J2000.0)から視位置への変換(太陽編) | トップページ | 天文計算のための天文用語集 - 1 »

2014年5月 8日 (木)

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...度)の場合を計算します。
001
じっさいに計算するとあっけなく収束してしまい離心近点角 u=1.01420 が求まります。
002


-----

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

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

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


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

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


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


  u=1.01420

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

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

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

« ICRS(J2000.0)から視位置への変換(太陽編) | トップページ | 天文計算のための天文用語集 - 1 »

天文計算」カテゴリの記事

コメント

I gotta favorite this website it seems very useful very helpful

コメントを書く

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.f.cocolog-nifty.com/t/trackback/1994983/56098823

この記事へのトラックバック一覧です: Excelで解くケプラーの方程式:

« ICRS(J2000.0)から視位置への変換(太陽編) | トップページ | 天文計算のための天文用語集 - 1 »

フォト

サイト内検索

  • 記事を探されるんでしたらこれがいちばん早くて確実です。私も使ってます (^^;; 検索窓が表示されるのにちょっと時間がかかるのはどうにかしてほしいです。

新着記事

リンク元別アクセス数

  • (アクセス元≒リンク元、原則PCのみ・ドメイン別、サイト内等除く)

人気記事ランキング

  • (原則PCのみ、直近2週間)
無料ブログはココログ