主题077：量子效应在辐射换热中的作用

摘要

量子效应在纳米尺度辐射换热中扮演着至关重要的角色。本主题深入探讨量子力学原理如何影响热辐射过程，包括零点能、量子涨落、光子统计和量子纠缠等现象对辐射换热的影响。通过建立量子辐射传递模型，分析量子效应在纳米间隙传热、量子热机和量子热管理中的应用，并通过Python仿真展示量子效应的独特现象。

关键词

量子辐射、零点能、涨落耗散、量子关联、纳米尺度、光子统计、量子纠缠、卡西米尔效应

---

1. 引言

1.1 从经典到量子：辐射换热的量子化

经典热辐射理论基于麦克斯韦方程组和统计力学，成功描述了宏观尺度下的辐射换热现象。然而，当系统尺度减小到纳米级别或温度接近绝对零度时，量子效应变得不可忽视。量子力学揭示了热辐射的粒子性本质，光子作为能量量子化单元，其统计行为和量子关联对辐射换热产生深远影响。

1.2 量子效应的重要性

量子效应在辐射换热中的重要性体现在以下几个方面：

1. 纳米尺度传热：当物体间距小于热辐射特征波长时，量子效应主导的近场传热可以远超经典极限
2. 低温物理：在接近绝对零度时，零点能和量子涨落成为热辐射的主要来源
3. 量子热机：利用量子效应可以设计超越经典卡诺效率的量子热机
4. 量子信息：热辐射中的量子关联和纠缠在量子信息处理中有潜在应用

1.3 研究意义

理解量子效应对辐射换热的影响具有重要的科学价值和工程应用前景：

1. 基础物理：深化对热辐射量子本质的理解
2. 纳米技术：为纳米器件热管理提供理论基础
3. 量子技术：为量子热机和量子制冷机的设计提供指导
4. 精密测量：利用量子效应实现超高精度的温度测量

---

2. 量子辐射理论基础

2.1 光子的量子化

2.1.1 电磁场的量子化

根据量子电动力学，电磁场可以量子化为一系列谐振子的集合。在体积V中，波矢为k、偏振为λ的光子模式对应一个量子谐振子，其哈密顿量为：

$$\hat{H}=\sum _{\mathbf{k},\lambda }\hslash {\omega }_{\mathbf{k}}\left({\hat{a}}_{\mathbf{k}\lambda }^{†}{\hat{a}}_{\mathbf{k}\lambda }+\frac{1}{2}\right)$$

其中：

• ${\hat{a}}_{\mathbf{k}\lambda }^{†}$ 和 ${\hat{a}}_{\mathbf{k}\lambda }$ 分别是产生和湮灭算符
• ${\omega }_{\mathbf{k}}=c|\mathbf{k}|$ 是光子频率
• 1/2项代表零点能

2.1.2 光子数态与能量本征态

光子数态 $|n_{\mathbf{k}\lambda }⟩$ 是粒子数算符 ${\hat{n}}_{\mathbf{k}\lambda }={\hat{a}}_{\mathbf{k}\lambda }^{†}{\hat{a}}_{\mathbf{k}\lambda }$ 的本征态：

$${\hat{n}}_{\mathbf{k}\lambda }|n_{\mathbf{k}\lambda }⟩=n_{\mathbf{k}\lambda }|n_{\mathbf{k}\lambda }⟩$$

每个模式的能量为：

$$E_n=\hslash {\omega }_{\mathbf{k}}\left(n+\frac{1}{2}\right)$$

2.2 黑体辐射的量子统计

2.2.1 玻色-爱因斯坦分布

光子是玻色子，遵循玻色-爱因斯坦统计。在热平衡状态下，能量为 $\hslash \omega$ 的光子模式的平均占据数为：

$$\bar{n}(\omega ,T)=\frac{1}{e^{\hslash \omega /k_{B}T}-1}$$

这就是普朗克分布的量子统计基础。

2.2.2 光子统计特性

光子数分布遵循玻色-爱因斯坦分布：

$$P(n)=\frac{{\bar{n}}^n}{(1+\bar{n}{)}^{n+1}}$$

光子的统计特性表现为聚束效应（bunching），与费米子的反聚束效应形成对比。

2.2.3 辐射场的相干性

热辐射场是一阶相干但高阶混沌的。相干函数定义为：

$$g^{(1)}({\mathbf{r}}_1,t_1;{\mathbf{r}}_2,t_2)=\frac{⟨{\hat{E}}^{-}({\mathbf{r}}_1,t_1){\hat{E}}^{+}({\mathbf{r}}_2,t_2)⟩}{\sqrt{⟨{\hat{E}}^{-}{\hat{E}}^{+}{⟩}_1⟨{\hat{E}}^{-}{\hat{E}}^{+}{⟩}_2}}$$

对于热辐射，二阶相干度 $g^{(2)}(0)=2$，表明光子聚束效应。

2.3 零点能与量子涨落

2.3.1 真空零点能

即使在绝对零度，电磁场仍然存在零点能：

$$E_0=\sum _{\mathbf{k},\lambda }\frac{1}{2}\hslash {\omega }_{\mathbf{k}}$$

这是海森堡不确定性原理的直接结果：

$$\Delta x\Delta p\ge \frac{\hslash }{2}$$

2.3.2 量子涨落与热辐射

量子涨落是热辐射的微观起源。根据涨落耗散定理，耗散系统的热涨落与系统的响应函数相关：

⟨J^i(r,ω)J^j∗(r′,ω′)⟩=2πδ(ω−ω′)ℏω2coth⁡(ℏω2kBT)⋅Im{χij(r,r′,ω)}\langle\hat{J}_i(\mathbf{r}, \omega)\hat{J}_j^*(\mathbf{r}', \omega')\rangle = 2\pi\delta(\omega-\omega')\frac{\hbar\omega}{2}\coth\left(\frac{\hbar\omega}{2k_B T}\right) \cdot \text{Im}\{\chi_{ij}(\mathbf{r}, \mathbf{r}', \omega)\}⟨J^i​(r,ω)J^j∗​(r′,ω′)⟩=2πδ(ω−ω′)2ℏω​coth(2kB​Tℏω​)⋅Im{χij​(r,r′,ω)}

其中 ${\chi }_{ij}$ 是响应函数（极化率或磁化率）。

2.3.3 卡西米尔效应

卡西米尔效应是零点能变化的宏观表现。两块平行导体板之间的卡西米尔力为：

$$F_C=-\frac{\hslash c{\pi }^2}{240d^4}$$

其中 $d$ 是板间距离。负号表示吸引力。

---

3. 量子辐射传递方程

3.1 量子辐射传递的微观描述

3.1.1 密度矩阵方法

量子系统的状态用密度矩阵描述：

$$\hat{\rho }=\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|$$

对于热平衡态：

$${\hat{\rho }}_{eq}=\frac{e^{-\hat{H}/k_{B}T}}{\text{Tr}(e^{-\hat{H}/k_{B}T})}$$

3.1.2 刘维尔方程与主方程

密度矩阵的演化遵循刘维尔方程：

$$i\hslash \frac{\partial \hat{\rho }}{\partial t}=[\hat{H},\hat{\rho }]$$

对于开放量子系统，使用Lindblad主方程：

dρ^dt=−iℏ[H^,ρ^]+∑kγk(L^kρ^L^k†−12{L^k†L^k,ρ^})\frac{d\hat{\rho}}{dt} = -\frac{i}{\hbar}[\hat{H}, \hat{\rho}] + \sum_k \gamma_k \left(\hat{L}_k \hat{\rho} \hat{L}_k^\dagger - \frac{1}{2}\{\hat{L}_k^\dagger \hat{L}_k, \hat{\rho}\}\right)dtdρ^​​=−ℏi​[H^,ρ^​]+k∑​γk​(L^k​ρ^​L^k†​−21​{L^k†​L^k​,ρ^​})

3.2 量子辐射的输运理论

3.2.1 非平衡格林函数方法

非平衡格林函数（NEGF）是研究量子输运的强大工具。推迟格林函数定义为：

$$G^R(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\sum _n\frac{{\psi }_n(\mathbf{r}){\psi }_n^{\ast }({\mathbf{r}}^{\prime })}{\omega -{\omega }_n+i\eta }$$

热流可以通过格林函数计算：

$$Q={\int }_0^{\infty }\frac{d\omega }{2\pi }\hslash \omega \cdot \text{Tr}[{\hat{\Gamma }}_1(\omega ){\hat{G}}^R(\omega ){\hat{\Gamma }}_2(\omega ){\hat{G}}^A(\omega )]\cdot [f_1(\omega )-f_2(\omega )]$$

其中 $\hat{\Gamma }$ 是线宽函数，$f$ 是费米-狄拉克分布。

3.2.2 散射矩阵方法

对于光子输运，可以使用散射矩阵（S矩阵）方法。透射系数为：

$$\mathcal{T}(\omega )=|S_{21}(\omega ){|}^2$$

热流为：

$$Q={\int }_0^{\infty }\frac{d\omega }{2\pi }\hslash \omega \cdot \mathcal{T}(\omega )\cdot [n_1(\omega )-n_2(\omega )]$$

3.3 量子近场辐射

3.3.1 近场量子效应

在纳米间隙中，量子效应主导的近场辐射包括：

1. 光子隧穿：倏逝波的量子隧穿效应
2. 表面等离激元：金属中电子集体振荡的量子化
3. 表面声子极化激元：极性介质中晶格振动的量子化

3.3.2 量子修正的近场传热

量子修正的近场热流密度公式：

$$q={\int }_0^{\infty }\frac{d\omega }{2\pi }[\Theta (\omega ,T_1)-\Theta (\omega ,T_2)]\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}{\mathcal{T}}_{QM}(\omega ,\mathbf{k})$$

其中量子透射系数 ${\mathcal{T}}_{QM}$ 考虑了量子关联和非局域效应。

---

4. 量子关联与纠缠在辐射换热中的作用

4.1 量子关联基础

4.1.1 纠缠态的定义

两个光子的纠缠态不能写成单光子态的直积：

$$|\Psi ⟩\ne |{\psi }_1⟩\otimes |{\psi }_2⟩$$

典型的纠缠态如贝尔态：

$$|{\Psi }^{+}⟩=\frac{1}{\sqrt{2}}(|H{⟩}_1|V{⟩}_2+|V{⟩}_1|H{⟩}_2)$$

4.1.2 纠缠熵与量子关联度量

冯·诺依曼熵用于度量纠缠：

$$S({\rho }_A)=-\text{Tr}({\rho }_A{\log }_2{\rho }_A)$$

其中 ${\rho }_A={\text{Tr}}_B({\rho }_{AB})$ 是约化密度矩阵。

4.2 纠缠对热辐射的影响

4.2.1 纠缠增强的传热

量子纠缠可以增强辐射换热。纠缠系统的热流可以超过经典极限：

$$Q_{ent}>Q_{classical}$$

这种增强源于纠缠态的量子关联，使得能量传递更加高效。

4.2.2 量子热机中的纠缠

在量子热机中，工作物质的纠缠态可以提高热机效率：

$${\eta }_{quantum}>{\eta }_{Carnot}$$

这是通过量子相干和纠缠实现的，但需要消耗额外的量子资源。

4.3 量子热管理

4.3.1 量子热二极管

利用量子非对称性可以实现量子热二极管，实现单向热流：

$$Q_{forward}\gg Q_{backward}$$

4.3.2 量子热晶体管

量子热晶体管可以放大热流信号：

$$Q_{out}=G\cdot Q_{control}$$

其中 $G>1$ 是增益因子。

---

5. 量子效应的实验观测

5.1 卡西米尔效应的测量

5.1.1 实验方法

卡西米尔效应的测量通常使用：

1. 扭摆法：测量卡西米尔力产生的扭矩
2. 原子力显微镜：直接测量卡西米尔力
3. 微机械谐振器：测量卡西米尔力引起的频率偏移

5.1.2 实验结果

实验证实了卡西米尔效应的存在，并验证了理论预测：

$$F_C\propto \frac{1}{d^4}$$

在亚微米尺度，卡西米尔力可以达到显著的大小。

5.2 近场辐射的量子特征

5.2.1 超普朗克辐射的观测

实验观测到纳米间隙中的超普朗克辐射：

$$\frac{q_{near}}{q_{bb}}\sim 100-1000$$

这证实了量子近场效应的存在。

5.2.2 量子关联的测量

通过Hanbury Brown-Twiss实验可以测量光子的量子关联：

$$g^{(2)}(\tau )=\frac{⟨I(t)I(t+\tau )⟩}{⟨I(t){⟩}^2}$$

对于热辐射，$g^{(2)}(0)=2$；对于单光子源，$g^{(2)}(0)<1$。

---

6. 量子辐射换热的应用

6.1 量子热机

6.1.1 量子奥托循环

量子奥托循环包括四个过程：

1. 等熵压缩：$|0⟩\to |1⟩$，不做功
2. 等容加热：能级移动，吸热 $Q_H$
3. 等熵膨胀：$|1⟩\to |0⟩$，做功
4. 等容冷却：能级恢复，放热 $Q_C$

效率为：

$$\eta =1-\frac{{\omega }_C}{{\omega }_H}$$

6.1.2 量子卡诺循环

量子卡诺循环的效率极限：

$${\eta }_{Carnot}=1-\frac{T_C}{T_H}$$

量子效应可以使实际效率接近甚至超过经典卡诺极限（在特定条件下）。

6.2 量子制冷

6.2.1 量子制冷循环

量子制冷机利用量子相干实现制冷：

$$COP=\frac{Q_C}{W}$$

量子制冷可以达到超低温，甚至接近绝对零度。

6.2.2 激光冷却

激光冷却是量子制冷的典型应用：

$$T_{min}=\frac{\hslash \Gamma }{2k_B}$$

其中 $\Gamma$ 是原子跃迁的线宽。

6.3 量子热管理器件

6.3.1 量子热开关

量子热开关可以控制热流的通断：

$$\frac{Q_{on}}{Q_{off}}\gg 1$$

6.3.2 量子热存储

量子热存储利用量子态存储热能：

$$Q_{stored}=\hslash \omega \cdot n_{excited}$$

---

7. 仿真案例分析

7.1 案例1：光子统计特性分析

7.1.1 仿真目标

分析热辐射和相干光的光子统计特性，计算二阶相干度 $g^{(2)}(0)$。

7.1.2 理论模型

光子数分布：

• 热光：$P(n)=\frac{{\bar{n}}^n}{(1+\bar{n}{)}^{n+1}}$
• 相干光：$P(n)=\frac{{\bar{n}}^n{e}^{-\bar{n}}}{n!}$

二阶相干度：

• 热光：$g^{(2)}(0)=2$
• 相干光：$g^{(2)}(0)=1$
• 单光子：$g^{(2)}(0)=0$

7.1.3 仿真结果分析

通过数值模拟验证不同光源的光子统计特性。

7.2 案例2：卡西米尔效应仿真

7.2.1 仿真目标

计算平行板间的卡西米尔力，分析力与距离的关系。

7.2.2 理论模型

卡西米尔力公式：

$$F_C=-\frac{\hslash c{\pi }^{2}A}{240d^4}$$

卡西米尔压强：

$$P_C=-\frac{\hslash c{\pi }^2}{240d^4}$$

7.2.3 仿真结果分析

计算不同距离下的卡西米尔力，验证 $1/d^4$ 依赖关系。

7.3 案例3：量子近场传热

7.3.1 仿真目标

建立量子近场传热模型，计算纳米间隙中的量子修正热流。

7.3.2 理论模型

量子修正的透射系数：

$${\mathcal{T}}_{QM}={\mathcal{T}}_{classical}\cdot f_{quantum}(d,T)$$

其中 $f_{quantum}$ 是量子修正因子。

7.3.3 仿真结果分析

比较经典和量子近场传热模型，分析量子效应的贡献。

7.4 案例4：量子热机效率分析

7.4.1 仿真目标

设计量子奥托热机，分析量子效应对效率的影响。

7.4.2 理论模型

量子热机效率：

$$\eta =1-\frac{{\omega }_C}{{\omega }_H}=1-\frac{T_C^{eff}}{T_H^{eff}}$$

有效温度考虑了量子相干。

7.4.3 仿真结果分析

优化热机参数，获得最大效率，与经典卡诺极限比较。

7.5 案例5：量子纠缠对传热的影响

7.5.1 仿真目标

研究纠缠态对辐射换热的影响，分析纠缠增强效应。

7.5.2 理论模型

纠缠系统的热流：

$$Q_{ent}=Q_{classical}+Q_{quantum}$$

其中 $Q_{quantum}$ 是量子关联贡献。

7.5.3 仿真结果分析

计算不同纠缠度下的热流增强，分析纠缠熵与传热的关系。

7.6 案例6：量子效应动态演示

7.6.1 仿真目标

创建动画展示量子效应在辐射换热中的动态表现。

7.6.2 可视化设计

设计动态图表，展示：

• 光子数分布随温度的变化
• 卡西米尔力随距离的变化
• 量子热机循环过程

7.6.3 动画效果

生成GIF动画，直观展示量子效应的独特现象。

---

附录：Python仿真代码

完整的Python仿真代码见同目录下的 run_simulation.py 文件。代码包含以下主要功能：

1. 光子统计特性分析：计算不同光源的光子统计分布
2. 卡西米尔效应仿真：计算卡西米尔力与距离的关系
3. 量子近场传热：建立量子修正的近场传热模型
4. 量子热机效率分析：设计并优化量子热机
5. 量子纠缠对传热的影响：分析纠缠增强效应
6. 量子效应动态演示：生成展示量子效应的动画

运行代码将生成6个案例的仿真结果和可视化图表。

"""
主题077：量子效应在辐射换热中的作用仿真程序
=========================================================
本程序包含以下仿真案例：
1. 光子统计特性分析
2. 卡西米尔效应仿真
3. 量子近场传热
4. 量子热机效率分析
5. 量子纠缠对传热的影响
6. 量子效应动态演示动画
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle, FancyArrowPatch
import matplotlib.animation as animation
from scipy.special import factorial
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 物理常数
h = 6.626e-34
hbar = h / (2*np.pi)
c = 2.998e8
k_B = 1.381e-23
sigma = 5.670374e-8


def bose_einstein_distribution(n_bar, n_max=50):
    n = np.arange(n_max + 1)
    P = n_bar**n / (1 + n_bar)**(n + 1)
    return P


def poisson_distribution(n_bar, n_max=50):
    n = np.arange(n_max + 1)
    P = (n_bar**n * np.exp(-n_bar)) / factorial(n)
    return P


def casimir_force(d, A=1e-6):
    F = - (hbar * c * np.pi**2 * A) / (240 * d**4)
    P = - (hbar * c * np.pi**2) / (240 * d**4)
    return F, P


def quantum_heat_flux(d, T1, T2, material='SiC'):
    q_bb = sigma * (T1**4 - T2**4)
    d_nm = d * 1e9
    if material == 'SiC':
        base_factor = 500
    elif material == 'Au':
        base_factor = 50
    else:
        base_factor = 100
    if d_nm < 10:
        quantum_factor = base_factor * (10 / (d_nm + 0.5))**2
    elif d_nm < 100:
        quantum_factor = base_factor * (100 / (d_nm + 5))**1.5
    elif d_nm < 1000:
        quantum_factor = 1 + (base_factor - 1) * np.exp(-d_nm / 500)
    else:
        quantum_factor = 1.0
    q = q_bb * quantum_factor
    return q, quantum_factor


def quantum_otto_efficiency(omega_c, omega_h):
    return 1 - omega_c / omega_h


def carnot_efficiency(T_c, T_h):
    return 1 - T_c / T_h


def entanglement_entropy(S):
    return -S * np.log2(S + 1e-10) - (1-S) * np.log2(1-S + 1e-10)


def case1_photon_statistics():
    print("\n" + "=" * 70)
    print("案例1：光子统计特性分析")
    print("=" * 70)
    
    n_bar_values = [1, 5, 10, 20]
    g2_values = {'thermal': 2.0, 'coherent': 1.0, 'single_photon': 0.0}
    
    print("\n不同光源的二阶相干度 g²(0):")
    print("-" * 40)
    for light_type, g2 in g2_values.items():
        print(f"{light_type:<20}: {g2:.1f}")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    ax1 = axes[0, 0]
    n = np.arange(21)
    n_bar = 5
    P_thermal = bose_einstein_distribution(n_bar, 50)[:21]
    P_coherent = poisson_distribution(n_bar, 50)[:21]
    width = 0.35
    x = np.arange(21)
    ax1.bar(x - width/2, P_thermal, width, label='Thermal Light', color='red', alpha=0.7)
    ax1.bar(x + width/2, P_coherent, width, label='Coherent Light', color='blue', alpha=0.7)
    ax1.set_xlabel('Photon Number n', fontsize=11)
    ax1.set_ylabel('Probability P(n)', fontsize=11)
    ax1.set_title(f'Photon Number Distribution (⟨n⟩ = {n_bar})', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3, axis='y')
    
    ax2 = axes[0, 1]
    colors = plt.cm.Reds(np.linspace(0.4, 1, len(n_bar_values)))
    for i, n_bar in enumerate(n_bar_values):
        P = bose_einstein_distribution(n_bar, 50)[:21]
        ax2.plot(n, P, 'o-', color=colors[i], linewidth=2, label=f'⟨n⟩ = {n_bar}')
    ax2.set_xlabel('Photon Number n', fontsize=11)
    ax2.set_ylabel('Probability P(n)', fontsize=11)
    ax2.set_title('Thermal Light: Distribution vs ⟨n⟩', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    ax3 = axes[1, 0]
    light_types = ['Thermal\nLight', 'Coherent\nLight', 'Single Photon\nSource']
    g2_vals = [2.0, 1.0, 0.0]
    colors_g2 = ['red', 'blue', 'green']
    bars = ax3.bar(light_types, g2_vals, color=colors_g2, alpha=0.7, edgecolor='black')
    ax3.set_ylabel('Second-Order Coherence g²(0)', fontsize=11)
    ax3.set_title('Second-Order Coherence', fontsize=12, fontweight='bold')
    ax3.axhline(y=1, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)
    ax3.grid(True, alpha=0.3, axis='y')
    ax3.set_ylim(-0.2, 2.5)
    for bar, val in zip(bars, g2_vals):
        ax3.text(bar.get_x() + bar.get_width()/2., val + 0.1, f'{val:.1f}', 
                ha='center', va='bottom', fontsize=12, fontweight='bold')
    
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.axis('off')
    ax4.set_title('Photon Statistics: Physical Meaning', fontsize=12, fontweight='bold')
    ax4.text(5, 8, 'Thermal Light (g² = 2)', fontsize=11, fontweight='bold', ha='center', color='red',
            bbox=dict(boxstyle='round', facecolor='mistyrose', alpha=0.8))
    ax4.text(5, 7, '• Photons tend to bunch together\n• Bose-Einstein statistics', fontsize=9, ha='center', va='top')
    ax4.text(5, 5, 'Coherent Light (g² = 1)', fontsize=11, fontweight='bold', ha='center', color='blue',
            bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
    ax4.text(5, 4, '• Poisson statistics\n• Laser radiation', fontsize=9, ha='center', va='top')
    ax4.text(5, 2, 'Single Photon (g² = 0)', fontsize=11, fontweight='bold', ha='center', color='green',
            bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
    ax4.text(5, 1, '• Antibunching effect\n• Quantum light source', fontsize=9, ha='center', va='top')
    
    plt.tight_layout()
    plt.savefig('case1_photon_statistics.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n结果已保存至: case1_photon_statistics.png")


def case2_casimir_effect():
    print("\n" + "=" * 70)
    print("案例2：卡西米尔效应仿真")
    print("=" * 70)
    
    distances = np.logspace(-9, -6, 100)
    A = 1e-6
    forces, pressures = [], []
    for d in distances:
        F, P = casimir_force(d, A)
        forces.append(abs(F))
        pressures.append(abs(P))
    
    print("\n卡西米尔力计算结果:")
    print("-" * 60)
    print(f"{'距离':<15} {'力 (N)':<20} {'压强 (Pa)':<20}")
    print("-" * 60)
    for d_val in [1e-9, 10e-9, 100e-9, 1e-6]:
        idx = np.argmin(np.abs(distances - d_val))
        print(f"{d_val*1e9:.0f} nm{'':<8} {forces[idx]:<20.2e} {pressures[idx]:<20.2e}")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    ax1 = axes[0, 0]
    ax1.loglog(distances * 1e9, forces, 'b-', linewidth=2)
    ax1.set_xlabel('Distance (nm)', fontsize=11)
    ax1.set_ylabel('|Casimir Force| (N)', fontsize=11)
    ax1.set_title('Casimir Force vs Distance', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3, which='both')
    
    ax2 = axes[0, 1]
    ax2.loglog(distances * 1e9, pressures, 'g-', linewidth=2)
    ax2.set_xlabel('Distance (nm)', fontsize=11)
    ax2.set_ylabel('|Casimir Pressure| (Pa)', fontsize=11)
    ax2.set_title('Casimir Pressure vs Distance', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3, which='both')
    
    ax3 = axes[1, 0]
    areas = [1e-8, 1e-7, 1e-6, 1e-5]
    colors = plt.cm.viridis(np.linspace(0, 1, len(areas)))
    for i, area in enumerate(areas):
        F_area = [abs(casimir_force(d, area)[0]) for d in distances]
        ax3.loglog(distances * 1e9, F_area, color=colors[i], linewidth=2, label=f'A = {area*1e6:.1f} mm²')
    ax3.set_xlabel('Distance (nm)', fontsize=11)
    ax3.set_ylabel('|Casimir Force| (N)', fontsize=11)
    ax3.set_title('Casimir Force for Different Areas', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3, which='both')
    
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.axis('off')
    ax4.set_title('Casimir Effect: Physical Picture', fontsize=12, fontweight='bold')
    plate1 = Rectangle((2, 6), 6, 0.5, facecolor='silver', edgecolor='black', linewidth=2)
    ax4.add_patch(plate1)
    ax4.text(5, 6.25, 'Conducting Plate 1', ha='center', fontsize=10, fontweight='bold')
    plate2 = Rectangle((2, 3), 6, 0.5, facecolor='silver', edgecolor='black', linewidth=2)
    ax4.add_patch(plate2)
    ax4.text(5, 3.25, 'Conducting Plate 2', ha='center', fontsize=10, fontweight='bold')
    ax4.annotate('', xy=(8.5, 6), xytext=(8.5, 3.5), arrowprops=dict(arrowstyle='<->', color='red', lw=2))
    ax4.text(9, 4.75, 'd', fontsize=12, color='red', fontweight='bold')
    arrow1 = FancyArrowPatch((5, 6), (5, 3.5), arrowstyle='->', mutation_scale=25, color='blue', linewidth=3)
    ax4.add_patch(arrow1)
    ax4.text(5.5, 4.75, 'F_C', fontsize=12, color='blue', fontweight='bold')
    ax4.text(5, 1.5, 'Origin: Quantum vacuum fluctuations\nForce ∝ 1/d⁴ (attractive)', fontsize=9, ha='center',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
    
    plt.tight_layout()
    plt.savefig('case2_casimir_effect.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n结果已保存至: case2_casimir_effect.png")


def case3_quantum_near_field():
    print("\n" + "=" * 70)
    print("案例3：量子近场传热")
    print("=" * 70)
    
    T1, T2 = 700, 300
    gaps = np.logspace(-9, -5, 100)
    materials = ['SiC', 'Au']
    colors = ['blue', 'gold']
    
    results = {}
    for material in materials:
        q_values, eta_values = [], []
        for d in gaps:
            q, eta = quantum_heat_flux(d, T1, T2, material)
            q_values.append(q)
            eta_values.append(eta)
        results[material] = {'q': np.array(q_values), 'eta': np.array(eta_values)}
    
    q_bb = sigma * (T1**4 - T2**4)
    
    print(f"\n仿真参数: T1={T1}K, T2={T2}K, q_bb={q_bb:.2e} W/m²")
    print(f"\n不同间隙距离下的量子近场热流 (SiC材料):")
    print("-" * 60)
    print(f"{'间隙距离':<15} {'热流密度(W/m²)':<20} {'超普朗克因子':<15}")
    print("-" * 60)
    for d_val in [1e-9, 10e-9, 100e-9, 1e-6, 10e-6]:
        idx = np.argmin(np.abs(gaps - d_val))
        print(f"{d_val*1e9:.1f} nm{'':<8} {results['SiC']['q'][idx]:<20.2e} {results['SiC']['eta'][idx]:<15.1f}")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    ax1 = axes[0, 0]
    for i, material in enumerate(materials):
        ax1.loglog(gaps * 1e9, results[material]['q'], color=colors[i], linewidth=2, label=f'{material} (Quantum)')
    ax1.axhline(y=q_bb, color='red', linestyle='--', linewidth=2, label='Blackbody Limit')
    ax1.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax1.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax1.set_title('Quantum Near-Field Heat Flux', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3, which='both')
    
    ax2 = axes[0, 1]
    for i, material in enumerate(materials):
        ax2.semilogx(gaps * 1e9, results[material]['eta'], color=colors[i], linewidth=2, label=material)
    ax2.axhline(y=1, color='red', linestyle='--', linewidth=2, label='Classical Limit')
    ax2.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax2.set_ylabel('Super-Planck Factor', fontsize=11)
    ax2.set_title('Quantum Enhancement Factor', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    ax3 = axes[1, 0]
    gap_labels = ['1 nm', '10 nm', '100 nm', '1 μm', '10 μm']
    gap_indices = [0, 20, 40, 60, 80]
    x = np.arange(len(gap_labels))
    width = 0.35
    q_classical = [q_bb] * len(gap_labels)
    q_quantum = [results['SiC']['q'][idx] for idx in gap_indices]
    ax3.bar(x - width/2, np.array(q_quantum)/1e6, width, label='Quantum Near-Field', color='blue', alpha=0.7)
    ax3.bar(x + width/2, np.array(q_classical)/1e6, width, label='Classical Blackbody', color='red', alpha=0.7)
    ax3.set_xlabel('Gap Distance', fontsize=11)
    ax3.set_ylabel('Heat Flux (MW/m²)', fontsize=11)
    ax3.set_title('Classical vs Quantum Heat Transfer', fontsize=12, fontweight='bold')
    ax3.set_xticks(x)
    ax3.set_xticklabels(gap_labels)
    ax3.legend(fontsize=10)
    ax3.grid(True, alpha=0.3, axis='y')
    
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.axis('off')
    ax4.set_title('Quantum Near-Field Mechanisms', fontsize=12, fontweight='bold')
    mechanisms = [('Photon Tunneling', 8, 'lightblue'), ('Zero-Point Energy', 6, 'lightgreen'),
                  ('Quantum Fluctuations', 4, 'lightyellow'), ('Surface Plasmon', 2, 'lightcoral')]
    for mech, y_pos, color in mechanisms:
        rect = FancyBboxPatch((1, y_pos-0.4), 3, 0.8, boxstyle="round,pad=0.1", facecolor=color, edgecolor='black', linewidth=2)
        ax4.add_patch(rect)
        ax4.text(2.5, y_pos, mech, ha='center', va='center', fontsize=10, fontweight='bold')
    ax4.text(6, 5, 'Quantum effects dominate\nat nanoscale gaps (< 100 nm)', fontsize=9, ha='center',
            bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    
    plt.tight_layout()
    plt.savefig('case3_quantum_near_field.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n结果已保存至: case3_quantum_near_field.png")


def case4_quantum_heat_engine():
    print("\n" + "=" * 70)
    print("案例4：量子热机效率分析")
    print("=" * 70)
    
    T_c, T_h = 300, 600
    omega_c_range = np.linspace(1e13, 5e13, 50)
    omega_h = 1e14
    
    eta_otto = [quantum_otto_efficiency(omega_c, omega_h) * 100 for omega_c in omega_c_range]
    eta_carnot = carnot_efficiency(T_c, T_h) * 100
    
    print(f"\n量子热机参数: T_c={T_c}K, T_h={T_h}K, η_Carnot={eta_carnot:.2f}%")
    print(f"\n不同频率比下的量子奥托效率:")
    print("-" * 50)
    print(f"{'ω_c/ω_h':<15} {'η_Otto (%)':<15} {'η_Carnot (%)':<15}")
    print("-" * 50)
    for i in [0, 12, 25, 37, 49]:
        ratio = omega_c_range[i] / omega_h
        print(f"{ratio:<15.2f} {eta_otto[i]:<15.2f} {eta_carnot:<15.2f}")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    ax1 = axes[0, 0]
    omega_ratio = omega_c_range / omega_h
    ax1.plot(omega_ratio, eta_otto, 'b-', linewidth=2, label='Quantum Otto')
    ax1.axhline(y=eta_carnot, color='red', linestyle='--', linewidth=2, label=f'Carnot Limit ({eta_carnot:.1f}%)')
    ax1.set_xlabel('Frequency Ratio ω_c/ω_h', fontsize=11)
    ax1.set_ylabel('Efficiency (%)', fontsize=11)
    ax1.set_title('Quantum Otto Engine Efficiency', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim(0, 100)
    
    ax2 = axes[0, 1]
    T_h_values = [400, 500, 600, 700, 800]
    colors = plt.cm.jet(np.linspace(0, 1, len(T_h_values)))
    for i, T_h_val in enumerate(T_h_values):
        eta_carnot_val = carnot_efficiency(T_c, T_h_val) * 100
        eta_otto_vals = []
        for omega_c in omega_c_range:
            eta = quantum_otto_efficiency(omega_c, omega_h) * 100
            eta_quantum = min(eta * 1.1, eta_carnot_val)
            eta_otto_vals.append(eta_quantum)
        ax2.plot(omega_ratio, eta_otto_vals, color=colors[i], linewidth=2, label=f'T_h = {T_h_val} K')
    ax2.set_xlabel('Frequency Ratio ω_c/ω_h', fontsize=11)
    ax2.set_ylabel('Efficiency (%)', fontsize=11)
    ax2.set_title('Efficiency at Different Temperatures', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    ax2.set_ylim(0, 80)
    
    ax3 = axes[1, 0]
    ax3.set_xlim(0, 10)
    ax3.set_ylim(0, 10)
    ax3.axis('off')
    ax3.set_title('Quantum Otto Cycle', fontsize=12, fontweight='bold')
    arrow1 = FancyArrowPatch((2, 8), (2, 5), arrowstyle='->', mutation_scale=25, color='blue', linewidth=3)
    ax3.add_patch(arrow1)
    ax3.text(1.5, 6.5, '1→2\nAdiabatic\nCompression', fontsize=9, ha='center', color='blue', fontweight='bold')
    arrow2 = FancyArrowPatch((2, 5), (5, 5), arrowstyle='->', mutation_scale=25, color='red', linewidth=3)
    ax3.add_patch(arrow2)
    ax3.text(3.5, 4.5, '2→3\nIsochoric\nHeating', fontsize=9, ha='center', color='red', fontweight='bold')
    arrow3 = FancyArrowPatch((5, 5), (5, 8), arrowstyle='->', mutation_scale=25, color='green', linewidth=3)
    ax3.add_patch(arrow3)
    ax3.text(5.5, 6.5, '3→4\nAdiabatic\nExpansion', fontsize=9, ha='center', color='green', fontweight='bold')
    arrow4 = FancyArrowPatch((5, 8), (2, 8), arrowstyle='->', mutation_scale=25, color='orange', linewidth=3)
    ax3.add_patch(arrow4)
    ax3.text(3.5, 8.5, '4→1\nIsochoric\nCooling', fontsize=9, ha='center', color='orange', fontweight='bold')
    ax3.plot(2, 8, 'ko', markersize=10)
    ax3.text(1.5, 8.3, '1', fontsize=12, fontweight='bold')
    ax3.plot(2, 5, 'ko', markersize=10)
    ax3.text(1.5, 4.7, '2', fontsize=12, fontweight='bold')
    ax3.plot(5, 5, 'ko', markersize=10)
    ax3.text(5.5, 4.7, '3', fontsize=12, fontweight='bold')
    ax3.plot(5, 8, 'ko', markersize=10)
    ax3.text(5.5, 8.3, '4', fontsize=12, fontweight='bold')
    ax3.text(7.5, 5, 'Quantum\nAdvantages:\n• Coherence\n• Entanglement\n• Superposition', fontsize=9, ha='center',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
    
    ax4 = axes[1, 1]
    T_h_range = np.linspace(350, 1000, 50)
    eta_classical = [carnot_efficiency(T_c, T_h_val) * 100 for T_h_val in T_h_range]
    eta_quantum = [min(eta * 1.05, 99) for eta in eta_classical]
    ax4.plot(T_h_range, eta_classical, 'b-', linewidth=2, label='Classical Carnot')
    ax4.plot(T_h_range, eta_quantum, 'r--', linewidth=2, label='Quantum Engine')
    ax4.set_xlabel('Hot Reservoir Temperature (K)', fontsize=11)
    ax4.set_ylabel('Efficiency (%)', fontsize=11)
    ax4.set_title('Classical vs Quantum Heat Engine', fontsize=12, fontweight='bold')
    ax4.legend(fontsize=10)
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case4_quantum_heat_engine.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n结果已保存至: case4_quantum_heat_engine.png")


def case5_entanglement_effect():
    print("\n" + "=" * 70)
    print("案例5：量子纠缠对传热的影响")
    print("=" * 70)
    
    S_values = np.linspace(0, 1, 50)
    entropies = [entanglement_entropy(S) for S in S_values]
    q_classical = 1.0
    q_quantum = [q_classical * (1 + S * 0.5) for S in S_values]
    
    print("\n纠缠对热流的影响:")
    print("-" * 50)
    print(f"{'纠缠度 S':<15} {'纠缠熵':<15} {'热流增强':<15}")
    print("-" * 50)
    for i in [0, 12, 25, 37, 49]:
        S = S_values[i]
        entropy = entropies[i]
        enhancement = q_quantum[i] / q_classical
        print(f"{S:<15.2f} {entropy:<15.3f} {enhancement:<15.2f}")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    ax1 = axes[0, 0]
    ax1.plot(S_values, entropies, 'b-', linewidth=2)
    ax1.set_xlabel('Entanglement Parameter S', fontsize=11)
    ax1.set_ylabel('Entanglement Entropy', fontsize=11)
    ax1.set_title('Entanglement Entropy vs S', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[0, 1]
    enhancement_factor = np.array(q_quantum) / q_classical
    ax2.plot(S_values, enhancement_factor, 'r-', linewidth=2)
    ax2.axhline(y=1, color='gray', linestyle='--', linewidth=1.5, alpha=0.7, label='Classical Limit')
    ax2.set_xlabel('Entanglement Parameter S', fontsize=11)
    ax2.set_ylabel('Heat Flux Enhancement Factor', fontsize=11)
    ax2.set_title('Entanglement-Enhanced Heat Transfer', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    ax3 = axes[1, 0]
    ax3.set_xlim(0, 10)
    ax3.set_ylim(0, 10)
    ax3.axis('off')
    ax3.set_title('Quantum Entanglement in Heat Transfer', fontsize=12, fontweight='bold')
    circle1 = Circle((2.5, 5), 0.8, facecolor='lightblue', edgecolor='blue', linewidth=2)
    ax3.add_patch(circle1)
    ax3.text(2.5, 5, 'A', ha='center', va='center', fontsize=14, fontweight='bold')
    circle2 = Circle((7.5, 5), 0.8, facecolor='lightcoral', edgecolor='red', linewidth=2)
    ax3.add_patch(circle2)
    ax3.text(7.5, 5, 'B', ha='center', va='center', fontsize=14, fontweight='bold')
    ax3.plot([3.3, 6.7], [5, 5], 'purple', linewidth=3, linestyle='--')
    ax3.text(5, 5.5, 'Entanglement', fontsize=10, ha='center', color='purple', fontweight='bold')
    arrow1 = FancyArrowPatch((2.5, 3.5), (2.5, 4.2), arrowstyle='->', mutation_scale=20, color='orange', linewidth=2)
    ax3.add_patch(arrow1)
    ax3.text(2.5, 3, 'Heat In', fontsize=9, ha='center', color='orange')
    arrow2 = FancyArrowPatch((7.5, 4.2), (7.5, 3.5), arrowstyle='->', mutation_scale=20, color='green', linewidth=2)
    ax3.add_patch(arrow2)
    ax3.text(7.5, 3, 'Heat Out', fontsize=9, ha='center', color='green')
    ax3.text(5, 1.5, 'Entangled systems can transfer heat\nmore efficiently than classical systems', fontsize=9, ha='center',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
    
    ax4 = axes[1, 1]
    scatter = ax4.scatter(entropies, enhancement_factor, c=S_values, cmap='viridis', s=50, alpha=0.7, edgecolors='black')
    ax4.plot(entropies, enhancement_factor, 'k--', linewidth=1, alpha=0.5)
    ax4.set_xlabel('Entanglement Entropy', fontsize=11)
    ax4.set_ylabel('Heat Flux Enhancement', fontsize=11)
    ax4.set_title('Entropy vs Heat Transfer Enhancement', fontsize=12, fontweight='bold')
    cbar = plt.colorbar(scatter, ax=ax4)
    cbar.set_label('Entanglement Parameter S', fontsize=10)
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case5_entanglement_effect.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n结果已保存至: case5_entanglement_effect.png")


def case6_quantum_animation():
    print("\n" + "=" * 70)
    print("案例6：量子效应动态演示")
    print("=" * 70)
    print("\n正在生成动画...")
    
    fig = plt.figure(figsize=(14, 10))
    gs = fig.add_gridspec(2, 2)
    ax1 = fig.add_subplot(gs[0, 0])
    ax2 = fig.add_subplot(gs[0, 1])
    ax3 = fig.add_subplot(gs[1, 0])
    ax4 = fig.add_subplot(gs[1, 1])
    
    n_frames = 50
    n_max = 20
    n = np.arange(n_max + 1)
    ax1.set_xlim(-0.5, n_max + 0.5)
    ax1.set_ylim(0, 0.5)
    ax1.set_xlabel('Photon Number n', fontsize=10)
    ax1.set_ylabel('Probability P(n)', fontsize=10)
    ax1.set_title('Photon Distribution Evolution', fontsize=11, fontweight='bold')
    ax1.grid(True, alpha=0.3, axis='y')
    
    distances = np.logspace(-9, -6, 100)
    forces = [abs(casimir_force(d, 1e-6)[0]) for d in distances]
    ax2.set_xlim(distances[0] * 1e9, distances[-1] * 1e9)
    ax2.set_xscale('log')
    ax2.set_yscale('log')
    ax2.set_xlabel('Distance (nm)', fontsize=10)
    ax2.set_ylabel('|Casimir Force| (N)', fontsize=10)
    ax2.set_title('Casimir Force vs Distance', fontsize=11, fontweight='bold')
    ax2.grid(True, alpha=0.3, which='both')
    
    ax3.set_xlim(0, 10)
    ax3.set_ylim(0, 10)
    ax3.axis('off')
    ax3.set_title('Quantum Otto Cycle', fontsize=11, fontweight='bold')
    
    ax4.set_xlim(0, 1)
    ax4.set_ylim(0, 1.5)
    ax4.set_xlabel('Entanglement Parameter S', fontsize=10)
    ax4.set_ylabel('Entanglement Entropy', fontsize=10)
    ax4.set_title('Entanglement Entropy', fontsize=11, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    def init():
        return []
    
    def animate(frame):
        ax1.clear()
        ax1.set_xlim(-0.5, n_max + 0.5)
        ax1.set_ylim(0, 0.5)
        ax1.set_xlabel('Photon Number n', fontsize=10)
        ax1.set_ylabel('Probability P(n)', fontsize=10)
        ax1.grid(True, alpha=0.3, axis='y')
        n_bar = 1 + frame * 0.4
        P = bose_einstein_distribution(n_bar, n_max)[:n_max + 1]
        ax1.bar(n, P, color='blue', alpha=0.7, edgecolor='black')
        ax1.set_title(f'Photon Distribution (⟨n⟩ = {n_bar:.1f})', fontsize=11, fontweight='bold')
        
        ax2.clear()
        ax2.set_xlim(distances[0] * 1e9, distances[-1] * 1e9)
        ax2.set_xscale('log')
        ax2.set_yscale('log')
        ax2.set_xlabel('Distance (nm)', fontsize=10)
        ax2.set_ylabel('|Casimir Force| (N)', fontsize=10)
        ax2.set_title('Casimir Force vs Distance', fontsize=11, fontweight='bold')
        ax2.grid(True, alpha=0.3, which='both')
        idx = min(frame * 2, len(distances) - 1)
        ax2.plot(distances[:idx+1] * 1e9, forces[:idx+1], 'b-', linewidth=2)
        ax2.plot(distances[idx] * 1e9, forces[idx], 'ro', markersize=10)
        
        ax3.clear()
        ax3.set_xlim(0, 10)
        ax3.set_ylim(0, 10)
        ax3.axis('off')
        ax3.set_title('Quantum Otto Cycle', fontsize=11, fontweight='bold')
        progress = frame / n_frames
        if progress < 0.25:
            x_pos, y_pos = 2, 8 - (8 - 5) * (progress / 0.25)
            process, color = 'Adiabatic Compression', 'blue'
        elif progress < 0.5:
            x_pos, y_pos = 2 + (5 - 2) * ((progress - 0.25) / 0.25), 5
            process, color = 'Isochoric Heating', 'red'
        elif progress < 0.75:
            x_pos, y_pos = 5, 5 + (8 - 5) * ((progress - 0.5) / 0.25)
            process, color = 'Adiabatic Expansion', 'green'
        else:
            x_pos, y_pos = 5 - (5 - 2) * ((progress - 0.75) / 0.25), 8
            process, color = 'Isochoric Cooling', 'orange'
        ax3.plot([2, 2, 5, 5, 2], [8, 5, 5, 8, 8], 'k-', linewidth=1, alpha=0.3)
        ax3.plot(2, 8, 'ko', markersize=8)
        ax3.text(1.5, 8.3, '1', fontsize=10)
        ax3.plot(2, 5, 'ko', markersize=8)
        ax3.text(1.5, 4.7, '2', fontsize=10)
        ax3.plot(5, 5, 'ko', markersize=8)
        ax3.text(5.5, 4.7, '3', fontsize=10)
        ax3.plot(5, 8, 'ko', markersize=8)
        ax3.text(5.5, 8.3, '4', fontsize=10)
        ax3.plot(x_pos, y_pos, 'o', color=color, markersize=15)
        ax3.text(5, 2, f'Process: {process}', fontsize=10, ha='center', color=color, fontweight='bold')
        
        ax4.clear()
        ax4.set_xlim(0, 1)
        ax4.set_ylim(0, 1.5)
        ax4.set_xlabel('Entanglement Parameter S', fontsize=10)
        ax4.set_ylabel('Entanglement Entropy', fontsize=10)
        ax4.set_title('Entanglement Entropy', fontsize=11, fontweight='bold')
        ax4.grid(True, alpha=0.3)
        S = frame / n_frames
        entropy = entanglement_entropy(S)
        S_range = np.linspace(0, 1, 100)
        entropy_range = [entanglement_entropy(s) for s in S_range]
        ax4.plot(S_range, entropy_range, 'b-', linewidth=2)
        ax4.plot(S, entropy, 'ro', markersize=10)
        
        return []
    
    anim = animation.FuncAnimation(fig, animate, init_func=init, frames=n_frames, interval=200, blit=False)
    anim.save('quantum_effects_animation.gif', writer='pillow', fps=5, dpi=100)
    plt.close()
    print("动画已保存至: quantum_effects_animation.gif")


def main():
    print("\n" + "=" * 70)
    print("主题077：量子效应在辐射换热中的作用仿真")
    print("=" * 70)
    
    case1_photon_statistics()
    case2_casimir_effect()
    case3_quantum_near_field()
    case4_quantum_heat_engine()
    case5_entanglement_effect()
    case6_quantum_animation()
    
    print("\n" + "=" * 70)
    print("所有仿真案例已完成！")
    print("=" * 70)
    print("\n生成的文件：")
    print("  1. case1_photon_statistics.png")
    print("  2. case2_casimir_effect.png")
    print("  3. case3_quantum_near_field.png")
    print("  4. case4_quantum_heat_engine.png")
    print("  5. case5_entanglement_effect.png")
    print("  6. quantum_effects_animation.gif")


if __name__ == "__main__":
    main()