大概是這樣;such and such

專案簡介

本專案示範如何使用 Python 及開源套件（Astropy、SciPy 等）實作天體力學（Celestial Mechanics），涵蓋

• 現代經典力學：行星運動、開普勒問題
• 廣義相對論修正：近星行星進動、時鐘延遲

範例程式碼基於 Astropy 與 SciPy 撰寫，並搭配 NumPy、Matplotlib 做數值計算與可視化。

---

環境安裝

# 建議使用 virtualenv / conda 建立隔離環境
pip install astropy scipy numpy matplotlib

---

理論背景

1. 經典力學（Newtonian Mechanics）

• 牛頓萬有引力定律：

  $$
  F = G \frac{m_1 m_2}{r^2}
  $$

• 行星的二體問題、開普勒三定律

• 開普勒方程式數值求解

2. 廣義相對論（General Relativity）修正

• 施瓦茲席爾德 (Schwarzschild) 度規下的行星軌道進動

  $$
  \Delta\omega \approx \frac{6\pi GM}{c^2 a(1 - e^2)}
  $$

• 時間延遲（Shapiro Delay）

---

目錄結構範例

celestial_mechanics/
│
├── data/                  # 觀測或模擬資料
├── src/
│   ├── classical.py       # 經典力學相關函式
│   ├── relativity.py      # 相對論修正函式
│   └── utils.py           # 公用工具函式
│
├── notebooks/             # Jupyter 筆記本示例
│   ├── kepler.ipynb
│   └── precession.ipynb
│
├── tests/                 # 單元測試
└── README.md

---

範例一：經典二體問題與開普勒方程式

# src/classical.py
import numpy as np
from scipy.optimize import newton

G = 6.67430e-11  # 萬有引力常數 [m^3 kg^-1 s^-2]

def kepler_equation(E, M, e):
    """開普勒方程式：E - e*sinE = M"""
    return E - e * np.sin(E) - M

def solve_kepler(M, e):
    """
    求解偏近點角 E
    M: 平均近點角 [rad]
    e: 離心率
    """
    return newton(kepler_equation, M, args=(M, e))

def orbital_elements_to_state(a, e, M, mu):
    """
    將軌道要素轉換為位置與速度向量（經典力學）
    a: 半長軸, e: 離心率, M: 平均近點角, mu = G*(m1+m2)
    """
    E = solve_kepler(M, e)
    # 球座標距離與真近點角
    r = a * (1 - e * np.cos(E))
    nu = 2 * np.arctan2(
        np.sqrt(1+e)*np.sin(E/2),
        np.sqrt(1-e)*np.cos(E/2)
    )
    # 平面座標
    x = r * np.cos(nu)
    y = r * np.sin(nu)
    # 速度計算略...
    return np.array([x, y, 0.0]), np.array([-np.sin(E), np.sqrt(1-e**2)*np.cos(E), 0.0]) * np.sqrt(mu*a) / r

---

範例二：廣義相對論進動修正

# src/relativity.py
import numpy as np

def pericenter_precession(a, e, M_central):
    """
    計算近日點進動角度 Δω（rad/週期）
    a: 半長軸 (m)
    e: 離心率
    M_central: 中心天體質量 (kg)
    """
    c = 299792458.0  # 光速 (m/s)
    return 6 * np.pi * G * M_central / (c**2 * a * (1 - e**2))

def shapiro_delay(r1, r2, M_central):
    """
    計算 Shapiro 時延 (秒)
    r1, r2: 訊號傳輸起訖距離向量長度 (m)
    """
    c = 299792458.0
    return 2 * G * M_central / c**3 * np.log((r1 + r2 + np.linalg.norm(r2-r1)) / (r1 + r2 - np.linalg.norm(r2-r1)))

---

Jupyter 筆記本示例

kepler.ipynb

1. 載入 src/classical.py
2. 設定地球–太陽系參數
3. 演示開普勒方程式求解
4. 畫出橢圓軌道位置

precession.ipynb

1. 載入 src/relativity.py
2. 以水星軌道為例，計算近日點進動
3. 與觀測值比較
4. 可視化進動角隨時間之變化

---

可視化範例

import matplotlib.pyplot as plt
from src.classical import orbital_elements_to_state
from src.relativity import pericenter_precession

# 假設參數：以水星為例
a = 5.79e10      # m
e = 0.2056
M_central = 1.9885e30  # 太陽質量 kg
mu = G * M_central

# 計算近日點進動
delta_omega = pericenter_precession(a, e, M_central)
print(f"每軌道進動角: {np.degrees(delta_omega):.6f} 度")

# 繪製一個軌道週期的軌跡
M_vals = np.linspace(0, 2*np.pi, 500)
coords = [orbital_elements_to_state(a, e, M, mu)[0] for M in M_vals]
xs, ys = zip(*[(c[0], c[1]) for c in coords])

plt.figure()
plt.plot(xs, ys)
plt.scatter([0], [0], color='orange', label='中央天體')
plt.axis('equal')
plt.legend()
plt.title("開普勒橢圓軌道示意圖")
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.show()

---

測試與持續整合

• 使用 pytest 撰寫單元測試，驗證 Kepler 解、近日點進動公式

• GitHub Actions 自動化：

  name: CI
  on: [push, pull_request]
  jobs:
    test:
      runs-on: ubuntu-latest
      steps:
      - uses: actions/checkout@v2
      - uses: actions/setup-python@v2
        with: python-version: '3.10'
      - run: pip install -r requirements.txt pytest
      - run: pytest --maxfail=1 --disable-warnings -q
  

---

未來拓展

• 多體問題：N-體模擬、Barnes–Hut 或螺旋法加速
• 數值積分器：Runge–Kutta、Symplectic integrator 比較
• 更深入相對論：利用 EinsteinPy 或自行實作 GR 數值积分
• 資料驅動：整合實際天文觀測資料，進行參數擬合

---

參考資料

• Astropy 官方文件：https://docs.astropy.org/
• SciPy Optimize 模組：https://docs.scipy.org/doc/scipy/reference/optimize.html
• 經典力學與廣義相對論教材
  https://github.com/manjunath5496/Physics-Books
  https://github.com/rgjha/PhysicsBooks