用 Codex 搭建黑洞模拟：天体物理学家的工作流拆解

"""
Schwarzschild 黑洞光线轨迹模拟
计算光子在 Schwarzschild 时空中的轨迹，可视化黑洞"阴影"边界

依赖：numpy, matplotlib
运行：python schwarzschild_raytrace.py
"""

import numpy as np
import matplotlib.pyplot as plt

# ---- 物理参数 ----
M = 1.0          # 黑洞质量（几何单位 G=c=1）
r_s = 2 * M      # Schwarzschild 半径（事件视界）

def effective_potential(r, b):
    """光线在 Schwarzschild 时空中的有效势
    b: impact parameter（碰撞参数）
    """
    return (1 - r_s / r) * (1 + b**2 / r**2)

def trace_ray(b, dr=0.01, r_start=50.0):
    """追踪一条碰撞参数为 b 的光线轨迹
    返回 (r, phi) 数组
    """
    r = r_start
    phi = 0.0
    trajectory = []

    # 光子临界碰撞参数：b_crit = 3*sqrt(3)*M ≈ 5.196
    b_crit = 3 * np.sqrt(3) * M

    while r > r_s * 1.01 and r < r_start * 1.5:
        # dphi/dr = b / (r^2 * sqrt(1 - b^2/r^2 * (1 - r_s/r)))
        denom_sq = 1 - (b**2 / r**2) * (1 - r_s / r)
        if denom_sq <= 0:
            # 光子被捕获，轨迹终止
            break
        dphi_dr = b / (r**2 * np.sqrt(denom_sq))
        phi += dphi_dr * dr
        r -= dr  # 光子向黑洞方向传播
        trajectory.append((r, phi))

    return np.array(trajectory) if trajectory else np.array([]).reshape(0, 2)

# ---- 批量追踪不同碰撞参数的光线 ----
b_values = np.linspace(3.0, 10.0, 40)
trajectories = {}

for b in b_values:
    traj = trace_ray(b)
    if len(traj) > 0:
        trajectories[b] = traj

# ---- 可视化 ----
fig, ax = plt.subplots(figsize=(8, 8))

# 绘制事件视界
theta_circle = np.linspace(0, 2 * np.pi, 100)
ax.plot(r_s * np.cos(theta_circle), r_s * np.sin(theta_circle),
        'k-', linewidth=2, label='事件视界')

# 绘制光线轨迹
for b, traj in trajectories.items():
    r_vals, phi_vals = traj[:, 0], traj[:, 1]
    x = r_vals * np.cos(phi_vals)
    y = r_vals * np.sin(phi_vals)
    color = 'red' if b < 3 * np.sqrt(3) * M else 'blue'
    ax.plot(x, y, color=color, alpha=0.6, linewidth=0.8)

# 光子球标注
r_photon = 3 * M
ax.plot(r_photon * np.cos(theta_circle), r_photon * np.sin(theta_circle),
        'g--', linewidth=1, label='光子球 r=3M')

ax.set_xlim(-15, 15)
ax.set_ylim(-15, 15)
ax.set_aspect('equal')
ax.set_xlabel('x (M)')
ax.set_ylabel('y (M)')
ax.set_title('Schwarzschild 黑洞光线轨迹')
ax.legend()
plt.tight_layout()
plt.savefig('blackhole_shadow.png', dpi=150)
plt.show()

print(f"临界碰撞参数 b_crit = {3*np.sqrt(3)*M:.3f} M")
print("红色轨迹：被黑洞捕获 | 蓝色轨迹：逃逸")

# 用 Codex CLI 生成 CUDA 加速版本的项目骨架
codex "将 schwarzschild_raytrace.py 改写为 CUDA 版本：
- 核心轨迹计算移到 GPU kernel
- 保持 Python 接口，用 PyCUDA 或 cupy 包装
- 添加批量光线并行追踪（100万条光线）
- 输出改为 numpy 数组而非逐条 list"

# 让 Codex 补充测试框架
codex "为 schwarzschild_raytrace.py 添加 pytest 测试：
- 验证 b_crit 附近光线的捕获/逃逸判定
- 验证远距离光线近似直线
- 验证有效势在光子球处极值"