莱斯公式 (Rice Formula) 穿越频率推导全过程可视化

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False

# 设定参数
a_threshold = 10.0; sigma_y = 2.0; sigma_ydot = 3.0; dt = 0.2

# 创建大图 (使用绝对坐标布局)
fig = plt.figure(figsize=(18, 14))

# 1. 最顶端大标题
fig.suptitle('莱斯公式 (Rice Formula) 穿越频率推导全过程可视化', 
             fontsize=24, fontweight='bold', y=0.96)

# --- 第一排 (向上微调) ---
ax1 = fig.add_axes([0.06, 0.68, 0.25, 0.17]) 
ax2 = fig.add_axes([0.38, 0.68, 0.25, 0.17]) 
ax3 = fig.add_axes([0.70, 0.68, 0.25, 0.17]) 

# --- 第二排 (向上微调) ---
ax4 = fig.add_axes([0.06, 0.41, 0.25, 0.17]) 
ax5 = fig.add_axes([0.38, 0.41, 0.25, 0.17]) 

# --- 右侧大图 ---
ax6 = fig.add_axes([0.70, 0.10, 0.25, 0.48]) 

# --- 底部总结文字 (高度缩小，位置下移，彻底避开上面的 x轴标签) ---
ax7 = fig.add_axes([0.06, 0.02, 0.57, 0.22])
ax7.axis('off')

# ==========================================
# 绘制图 1
# ==========================================
y_vals = np.linspace(5, 15, 500)
p_y = (1 / (np.sqrt(2 * np.pi) * sigma_y)) * np.exp(-((y_vals - 10)**2) / (2 * sigma_y**2))
p_a = (1 / (np.sqrt(2 * np.pi) * sigma_y)) * np.exp(-((a_threshold - 10)**2) / (2 * sigma_y**2))

ax1.plot(y_vals, p_y, 'k-', linewidth=2, label='位置概率密度 p(y)')
ax1.axvline(a_threshold, color='r', linestyle='--', label='阈值 a')
for i, v in enumerate([0.5, 1.0, 2.5]):
    w = v * dt
    rect = Rectangle((a_threshold - w, 0), w, p_a, edgecolor='b', facecolor=['lightblue','skyblue','dodgerblue'][i], alpha=0.7)
    ax1.add_patch(rect)
ax1.annotate('窗口厚度 = v·Δt', xy=(9.5, p_a * 0.5), xytext=(7.2, p_a * 1.1),
             arrowprops=dict(facecolor='black', shrink=0.05), fontsize=12, bbox=dict(facecolor='white', alpha=0.8, edgecolor='none'))
ax1.set_title('① 泰勒展开：构造空间窗口', fontsize=15, pad=10)
ax1.set_xlabel('位置 y\n公式: 窗口厚度 = v · Δt', fontsize=12)
ax1.set_ylabel('概率密度 p(y)', fontsize=12)
ax1.set_xlim(7, 13); ax1.set_ylim(0, max(p_y) * 1.4)
ax1.legend(loc='upper right', fontsize=10)

# ==========================================
# 绘制图 2
# ==========================================
ax2.plot(y_vals, p_y, 'k-', linewidth=2)
ax2.axvline(a_threshold, color='r', linestyle='--')
w_demo = 2.0 * dt
rect = Rectangle((a_threshold - w_demo, 0), w_demo, p_a, edgecolor='purple', facecolor='plum', alpha=0.8)
ax2.add_patch(rect)
ax2.annotate('面积 = 概率 P', xy=(a_threshold - w_demo/2, p_a/2), xytext=(7.2, p_a*0.8),
             arrowprops=dict(facecolor='black', shrink=0.05), fontsize=12)
ax2.text(10.2, p_a * 0.4, "频率 = P / Δt\n= p(a)·v·Δt / Δt\n= p(a)·v", fontsize=13, bbox=dict(facecolor='white', alpha=0.9, edgecolor='black'))
ax2.set_title('② 除以时间：转化为频率', fontsize=15, pad=10)
ax2.set_xlabel('位置 y\n公式: 频率 = 概率 / Δt = p(a) · v', fontsize=12)
ax2.set_xlim(7, 13); ax2.set_ylim(0, max(p_y) * 1.4)

# ==========================================
# 绘制图 3
# ==========================================
ydot_vals = np.linspace(0, 10, 500)
p_ydot = np.sqrt(2/np.pi)/sigma_ydot * np.exp(-(ydot_vals**2) / (2 * sigma_ydot**2))
ax3.plot(ydot_vals, p_ydot, 'g-', linewidth=2, label='速度密度 p(v)')
ax3.fill_between(ydot_vals, 0, p_ydot, color='lightgreen', alpha=0.3)
v_pt, dv = 2.5, 0.5
idx1, idx2 = np.searchsorted(ydot_vals, [v_pt - dv/2, v_pt + dv/2])
ax3.fill_between(ydot_vals[idx1:idx2], 0, p_ydot[idx1:idx2], color='green', alpha=0.8)
ax3.text(v_pt, p_ydot[idx1]+0.02, '微小概率\np(v)·dv', ha='center', fontsize=12, color='darkgreen')
ax3.set_title('③ 速度随机：引入 p(v)', fontsize=15, pad=10)
ax3.set_xlabel('穿越速度 v\n公式: 微小概率 dP = p(v) · dv', fontsize=12)
ax3.set_ylabel('概率密度', fontsize=12)
ax3.legend(fontsize=10)

# ==========================================
# 绘制图 4
# ==========================================
integrand = p_ydot * ydot_vals
ax4.plot(ydot_vals, p_ydot, 'g--', linewidth=1, alpha=0.5, label='p(v)')
ax4.plot(ydot_vals, integrand, 'm-', linewidth=2, label='p(v)·v')
ax4.fill_between(ydot_vals, 0, integrand, color='m', alpha=0.2)
integral_val = np.trapezoid(integrand, ydot_vals)
ax4.text(3.5, np.max(integrand)*0.6, f'曲线下面积 = 速度期望\nE[v] ≈ {integral_val:.2f}', fontsize=12, bbox=dict(facecolor='white'))
ax4.set_title('④ 速度加权：计算频率贡献', fontsize=15, pad=10)
ax4.set_xlabel('穿越速度 v\n公式: 频率贡献密度 = p(v) · v', fontsize=12)
ax4.set_ylabel('贡献密度', fontsize=12)
ax4.legend(fontsize=10)

# ==========================================
# 绘制图 5
# ==========================================
v_bins = np.linspace(0, 10, 20)
p_ydot_bins = np.sqrt(2/np.pi)/sigma_ydot * np.exp(-(v_bins**2) / (2 * sigma_ydot**2))
integrand_bins = p_ydot_bins * v_bins
ax5.bar(v_bins, integrand_bins, width=0.4, color='orchid', edgecolor='purple', alpha=0.7)
ax5.plot(ydot_vals, integrand, 'm-', linewidth=2, alpha=0.5)
ax5.text(4.5, np.max(integrand)*0.8, '积分: 无穷多个微小面积之和', fontsize=12, ha='center', bbox=dict(facecolor='white', alpha=0.9))
ax5.set_title('⑤ 积分累加：求期望速度', fontsize=15, pad=10)
ax5.set_xlabel('穿越速度 v\n公式: 期望速度 E[v] = ∫ p(v)·v dv', fontsize=12)
ax5.set_ylabel('面积微元', fontsize=12)

# ==========================================
# 绘制图 6 (右侧大图)
# ==========================================
ax6.set_xticks([]); ax6.set_yticks([])
for spine in ax6.spines.values(): spine.set_visible(False)

ax6.text(0.5, 0.9, "总穿越频率", fontsize=22, ha='center', va='center', fontweight='bold', color='navy')
ax6.text(0.25, 0.65, "位置概率密度\np(a)", fontsize=16, ha='center', va='center', bbox=dict(boxstyle="round,pad=0.5", fc="lightblue", ec="k", lw=2))
ax6.text(0.5, 0.65, "×", fontsize=24, ha='center', va='center')
ax6.text(0.75, 0.65, "速度期望积分\n∫ p(v)·v dv", fontsize=16, ha='center', va='center', bbox=dict(boxstyle="round,pad=0.5", fc="plum", ec="k", lw=2))
result_val = p_a * integral_val
ax6.text(0.5, 0.4, f"最终数值计算结果:\n {p_a:.3f} (1/m) × {integral_val:.3f} (m/s)\n = {result_val:.4f} (次/秒)", fontsize=16, ha='center', va='center', color='darkred', bbox=dict(boxstyle="round,pad=0.6", fc="lightyellow", ec="brown", lw=2))
ax6.text(0.5, 0.15, "核心物理结论:\n穿越频率 = 空间密度 × 运动速度", fontsize=17, ha='center', va='center', color='blue', fontweight='bold')
ax6.set_title('⑥ 终极合成：莱斯公式特例', fontsize=15, pad=10)
ax6.set_xlabel('公式: v_a+ = p(a) · ∫ p(v)·v dv', fontsize=14)

# ==========================================
# 绘制图 7 (底部总结) - 缩小了白框的垂直 padding
# ==========================================
summary_text = (
    "推导过程总结：\n"
    "1. 泰勒展开将极短时间 Δt 转化为了空间窗口厚度 (v·Δt)。\n"
    "2. 概率密度函数 p(a) 乘以该空间厚度，得到了发生穿越的概率。\n"
    "3. 根据定义，概率除以时间 Δt 即可得到频率，此时 Δt 完美约掉，留下速度 v 作为权重。\n"
    "4. 因为真实速度是随机的，需要对速度分布 p(v) 进行积分，求加权平均值。\n"
    "5. 因为在独立假设下，位置密度 p(a) 是常数，可以提到积分号外。\n"
    "量纲守恒验证：(1/米) × (米/秒) = 1/秒 (Hz)"
)
ax7.text(0.02, 0.45, summary_text, fontsize=14, va='center', linespacing=1.8, 
         bbox=dict(facecolor='#f8f9fa', alpha=0.95, edgecolor='#6c757d', boxstyle='round,pad=0.8'))

plt.show()