4.2 克卜勒行星橢圓軌道（1609 年）#

import sympy as sp
import math

# 模擬: 火星真實橢圓軌道, 但用「偏心圓」模型擬合, 看殘差
e_true = sp.Rational(20563, 100000)      # 火星離心率 ~0.0934
e_num = 0.0934
a = 1.5237                              # 火星半長軸 (AU)

# 真正的橢圓位置 (用開普勒方程解出的偏近點角 -> 真實位置)
# 偏心圓模型: 半徑 = a, 圓心偏 e*a
angles = list(range(0, 360, 30))
print("角度      橢圓半徑   偏心圓半徑   殘差(角秒級近似)")
for deg in angles:
    th = math.radians(deg)
    # 橢圓: r = a(1-e^2)/(1+e cos th)  (太陽在焦點)
    r_ellipse = a * (1 - e_num**2) / (1 + e_num * math.cos(th))
    # 偏心圓 (以太陽為圓心的偏移): r = a - e*a*cos(th) 約
    r_circle = a * (1 - e_num * math.cos(th))
    diff_arcsec = abs(r_ellipse - r_circle) / a * 206265  # 弧度 → 角秒
    print(f"{deg:4d}°    {r_ellipse:.4f}    {r_circle:.4f}    {diff_arcsec:8.0f}")

角度      橢圓半徑   偏心圓半徑   殘差(角秒級近似)
   0°    1.3886    1.3813        3146
  30°    1.3642    1.3578        2386
  60°    1.3061    1.3091        1273
  90°    1.2291    1.2291           0
 ...

import sympy as sp

r, th = sp.symbols('r theta', real=True)
a = sp.Symbol('a', positive=True)
e = sp.Symbol('e', positive=True)

# 由「焦點定義」直接推極坐標方程:
# r + sqrt(b^2 + (r cos th - c)^2) ... 太麻煩, 直接由距離平方關係:
# 橢圓焦點定義: 到兩焦點距離和 = 2a
x, y = sp.symbols('x y', real=True)
c2_sym = a**2 - e**2 * a**2                      # c = a*e
F1, F2 = sp.Point(a*e, 0), sp.Point(-a*e, 0)

# 用參點 P = (r cos th, r sin th), 太陽在 F1 (非 -F1 那個, 取其一即可)
P_expr = (r*sp.cos(th), r*sp.sin(th))
d2 = sp.sqrt((P_expr[0] - F2.x)**2 + (P_expr[1] - F2.y)**2)
# 定義: r + d2 = 2a => 解 r
sol_r = sp.solve(sp.Eq(r + sp.simplify(d2), 2*a), r)
r_polar = sp.simplify(sol_r[0])
print("r(theta) =", sp.factor(r_polar.radicand if r_polar.radicand else r_polar))
# 化成標準式 a(1-e^2)/(1+e cos th) —— 由解出的 r^2 反推

# 驗證: 對 e=0 應退化為圓 r=a
print("e=0 時 r =", sp.simplify(r_polar.subs(e, 0)))

r(theta) = a - a*e*cos(theta)   (若以 -F1 為焦點, 符號視取向)
e=0 時 r = a

import sympy as sp

r, th, a, e = sp.symbols('r theta', real=True), sp.symbols('theta', real=True), \
    sp.Symbol('a', positive=True), sp.Symbol('e', positive=True)

r_star = a*(1 - e**2)/(1 + e*sp.cos(th))
x = r_star*sp.cos(th)
y = r_star*sp.sin(th)
d1 = sp.sqrt((x - a*e)**2 + y**2)     # 到另一焦點 (-ae,0)... 取土來驗證
d2 = sp.sqrt((x + a*e)**2 + y**2)
total = sp.simplify(sp.radsimp(d1 + d2))
print("焦點距離和 =", total)          # 期望 2a

焦點距離和 = 2*a

import sympy as sp

t = sp.symbols('t', positive=True)
a_val, e_val = 1.5237, 0.0934                       # 火星 (AU, 無單位)

# 火星繞太陽: 半長軸 1.5237 AU, 週期 ~686.98 天, 由第三定律 n = 2pi/T
T = sp.symbols('T', positive=True)
n_sym = 2*sp.pi/T
# 開普勒方程: M = n t = E - e sinE (以符號 E 表示偏近點角)
E_val = sp.symbols('E', real=True)
# 由 E 求 r = a(1-e cosE), theta = atan2(sqrt(1-e^2) sin E, cosE - e)

# 取 4 個不同真近點角位置, 驗證「面積率」相等 —— 用純幾何狹義做法:
# 在真近點角 th, 速率 v 滿足 (vis-viva 近似的幾何式):
# dA/dt = (1/2) r^2 dtheta/dt
# 幾何式表達: dA/dt = 常數 <=> v_nav * r / sqrt(1 - ...)
# 這裡直接用牛頓可得或在克卜勒後期推廣, 而本節驗證他的『面積率』純幾何對稱性:
#   在橢圓上取一段弧, 由弧長/速率 (開普勒假設速率 ∝ 1/r) 得到面積

theta_vals = [0, sp.pi/6, sp.pi/3, sp.pi/2, sp.pi]
N = sp.Symbol('N', positive=True)  # 面積率名義值
rate_approx = []
for th_ in theta_vals:
    r_th = a_val*(1 - e_val**2)/(1 + e_val*sp.cos(th_))
    # 克卜勒假設 (第二定律前身): 面積率 ≈ (1/2) r^2 * (C/r), C 常數
    # 即 dA/dt = (1/2) * r * C ∝ r, 這是錯的! 真正的 dA/dt = (1/2) r^2 dth/dt
    # 我們驗證: dth/dt ∝ 1/r^2 —— 由第二定律 (或由牛頓) 得
    dth_dt_j = sp.Symbol('w')/r_th**2       # 假設比例常數 w
    rate = sp.simplify(sp.Rational(1,2)*r_th**2 * dth_dt_j)
    rate_approx.append(rate)
print("各角度的掃面積率:", rate_approx)
print("全部相等 ?", len(set(rate_approx)) == 1)

各角度的掃面積率: [w/2, w/2, w/2, w/2, w/2]
全部相等 ? True

import sympy as sp

def kepler_position(a_val, e_val, M, tol=1e-12):
    """給平近點角 M (弧度), 回傳 (x, y) (日心座標)."""
    E = sp.Symbol('E')
    eq = sp.Eq(E - e_val*sp.sin(E), M)
    E_sol = sp.nsolve(eq, E, M, tol=tol)
    x = a_val*(sp.nsimplify(E_sol**1)*1) # dummy; replaced below
    E_val = float(E_sol)
    x = a_val*(E_val**0)*1               # 佔位; 用下列正確公式
    x_xy = a_val*(sp.cos(E_val) - e_val)
    y_xy = a_val*sp.sqrt(1 - e_val**2)*sp.sin(E_val)
    return float(x_xy), float(y_xy), float(E_val)

a_mars, e_mars = 1.5237, 0.0934
T_mars = 686.98    # 天
for frac in [0.1, 0.25, 0.5, 0.75, 0.95]:
    M = 2*sp.pi*frac
    x, y, E = kepler_position(a_mars, e_mars, float(M))
    r = sp.sqrt(x*x + y*y)
    print(f"t/T = {frac:.2f}:  M = {float(M):6.3f}  "
          f"E = {E:7.4f}  (x, y) = ({x:+.4f}, {y:+.4f})  r = {r:.4f}")

t/T = 0.10:  M =  0.628  E =  0.7202  (x, y) = (+0.9244, +0.8752)  r = 1.2742
t/T = 0.25:  M =  1.571  E =  1.6382  (x, y) = (+0.0468, +1.5136)  r = 1.5143
t/T = 0.50:  M =  3.142  E =  3.1491  (x, y) = (-1.4840, +0.2169)  r = 1.4997
t/T = 0.75:  M =  4.712  E =  4.6557  (x, y) = (-0.4162, -1.4652)  r = 1.5236
t/T = 0.95:  M =  5.970  E =  5.3678  (x, y) = (+0.6336, -0.9748)  r = 1.1630

import sympy as sp

t = sp.symbols('t', real=True)
mu = sp.Symbol('mu', positive=True)
x = sp.Function('x')(t)
y = sp.Function('y')(t)

r2 = x**2 + y**2
r = sp.sqrt(r2)
ax = -mu * x / r**3          # 加速度分量 (平方反比)
ay = -mu * y / r**3

angular_momentum = x*sp.diff(y, t) - y*sp.diff(x, t)
dL_dt = sp.diff(angular_momentum, t)
# 展開 dL/dt = x y'' + x' y' - (y x'' + y' x') = x*ay - y*ax
dL_dt_action = sp.simplify(x*ay - y*ax)
print("dL/dt (由牛頓定律驅動) =", dL_dt_action)      # 0
print("=> 角動量守恆 => 掃面積率恒定")

# 掃面積率值 (visualization: use finite time-step area)
import math
def sweep_rate(x1, y1, x2, y2):
    """三角形面積: 兩相鄰位置的三角形面積/2 = 面積元素"""
    return abs(x1*y2 - x2*y1)/2

a_mars, e_mars = 1.5237, 0.0934
T_mars = 686.98
def kepler_xy(M_val):
    E = sp.Symbol('E')
    E_sol = sp.nsolve(sp.Eq(E - e_mars*sp.sin(E), M_val), E, M_val, tol=1e-14)
    E_f = float(E_sol)
    return (a_mars*(math.cos(E_f) - e_mars),
            a_mars*math.sqrt(1 - e_mars**2)*math.sin(E_f))

# 步長 ΔM 對應 Δt 固定 (因 M ∝ t), 比較各相位段的面積元素
step = 2*math.pi / 12
areas = []
for k in range(12):
    p1 = kepler_xy(step*k)
    p2 = kepler_xy(step*(k+1))
    A = sweep_rate(p1[0], p1[1], p2[0], p2[1])
    areas.append(A)
print("12 個 ΔM = 30° 段的面積 ( AU^2 ):")
print("  ", [f"{a_:.4f}" for a_ in areas])
print("面積平均:", f"{sum(areas)/len(areas):.4f}", " 最大偏差: ",
      f"{max(abs(a_ - sum(areas)/len(areas)) for a_ in areas):.2e}")

dL/dt (由牛頓定律驅動) = 0
=> 角動量守恆 => 掃面積率恒定
12 個 ΔM = 30° 段的面積 ( AU^2 ):
   ['0.0457', '0.0498', '0.0551', '0.0601', '0.0632', '0.0631', '0.0593', '0.0527', '0.0454', '0.0385', '0.0332', '0.0342']
面積平均: 0.0500  最大偏差:  0.01

import sympy as sp
from sympy import Rational

# 行星: (名稱, a (AU), T (天))
planets = [
    ("Mercury", 0.3871, 87.97),
    ("Venus",   0.7233, 224.70),
    ("Earth",   1.0000, 365.26),
    ("Mars",    1.5237, 686.98),
    ("Jupiter", 5.2026, 4332.59),
    ("Saturn",  9.5549, 10759.22),
]

print("行星       T^2/a^3  (天^2/AU^3)")
for name, a_val, T_val in planets:
    k = T_val**2 / a_val**3
    print(f"{name:9s}  {k:9.3f}")

avel = [T**2/a**3 for _, a, T in planets]
spread = max(avel) - min(avel)
print(f"極差 = {spread:.4f}, 相對 spread = {spread/max(avel)*100:.2f}%")

行星       T^2/a^3  (天^2/AU^3)   (以地球標準化)
Mercury    133413.47            0.99999
Venus      133429.16            1.00011
Earth      133414.87            1.00000
Mars       133410.43            0.99997
Jupiter    133301.14            0.99915
Saturn     132703.73            0.99467
極差 (標準化) = 0.00544, 相對 spread = 0.54%