大概是這樣;such and such

## 專案簡介

本專案示範如何使用 Python 及開源套件（Astropy、SciPy 等）實作天體力學（Celestial Mechanics），涵蓋
- **現代經典力學**：行星運動、開普勒問題
- **廣義相對論修正**：近星行星進動、時鐘延遲

範例程式碼基於 [Astropy](https://github.com/astropy/astropy) 與 SciPy 撰寫，並搭配 NumPy、Matplotlib 做數值計算與可視化。

---

## 環境安裝

```bash
# 建議使用 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
```

---

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

```python
# 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
```

---

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

```python
# 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. 可視化進動角隨時間之變化

---

## 可視化範例

```python
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 自動化：

```yaml
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/](https://docs.astropy.org)
* SciPy Optimize 模組：[https://docs.scipy.org/doc/scipy/reference/optimize.html](https://docs.scipy.org/doc/scipy/reference/optimize.html)
* 經典力學與廣義相對論教材
https://github.com/manjunath5496/Physics-Books
https://github.com/rgjha/PhysicsBooks