2つの天体が互いに万有引力で運動することを2体問題といいます。太陽系の惑星の運動など，いわゆるケプラーの法則に従う運動で，円ではなく楕円軌道を描くことが大きな発見だったわけです。この運動を再現する（天体の位置を計算で求める）方法は，現在ではコンピューターの数値計算で（中心力を万有引力として計算することが）可能[1]ですが，暦計算的には今でもケプラーの時代と同じ古典的な方法が用いられています。公転する天体の楕円軌道に外接する円への垂線と天体が交わる点$Q’$と，円の中心からの角度（離心近点角）$u$と，楕円の焦点と天体の真近点角$f$の関係式から位置を求めるものです。これは星座早見プログラムでも使っている，公転運動の平均近点離角$nt$から離心近点角$u$を求めるケプラー方程式によって解を求める方法です。

天文学事典https://astro-dic.jp/による

楕円の性質として重要なのが，離心率$e$です。楕円には中心$O$と，２つの焦点$F$（$O$に対象）があります。離心率$e$は，この図で$\frac{OF}{OP}$の比で，０から１までの値をとり，楕円のつぶれ具合（0は真円で，1はつぶれた状態）を示します。楕円の大きさと形は，長半径$OP$の長さ（$a$で表わす）と$e$で決まります。高校の数学では，ここまでですが，平面座標を使うより半径と角度で表わす極座標で楕円を表わすと，$r=\frac{l}{1 + e \cos \theta}$で表わされます。ここで$\theta$は$f$で，$l$は，図の$F$から$y$軸と楕円の交点までの長さで，半直弦と呼ばれます。ちょっと説明が大変になるのですが，これらの図形と角度との関係を調べていくと，さらに，運動する惑星と中心天体を結ぶ動径$r$は，$r=a(1-e\cos u)$で表わされ，$u$と$f$の関係式も求まります[2]。そして，離心近点角と呼ばれる$u$ですが，惑星の時間にともなう平均運動での位置，平均近点離角$nt$に対して，$nt=u-e\sin u$という関係式が成り立ちます。これをケプラー方程式と言い，この解としての$u$の値が分かれば，$r$，と$f$が求められ惑星（天体）の位置が計算できます。この方程式を解く計算は，はじめに$nt$の値を$u$の近似値として繰り込み計算を行うことで求まります。ケプラーの時代から行われていた方法で，今でこそパソコンであっという間にできますが，ケプラーは，このような方法を駆使して過去の惑星の位置を手計算で求め，聖書にあるベツレヘムの星を紀元前7年の木星と土星の会合現象だと指摘したことに驚かされます。それから，ニュートンの「プリンキピア」が微分や積分を使わず，幾何学的に運動を説明していると言うことが，とても複雑とはいえ，説得力があることだなと思いました。

参考にしました。[1]https://cattech-lab.com/science-tools/simulation-lecture-2-5/
　　　　　　　　[2]http://fnorio.com/0158Kepler_equation/Kepler_equation.html
　　　　　　　　[3]長沢工(1985)：増補版　天体の位置計算，地人書館

離心率0.8のとき
離心率0.4のとき
'''ケプラーの楕円運動をアニメーションで示すプログラム
２０２１年3月5日作成'''

import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def Kepler_EQ(ec,ma):       # ケプラー方程式の解を求める関数
    ex = ma
    ss = ma
    while abs(ss) > 1.0e-8:
        ss = ex - ec * math.sin(ex) - ma
        de = ss / (1 - ec * math.cos(ex))
        ex = ex - de
    return ex

# 軌道要素
MU = math.radians(1.0)   # 日日平均運動 地球の場合は0.9856076°
AX = 1.0    # 軌道長半径  地球は１天文単位
EC = 0.8    # 離心率

# AxesとFigureを設定
fig, ax = plt.subplots(figsize = (9.0, 6.0))
ax.set_xlim(-2.0, 1.0)
ax.set_ylim(-1.0, 1.0)
ax.set_xlabel("x", fontsize = 15)
ax.set_ylabel("y", fontsize = 15)
ax.axhline(y=0)
ax.axvline(x=0)
ax.set_title("AX=1.0  EC=0.8")
p = []

def dot_Ellips(frame):
    MA = MU * frame    # MA:平均近点角（平均運動の角度　=　フレーム値で調整）
    EX = Kepler_EQ(EC, MA)   # 離心近点角EXをケプラー方程式で求める
    T = math.sqrt((1 + EC) / (1 - EC)) * math.tan(EX / 2)
    f = 2 * math.atan(T)   # 真近点角ｆ
    r = (1 - EC * EC) / (1 + EC * math.cos(f))    # ｒ:動径
    X = r * math.cos(f)  # 近点をｘ座標とした位置
    Y = r * math.sin(f)
    if len(p) > 0:    # 前の表示を消す処理
        item = p.pop(0)
        ax.lines.remove(item)   #
    p.append(ax.plot(X, Y, "bo", c="blue")[0])

ani = animation.FuncAnimation(
      fig,  # Figureオブジェクト
      dot_Ellips,  # グラフ描画関数
      frames = np.arange(0,360,10),  # フレームを設定 360°を10°毎に表示
      interval = 40,  # 更新間隔(ms)
      repeat = True,  # 描画を繰り返す
      )
ani.save("c04.gif", writer="imagemagick")
plt.show()