背景
go言語と宇宙工学の勉強のために、ケプラーの方程式を反復計算法（ニュートン ラフソン法）で解いてみます。
ケプラー方程式
まずケプラー方程式は、離心近点離角をE、平均近点離角をM、離心率e とした時に、
で表される方程式です。
地球を焦点とした衛星の楕円軌道を考えたときに以下の図でMとEを表せます。aは半長軸になります。
この方程式は普通の計算では解けず、ニュートン-ラフソン法のような反復方程式などによって解けます。
さて、この方程式を、
と定義すると、
と書くことができ、反復計算ができる形になります。（細かい導出方法はすっ飛ばしてます。参考文献の本に詳細に記述されています。）
golangで実装
では、ここからは実装していきます。
定数の設定
まずは定数を設定します。
// 地球の重力定数
const gravitationalConstant = 398600.4平均近点離角Mを計算する関数
MeanAnomalyという名前で関数を作成します。
func MeanAnomaly(a, deltaT float64)float64{
return math.Sqrt(gravitationalConstant/math.Pow(a, 3)) * deltaT
}aは半長軸(km)を表しており、deltaTは 元期 から求めたい瞬間の時刻の経過時間(秒)になります。
この関数によって平均点離角（radian）をreturnしています。
ケプラー方程式をクロージャとして定義
ケプラー方程式をクロージャとして定義し、後ほど反復法の関数にて利用していきます。
// ケプラー方程式
func KeplerEquation(e, M float64) func(Ebefore float64) float64{
return func(Ebefore float64)float64{
FE := Ebefore - e*math.Sin(Ebefore) - M
Eafter := Ebefore - FE/(1-e*math.Cos(Ebefore))
return Eafter
}
}それぞれの引数は以下を表しています。
Ebefore：n回目の離心近点離角
M：平均近点離角
e：離心率
Eafter：n+1回目の離心近点離角
これらの引数を入力して宣言することで、return文に書いた関数を都度実行できるようになります。
反復計算法（ニュートン ラフソン法）
では上記で実装したケプラー方程式のクロージャをニュートン-ラフソン法で計算する関数を実装します。
許容誤差(ae)を引数として指定し、計算誤差が許容誤差よりも小さくなったら、forループを抜け、計算結果をreturnするようにしています。
// ケプラー方程式をニュートン-ラフソン法で解く
func NewtonRaphson(e, before, ae float64) float64 {
// ae: allowable error 許容誤差
equation := KeplerEquation(e, before)
var after float64
err := 100.0 // 誤差の初期化（0だとforが回らないため）
count := 0
for ; err > ae; {
after = equation(before)
err = math.Abs(after - before)
fmt.Println("誤差: ", err)
before = after
count = count + 1
fmt.Println("繰り返し：", count, "回目")
fmt.Println("----")
}
return after
}main関数
ここまでで、必要な関数の定義は終わりです。
あとはmain関数で実際に使ってみて、正しい計算結果が得られるか確認してみます。
半長軸aを2000km、経過時間deltaTを3600秒、離心率eを0.75、許容誤差を0.0001として計算してみます。
func main() {
fmt.Println("ケプラー方程式を使って離心近点離角Eを求めます：")
M := MeanAnomaly(2000, 3600)
Mdigree := int(M * 180 / math.Pi) % 360
E := NewtonRaphson(0.75, M, 0.0001)
Edigree := int(E * 180 / math.Pi) % 360
Mrad := math.Round(M*100)/100
Erad := math.Round(E*100)/100
fmt.Printf("\n平均近点離角Mは %vrad、%v° です。", Mrad, Mdigree)
fmt.Printf("\n離心近点離角Eは %vrad、%v° です。", Erad, Edigree)
}出力結果は
ケプラー方程式を使って離心近点離角Eを求めます：
誤差: 0.7393437287609679
繰り返し： 1 回目
----
誤差: 0.1667982273838291
繰り返し： 2 回目
----
誤差: 0.016917352484746573
繰り返し： 3 回目
----
誤差: 0.00016185780524580196
繰り返し： 4 回目
----
誤差: 1.4669140568912553e-08
繰り返し： 5 回目
----平均近点離角Mは 25.41rad、15° です。
離心近点離角Eは 25.97rad、47° です。
となります。
下記のサイトによる計算結果と比較してみると、少しの計算誤差はあるものの、ほぼ正しく計算できていることが分かります。
コード全体
以下にコード全体を残しておきます。
package main
import (
"fmt"
"math"
)
// 地球の重力定数
const gravitationalConstant = 398600.4
// 平均近点離角
func MeanAnomaly(a, deltaT float64)float64{
return math.Sqrt(gravitationalConstant/math.Pow(a, 3)) * deltaT
}
// ケプラー方程式
func KeplerEquation(e, M float64) func(Ebefore float64) float64{
// Ebefore: n回目の離心近点離角
// M: 平均近点離角
// e: 離心率
// Eafter: n+1回目の離心近点離角
return func(Ebefore float64)float64{
FE := Ebefore - e*math.Sin(Ebefore) - M
Eafter := Ebefore - FE/(1-e*math.Cos(Ebefore))
return Eafter
}
}
// ケプラー方程式をニュートン-ラフソン法で解く
func NewtonRaphson(e, before, a float64) float64 {
// M: 平均近点離角
// e: 離心率
// a: allowable error 許容誤差
equation := KeplerEquation(e, before)
var after float64
err := 100.0 // 許容誤差の初期化（0だとforが回らないため）
count := 0
for ; err > a; {
after = equation(before)
err = math.Abs(after - before)
fmt.Println("誤差: ", err)
before = after
count = count + 1
fmt.Println("繰り返し：", count, "回目")
fmt.Println("----")
}
return after
}
func main() {
fmt.Println("ケプラー方程式を使って離心近点離角Eを求めます：")
M := MeanAnomaly(2000, 3600)
Mdigree := int(M * 180 / math.Pi) % 360
E := NewtonRaphson(0.75, M, 0.0001)
Edigree := int(E * 180 / math.Pi) % 360
Mrad := math.Round(M*100)/100
Erad := math.Round(E*100)/100
fmt.Printf("\n平均近点離角Mは %vrad、%v° です。", Mrad, Mdigree)
fmt.Printf("\n離心近点離角Eは %vrad、%v° です。", Erad, Edigree)
}参考文献
以下の本で理論を参考にしました。
https://www.amazon.co.jp/exec/obidos/ASIN/4768704395/detethcotaf-22/
