主题063：各向异性散射介质辐射

摘要

各向异性散射是参与性介质辐射传递中的重要物理现象，广泛存在于大气、生物组织、工业悬浮液等实际系统中。本主题深入研究各向异性散射对辐射传递的影响机制，系统讲解散射相函数的理论基础、数学描述和数值处理方法。通过建立一般散射相函数的辐射传递模型，分析不对称因子、前向散射和后向散射对辐射场分布的影响规律。结合四个完整的Python仿真案例，包括大气气溶胶散射、生物组织光学、工业颗粒系统和相函数反演，帮助读者全面掌握各向异性散射介质的辐射换热分析方法。

关键词

各向异性散射、相函数、不对称因子、前向散射、后向散射、Henyey-Greenstein相函数、Legendre展开、辐射传递方程

---

1. 引言

1.1 各向异性散射的物理意义

在辐射传递过程中，当光子与介质中的颗粒或分子相互作用时，会发生散射现象。散射的方向分布特性直接决定了辐射能在介质中的传播路径和能量分布。根据散射方向分布的不同，可以将散射分为：

各向同性散射：散射光在所有方向上均匀分布，这是一种理想化的简化模型。

各向异性散射：散射光在不同方向上具有不同的强度分布，这是自然界和工程实际中更普遍的情况。

各向异性散射的物理机制包括：

• 瑞利散射：当散射颗粒尺寸远小于波长时发生，如大气分子对阳光的散射
• 米氏散射：当散射颗粒尺寸与波长相当或更大时发生，如气溶胶、水滴的散射
• 几何光学散射：当颗粒尺寸远大于波长时，可以用几何光学近似处理

1.2 各向异性散射的工程应用背景

各向异性散射现象在众多工程领域具有重要应用：

大气科学：气溶胶和云滴的散射特性直接影响地球辐射平衡和气候变化。前向散射主导的气溶胶会增强地表太阳辐射，而后向散射则会增加大气反照率。

生物医学光学：人体组织中的细胞、胶原蛋白等结构产生强烈的各向异性散射，前向散射特性使得光能够在组织中传播较远距离，这是光学成像和光动力治疗的基础。

工业过程：燃烧系统中的煤粉、喷雾干燥中的液滴、结晶过程中的颗粒等都表现出明显的各向异性散射特性，影响反应器内的辐射传热效率。

遥感技术：卫星遥感依赖于大气和地表散射特性的精确建模，各向异性散射的准确描述是反演算法的基础。

1.3 本主题的学习目标

通过本主题的学习，读者将能够：

1. 理解各向异性散射的物理机制和数学描述方法
2. 掌握常用散射相函数模型及其适用范围
3. 建立各向异性散射介质的辐射传递方程
4. 运用数值方法求解各向异性散射辐射问题
5. 分析不对称因子对辐射传递的影响规律
6. 通过Python仿真实现各向异性散射辐射计算

---

2. 散射理论基础

2.1 散射的基本概念

2.1.1 散射的定义与分类

散射是指电磁波与物质相互作用时，改变传播方向而不改变频率的现象。根据散射前后能量和动量的变化，可以分为：

弹性散射：散射前后光子能量不变，仅改变传播方向。包括瑞利散射、米氏散射等。

非弹性散射：散射前后光子能量发生变化。包括拉曼散射、康普顿散射等。在热辐射分析中，通常只考虑弹性散射。

根据散射方向特性，又可分为：

各向同性散射：散射强度与散射角无关，在各个方向上均匀分布。

各向异性散射：散射强度随散射角变化，具有明显的方向偏好性。

2.1.2 散射的微观机制

从微观角度看，散射是入射电磁波与物质中电荷相互作用的结果：

电子极化：入射电磁波的电场使原子或分子中的电子云发生位移，形成振荡电偶极子。这个振荡偶极子成为次级波源，向各个方向辐射电磁波。

感应偶极矩：感应偶极矩 $\mathbf{p}$ 与入射电场 ${\mathbf{E}}_0$ 的关系为：

$$\mathbf{p}=\alpha {\mathbf{E}}_0$$

其中 $\alpha$ 是极化率，取决于物质的电子结构。

散射强度：单个粒子的散射强度与极化率的平方成正比：

$$I_s\propto |\alpha {|}^2{I}_0$$

其中 $I_0$ 是入射光强度。

2.2 散射相函数

2.2.1 相函数的定义

散射相函数 $\Phi (\Theta )$ 描述了散射能量随散射角 $\Theta$ 的分布规律，定义为：

$$\Phi (\Theta )=\frac{1}{{\sigma }_s}\frac{d{\sigma }_s}{d\Omega }$$

其中：

• $\Theta$ 是散射角（入射方向与散射方向之间的夹角）
• ${\sigma }_s$ 是总散射截面
• $d{\sigma }_s/d\Omega$ 是微分散射截面
• $d\Omega$ 是立体角元

相函数满足归一化条件：

$$\frac{1}{4\pi }{\int }_{4\pi }\Phi (\Theta )d\Omega =1$$

在球坐标系中，$d\Omega =\sin \Theta d\Theta d\varphi$，因此：

$$\frac{1}{2}{\int }_0^{\pi }\Phi (\Theta )\sin \Theta d\Theta =1$$

2.2.2 各向同性散射相函数

对于各向同性散射，相函数为常数：

$${\Phi }_{iso}(\Theta )=1$$

这表示散射能量在所有方向上均匀分布。

2.2.3 各向异性散射的特征

各向异性散射相函数具有以下特征：

前向散射峰：当 $\Theta \approx 0$ 时，$\Phi (\Theta )$ 出现峰值，表示大部分散射能量集中在入射方向附近。

后向散射峰：当 $\Theta \approx \pi$ 时，$\Phi (\Theta )$ 可能出现次峰值，表示部分能量被反向散射。

侧向散射：在 $\Theta =\pi /2$ 附近的散射特性。

2.3 不对称因子

2.3.1 不对称因子的定义

不对称因子 $g$（也称为各向异性因子或平均余弦）是描述散射前向性或后向性的重要参数，定义为散射角余弦的平均值：

$$g=⟨\cos \Theta ⟩=\frac{1}{4\pi }{\int }_{4\pi }\Phi (\Theta )\cos \Theta d\Omega$$

在球坐标系中：

$$g=\frac{1}{2}{\int }_0^{\pi }\Phi (\Theta )\cos \Theta \sin \Theta d\Theta$$

2.3.2 不对称因子的物理意义

不对称因子的取值范围和物理意义：

• $g=0$：各向同性散射，前向和后向散射对称
• $g>0$：前向散射占优，$g$ 越接近1，前向散射越强
• $g<0$：后向散射占优，$g$ 越接近-1，后向散射越强

典型值：

• 瑞利散射：$g=0$（各向同性）
• 大气气溶胶：$g\approx 0.6-0.8$（强前向散射）
• 生物组织：$g\approx 0.8-0.95$（极强前向散射）
• 某些人工材料：$g<0$（后向散射增强）

2.3.3 传输平均自由程

不对称因子与辐射传输密切相关。定义传输平均自由程 $l_{tr}$：

$$l_{tr}=\frac{l_s}{1-g}$$

其中 $l_s=1/{\sigma }_s$ 是散射平均自由程。

当 $g\to 1$（强前向散射）时，$l_{tr}\gg l_s$，表示光子虽然频繁散射，但传播方向改变很小，能够有效传输较长距离。

当 $g=0$（各向同性散射）时，$l_{tr}=l_s$。

当 $g<0$（后向散射）时，$l_{tr}<l_s$，表示光子传播受到更强的阻碍。

---

3. 常用散射相函数模型

3.1 Henyey-Greenstein相函数

3.1.1 数学表达式

Henyey-Greenstein（HG）相函数是描述各向异性散射最常用的经验模型，由Henyey和Greenstein于1941年提出，用于描述星际尘埃的散射：

$${\Phi }_{HG}(\Theta )=\frac{1-g^2}{(1+g^2-2g\cos \Theta {)}^{3/2}}$$

其中 $g$ 是不对称因子，取值范围为 $-1<g<1$。

3.1.2 特性分析

HG相函数具有以下重要特性：

归一化性：
$$\frac{1}{2}{\int }_0^{\pi }{\Phi }_{HG}(\Theta )\sin \Theta d\Theta =1$$

不对称因子：
$$\frac{1}{2}{\int }_0^{\pi }{\Phi }_{HG}(\Theta )\cos \Theta \sin \Theta d\Theta =g$$

极限情况：

• 当 $g\to 0$ 时，${\Phi }_{HG}(\Theta )\to 1$（各向同性散射）
• 当 $g\to 1$ 时，${\Phi }_{HG}(\Theta )$ 在 $\Theta =0$ 处出现尖锐峰值（纯前向散射）
• 当 $g\to -1$ 时，${\Phi }_{HG}(\Theta )$ 在 $\Theta =\pi$ 处出现尖锐峰值（纯后向散射）

3.1.3 优缺点

优点：

• 数学形式简单，只有一个参数 $g$
• 能够描述从各向同性到强各向异性的各种散射情况
• 计算效率高，便于数值实现
• 与许多实际系统的散射特性吻合较好

缺点：

• 不能同时描述前向和后向散射峰（对于某些颗粒系统，米氏散射理论预测同时存在前向和后向峰值）
• 在 $g$ 接近 $\pm 1$ 时，前向或后向峰值过于尖锐，与实际情况有偏差
• 缺乏明确的物理基础，是经验公式

3.1.4 应用实例

HG相函数广泛应用于：

• 大气辐射传输模型（如MODTRAN、SBDART）
• 生物医学光学（蒙特卡洛光传输模拟）
• 海洋光学
• 天体物理
• 计算机图形学

3.2 双项Henyey-Greenstein相函数

3.2.1 数学表达式

为了同时描述前向和后向散射峰，可以使用双项HG相函数：

$${\Phi }_{DHG}(\Theta )=\alpha {\Phi }_{HG}(\Theta ,g_1)+(1-\alpha ){\Phi }_{HG}(\Theta ,g_2)$$

其中：

• $\alpha$ 是权重系数，$0\le \alpha \le 1$
• $g_1$ 通常是正值（描述前向散射）
• $g_2$ 通常是负值（描述后向散射）

3.2.2 参数确定

双项HG相函数的有效不对称因子为：

$$g_{eff}=\alpha g_1+(1-\alpha )g_2$$

通过调整 $\alpha$、$g_1$、$g_2$ 三个参数，可以更灵活地拟合实际散射数据。

3.3 Legendre多项式展开

3.3.1 展开公式

相函数可以用Legendre多项式展开：

$$\Phi (\Theta )=\sum _{l=0}^{\infty }(2l+1)a_l{P}_l(\cos \Theta )$$

其中：

• $P_l(\cos \Theta )$ 是第 $l$ 阶Legendre多项式
• $a_l$ 是展开系数

展开系数由下式计算：

$$a_l=\frac{1}{2}{\int }_0^{\pi }\Phi (\Theta )P_l(\cos \Theta )\sin \Theta d\Theta$$

3.3.2 截断近似

在实际计算中，Legendre级数需要截断：

$$\Phi (\Theta )\approx \sum _{l=0}^L(2l+1)a_l{P}_l(\cos \Theta )$$

截断阶数 $L$ 的选择取决于相函数的复杂程度：

• 各向同性散射：$L=0$
• HG相函数：通常需要 $L\ge 8$ 才能较好近似
• 复杂相函数（同时具有前向和后向峰）：可能需要 $L\ge 20$

3.3.3 递推关系

Legendre多项式满足递推关系：

$$(l+1)P_{l+1}(\mu )=(2l+1)\mu P_l(\mu )-l{P}_{l-1}(\mu )$$

其中 $\mu =\cos \Theta$，初始条件为 $P_0(\mu )=1$，$P_1(\mu )=\mu$。

前几项Legendre多项式：

• $P_0(\mu )=1$
• $P_1(\mu )=\mu$
• $P_2(\mu )=\frac{1}{2}(3{\mu }^2-1)$
• $P_3(\mu )=\frac{1}{2}(5{\mu }^3-3\mu )$
• $P_4(\mu )=\frac{1}{8}(35{\mu }^4-30{\mu }^2+3)$

3.4 米氏散射理论

3.4.1 理论基础

米氏散射理论（Mie theory）是描述均匀球形颗粒电磁散射的精确理论，由Gustav Mie于1908年提出。该理论基于麦克斯韦方程组，通过求解球坐标系中的波动方程得到散射场分布。

3.4.2 尺寸参数

米氏散射由尺寸参数 $x$ 表征：

$$x=\frac{2\pi r}{\lambda }$$

其中：

• $r$ 是颗粒半径
• $\lambda$ 是入射光波长

根据尺寸参数的不同：

• $x\ll 1$：瑞利散射区
• $x\approx 1$：米氏散射区
• $x\gg 1$：几何光学区

3.4.3 相函数计算

米氏散射的相函数可以通过散射振幅函数 $S_1(\Theta )$ 和 $S_2(\Theta )$ 计算：

$$\Phi (\Theta )=\frac{|S_1(\Theta ){|}^2+|S_2(\Theta ){|}^2}{2}$$

散射振幅函数由无穷级数表示：

$$S_1(\Theta )=\sum _{n=1}^{\infty }\frac{2n+1}{n(n+1)}[a_n{\pi }_n(\cos \Theta )+b_n{\tau }_n(\cos \Theta )]$$

$$S_2(\Theta )=\sum _{n=1}^{\infty }\frac{2n+1}{n(n+1)}[a_n{\tau }_n(\cos \Theta )+b_n{\pi }_n(\cos \Theta )]$$

其中 $a_n$ 和 $b_n$ 是米氏系数，${\pi }_n$ 和 ${\tau }_n$ 是角函数。

3.4.4 米氏散射的特征

米氏散射相函数的典型特征：

前向散射峰：无论颗粒大小，米氏散射总是表现出强烈的前向散射峰。

后向散射峰：对于某些折射率和尺寸参数，会出现后向散射增强。

振荡结构：相函数随散射角呈现复杂的振荡结构，这是干涉效应的结果。

彩虹现象：在特定角度（约138°）出现强度峰值，这是光学中的彩虹现象。

3.5 其他相函数模型

3.5.1 Rayleigh相函数

对于尺寸远小于波长的颗粒（瑞利散射）：

$${\Phi }_{Rayleigh}(\Theta )=\frac{3}{4}(1+{\cos }^2\Theta )$$

这是对称的，前向和后向散射相等，不对称因子 $g=0$。

3.5.2 线性各向异性相函数

最简单的各向异性模型：

$${\Phi }_{linear}(\Theta )=1+g\cos \Theta$$

要求 $|g|\le 1$ 以保证相函数非负。

3.5.3 Delta-Eddington相函数

用于强前向散射的近似：

$${\Phi }_{DE}(\Theta )=2f\delta (1-\cos \Theta )+(1-f)(1+3g^{\prime }\cos \Theta )$$

其中 $f$ 是前向散射峰值的比例，$g^{\prime }$ 是调整后的不对称因子。

---

4. 各向异性散射的辐射传递方程

4.1 辐射传递方程的一般形式

考虑各向异性散射的稳态辐射传递方程为：

$$\frac{dI(\mathbf{r},\hat{\mathbf{s}})}{ds}=-\beta I(\mathbf{r},\hat{\mathbf{s}})+\kappa I_b(\mathbf{r})+\frac{{\sigma }_s}{4\pi }{\int }_{4\pi }\Phi ({\hat{\mathbf{s}}}^{\prime },\hat{\mathbf{s}})I(\mathbf{r},{\hat{\mathbf{s}}}^{\prime })d{\Omega }^{\prime }$$

其中：

• $I(\mathbf{r},\hat{\mathbf{s}})$ 是辐射强度
• $\beta =\kappa +{\sigma }_s$ 是消光系数
• $\kappa$ 是吸收系数
• ${\sigma }_s$ 是散射系数
• $\Phi ({\hat{\mathbf{s}}}^{\prime },\hat{\mathbf{s}})$ 是散射相函数
• 积分项表示来自所有方向的散射贡献

4.2 相函数的坐标变换

在实际计算中，相函数通常表示为入射方向 ${\hat{\mathbf{s}}}^{\prime }$ 和散射方向 $\hat{\mathbf{s}}$ 之间夹角的函数。设两个方向的球坐标分别为 $({\theta }^{\prime },{\varphi }^{\prime })$ 和 $(\theta ,\varphi )$，则散射角的余弦为：

$$\cos \Theta ={\hat{\mathbf{s}}}^{\prime }\cdot \hat{\mathbf{s}}=\cos {\theta }^{\prime }\cos \theta +\sin {\theta }^{\prime }\sin \theta \cos ({\varphi }^{\prime }-\varphi )$$

因此，相函数可以写成：

$$\Phi ({\hat{\mathbf{s}}}^{\prime },\hat{\mathbf{s}})=\Phi (\cos \Theta )=\Phi (\cos {\theta }^{\prime }\cos \theta +\sin {\theta }^{\prime }\sin \theta \cos ({\varphi }^{\prime }-\varphi ))$$

4.3 离散坐标法处理

4.3.1 角度离散

采用离散坐标法时，将立体角离散为有限个方向。设离散方向为 ${\hat{\mathbf{s}}}_j$（$j=1,2,...,M$），对应的权重为 $w_j$，满足：

$$\sum _{j=1}^M{w}_j=4\pi$$

辐射传递方程在每个离散方向上变为：

$$\frac{d{I}_j}{ds}=-\beta I_j+\kappa I_b+\frac{{\sigma }_s}{4\pi }\sum _{j^{\prime }=1}^M{w}_{j^{\prime }}{\Phi }_{j^{\prime }j}{I}_{j^{\prime }}$$

其中 ${\Phi }_{j^{\prime }j}=\Phi ({\hat{\mathbf{s}}}_{j^{\prime }},{\hat{\mathbf{s}}}_j)$。

4.3.2 散射项矩阵表示

定义散射矩阵 $\mathbf{\Phi }$，其元素为 ${\Phi }_{j^{\prime }j}$，则散射项可以写成矩阵形式：

$$\mathbf{S}=\frac{{\sigma }_s}{4\pi }\mathbf{\Phi }\cdot \mathbf{W}\cdot \mathbf{I}$$

其中 $\mathbf{W}$ 是权重对角矩阵，$\mathbf{I}$ 是辐射强度向量。

4.3.3 Legendre展开的应用

利用Legendre展开可以简化散射矩阵的计算：

$${\Phi }_{j^{\prime }j}=\sum _{l=0}^L(2l+1)a_l{P}_l(\cos {\Theta }_{j^{\prime }j})$$

其中 $\cos {\Theta }_{j^{\prime }j}={\hat{\mathbf{s}}}_{j^{\prime }}\cdot {\hat{\mathbf{s}}}_j$。

利用Legendre多项式的加法定理：

$$P_l(\cos {\Theta }_{j^{\prime }j})=\frac{4\pi }{2l+1}\sum _{m=-l}^l{Y}_{lm}({\hat{\mathbf{s}}}_{j^{\prime }})Y_{lm}^{\ast }({\hat{\mathbf{s}}}_j)$$

其中 $Y_{lm}$ 是球谐函数。这可以显著简化散射矩阵的计算和存储。

4.4 球谐函数法

4.4.1 辐射强度的展开

将球谐函数法（$P_N$ 近似）应用于各向异性散射问题时，将辐射强度展开为球谐函数的级数：

$$I(\mathbf{r},\hat{\mathbf{s}})=\sum _{l=0}^N\sum _{m=-l}^l{I}_{lm}(\mathbf{r})Y_{lm}(\hat{\mathbf{s}})$$

4.4.2 散射项的简化

利用球谐函数的正交性和Legendre展开，散射项可以大大简化：

$${\int }_{4\pi }\Phi ({\hat{\mathbf{s}}}^{\prime },\hat{\mathbf{s}})I(\mathbf{r},{\hat{\mathbf{s}}}^{\prime })d{\Omega }^{\prime }=4\pi \sum _{l=0}^N\sum _{m=-l}^l{a}_l{I}_{lm}(\mathbf{r})Y_{lm}(\hat{\mathbf{s}})$$

这使得各向异性散射的辐射传递方程转化为一组关于展开系数 $I_{lm}$ 的耦合偏微分方程。

4.4.3 $P_1$ 近似

在$P_1$近似下（只保留 $l=0,1$ 项），辐射强度表示为：

$$I(\mathbf{r},\hat{\mathbf{s}})\approx \frac{1}{4\pi }G(\mathbf{r})+\frac{3}{4\pi }\mathbf{q}(\mathbf{r})\cdot \hat{\mathbf{s}}$$

其中 $G$ 是入射辐射，$\mathbf{q}$ 是辐射热流。散射项简化为：

$$\frac{{\sigma }_s}{4\pi }{\int }_{4\pi }\Phi ({\hat{\mathbf{s}}}^{\prime },\hat{\mathbf{s}})I(\mathbf{r},{\hat{\mathbf{s}}}^{\prime })d{\Omega }^{\prime }=\frac{{\sigma }_s}{4\pi }G(\mathbf{r})+\frac{3g{\sigma }_s}{4\pi }\mathbf{q}(\mathbf{r})\cdot \hat{\mathbf{s}}$$

可以看到，不对称因子 $g$ 直接出现在辐射热流的系数中。

4.5 扩散近似

4.5.1 扩散方程推导

在光学厚介质中，可以采用扩散近似。对于各向异性散射，扩散系数 $D$ 为：

$$D=\frac{1}{3[\kappa +{\sigma }_s(1-g)]}=\frac{1}{3{\kappa }_{tr}}$$

其中 ${\kappa }_{tr}=\kappa +{\sigma }_s(1-g)$ 是传输消光系数。

4.5.2 Rosseland平均

在考虑波谱依赖时，使用Rosseland平均消光系数：

$$\frac{1}{{\kappa }_R}=\frac{{\int }_0^{\infty }\frac{1}{{\kappa }_{\lambda }+{\sigma }_{s,\lambda }(1-g_{\lambda })}\frac{d{I}_{b\lambda }}{dT}d\lambda }{{\int }_0^{\infty }\frac{d{I}_{b\lambda }}{dT}d\lambda }$$

---

5. 数值求解方法

5.1 蒙特卡洛方法

5.1.1 基本原理

蒙特卡洛方法通过追踪大量光子的随机行进来求解辐射传递问题。对于各向异性散射，关键步骤是：

1. 光子发射：根据源项分布随机产生光子
2. 自由程采样：根据消光系数采样光子传输距离
3. 散射方向采样：根据相函数采样散射后的方向
4. 权重更新：考虑吸收和散射的概率

5.1.2 散射方向采样

对于HG相函数，散射角可以通过解析方法采样。设 $\xi$ 是 $[0,1]$ 区间均匀分布的随机数，则：

$$\cos \Theta =\frac{1}{2g}[1+g^2-(\frac{1-g^2}{1-g+2g\xi }{)}^2]$$

当 $g=0$ 时，退化为 $\cos \Theta =2\xi -1$（各向同性散射）。

方位角 $\varphi$ 在各向同性介质中是均匀分布的：

$$\varphi =2\pi {\xi }^{\prime }$$

5.1.3 方向余弦更新

设散射前的方向余弦为 $(u,v,w)$，散射后的方向余弦 $(u^{\prime },v^{\prime },w^{\prime })$ 可以通过以下公式计算：

$$u^{\prime }=\frac{\sin \Theta }{\sqrt{1-w^2}}(uw\cos \varphi -v\sin \varphi )+u\cos \Theta$$

$$v^{\prime }=\frac{\sin \Theta }{\sqrt{1-w^2}}(vw\cos \varphi +u\sin \varphi )+v\cos \Theta$$

$$w^{\prime }=-\sin \Theta \sqrt{1-w^2}\cos \varphi +w\cos \Theta$$

当 $|w|\approx 1$ 时需要特殊处理。

5.2 有限差分法

5.2.1 空间离散

将计算域划分为网格，在每个网格上求解离散方向的辐射传递方程。对于一维平面几何：

$${\mu }_j\frac{d{I}_j}{dz}=-\beta I_j+\kappa I_b+\frac{{\sigma }_s}{4\pi }\sum _{j^{\prime }=1}^M{w}_{j^{\prime }}{\Phi }_{j^{\prime }j}{I}_{j^{\prime }}$$

5.2.2 迭代求解

采用源项迭代法：

1. 初始化辐射强度分布
2. 计算散射源项
3. 求解每个方向的辐射传递方程
4. 更新辐射强度
5. 检查收敛，若不收敛返回步骤2

5.2.3 收敛加速

对于强散射介质（光学厚、高反照率），收敛可能很慢。可以采用：

扩散综合加速（DSA）：利用扩散方程加速收敛

** krylov 子空间方法**：如GMRES加速迭代

5.3 有限体积法

5.3.1 控制体积分

在每个控制体和控制角上积分辐射传递方程：

$${\int }_{\Delta V}{\int }_{\Delta {\Omega }_j}\hat{\mathbf{s}}\cdot \nabla Id\Omega dV={\int }_{\Delta V}{\int }_{\Delta {\Omega }_j}[-\beta I+\kappa I_b+S_s]d\Omega dV$$

其中 $S_s$ 是散射源项。

5.3.2 通量计算

利用散度定理将左边转化为边界通量：

$$\sum _{faces}{A}_f{\int }_{\Delta {\Omega }_j}(\hat{\mathbf{s}}\cdot {\hat{\mathbf{n}}}_f)I_{f}d\Omega$$

其中 ${\hat{\mathbf{n}}}_f$ 是面法向，$I_f$ 是界面辐射强度。

5.4 离散传递法

5.4.1 方法原理

离散传递法（DTM）将辐射传递方程转化为沿特征线的常微分方程：

$$\frac{dI}{ds}+\beta I=\kappa I_b+\frac{{\sigma }_s}{4\pi }\sum _{j^{\prime }}{w}_{j^{\prime }}{\Phi }_{j^{\prime }j}{I}_{j^{\prime }}$$

5.4.2 求解过程

1. 选择一组离散方向（射线）
2. 沿每条射线追踪，计算与网格的交点
3. 在每个网格内求解辐射传递方程
4. 累加所有射线的贡献

---

6. 各向异性散射的影响分析

6.1 对辐射传输的影响

6.1.1 有效光学厚度

各向异性散射改变了辐射传输的有效光学厚度。定义传输光学厚度：

$${\tau }_{tr}={\int }_0^L[\kappa +{\sigma }_s(1-g)]ds$$

当 $g>0$（前向散射）时，${\tau }_{tr}<{\tau }_{ext}$（消光光学厚度），辐射能够传输更远距离。

当 $g<0$（后向散射）时，${\tau }_{tr}>{\tau }_{ext}$，辐射传输受到更强阻碍。

6.1.2 穿透深度

对于半无限大介质，辐射的穿透深度与传输平均自由程相关：

$$\delta \approx \frac{1}{\sqrt{3\kappa [\kappa +{\sigma }_s(1-g)]}}$$

强前向散射（$g\to 1$）显著增加穿透深度。

6.2 对温度分布的影响

6.2.1 辐射平衡温度

在辐射平衡条件下，能量方程为：

$$\nabla \cdot {\mathbf{q}}_r=0$$

各向异性散射通过改变辐射热流 ${\mathbf{q}}_r$ 来影响温度分布。

6.2.2 边界层效应

在边界附近，各向异性散射产生显著影响：

• 前向散射增强：边界层温度梯度减小
• 后向散射增强：边界层温度梯度增大

6.3 对辐射热流的影响

6.3.1 等效导热系数

在扩散极限下，可以定义等效辐射导热系数：

$$k_r=\frac{16\sigma T^3}{3[\kappa +{\sigma }_s(1-g)]}$$

前向散射（$g>0$）增大等效导热系数。

6.3.2 热流方向性

各向异性散射导致辐射热流呈现方向依赖性，在复杂几何中尤为重要。

---

7. 工程应用案例

7.1 大气气溶胶散射

大气气溶胶（尘埃、烟雾、海盐等）表现出强烈的各向异性散射特性，不对称因子通常在0.6-0.8之间。这种强前向散射特性对卫星遥感、太阳辐射传输和大气能见度有重要影响。

7.2 生物组织光学

人体组织中的细胞、胶原蛋白纤维等结构产生极强的前向散射，$g$ 值可达0.9以上。这种特性使得光能够在组织中传播数毫米，是光学相干断层扫描（OCT）和光动力治疗的基础。

7.3 燃烧系统

煤粉燃烧、喷雾燃烧等过程中，燃料颗粒的散射特性影响炉膛内的辐射传热分布。颗粒尺寸分布和折射率决定了散射的各向异性程度。

7.4 海洋光学

海水中的浮游生物、悬浮颗粒产生各向异性散射，影响海洋的光学特性和初级生产力分布。

---

8. 仿真案例分析

8.1 案例1：Henyey-Greenstein相函数特性分析

本案例通过Python仿真分析HG相函数随不对称因子的变化规律，可视化不同 $g$ 值下的散射特性。

8.2 案例2：大气辐射传输中的各向异性散射

模拟太阳光穿过大气层的过程，比较各向同性和各向异性散射对地表辐射通量的影响。

8.3 案例3：生物组织中的光传输

使用蒙特卡洛方法模拟光子在生物组织中的传输，分析前向散射对光穿透深度的影响。

8.4 案例4：工业颗粒系统的散射相函数反演

根据散射强度测量数据，利用优化算法反演颗粒系统的等效相函数参数。

---

9. 发展趋势与展望

9.1 多尺度建模

结合分子动力学、米氏理论和辐射传递方程，建立从微观到宏观的多尺度散射模型。

9.2 机器学习应用

利用神经网络学习复杂相函数的映射关系，加速辐射传递计算。

9.3 实时计算

开发GPU加速算法，实现复杂各向异性散射系统的实时辐射模拟。

9.4 实验验证

结合先进的光学测量技术（如动态光散射、角分辨散射仪），验证和改进理论模型。

---

10. 总结

各向异性散射是辐射传递中的核心物理现象，对大气科学、生物医学、工业过程等众多领域具有重要影响。本主题系统介绍了各向异性散射的理论基础、数学描述和数值方法：

1. 理论基础：从散射的微观机制出发，建立了相函数和不对称因子的概念框架

2. 相函数模型：详细介绍了HG相函数、Legendre展开、米氏理论等常用模型

3. 辐射传递方程：推导了考虑各向异性散射的辐射传递方程及其简化形式

4. 数值方法：讲解了蒙特卡洛、有限差分、有限体积等求解方法

5. 影响分析：定量分析了各向异性散射对辐射传输、温度分布和热流的影响

6. 仿真案例：通过四个完整的Python案例，展示了各向异性散射辐射问题的求解过程

掌握各向异性散射的物理机制和数值处理方法，对于准确预测参与性介质中的辐射传递行为、优化辐射换热系统设计具有重要意义。

---

参考文献

1. Chandrasekhar, S. (1960). Radiative Transfer. Dover Publications.

2. Modest, M. F. (2013). Radiative Heat Transfer (3rd ed.). Academic Press.

3. Henyey, L. G., & Greenstein, J. L. (1941). Diffuse radiation in the Galaxy. Astrophysical Journal, 93, 70-83.

4. Bohren, C. F., & Huffman, D. R. (1983). Absorption and Scattering of Light by Small Particles. Wiley.

5. Ishimaru, A. (1978). Wave Propagation and Scattering in Random Media. Academic Press.

6. Welch, A. J., & van Gemert, M. J. C. (2011). Optical-Thermal Response of Laser-Irradiated Tissue (2nd ed.). Springer.

7. Liou, K. N. (2002). An Introduction to Atmospheric Radiation (2nd ed.). Academic Press.

8. Marshak, A., & Davis, A. B. (2005). 3D Radiative Transfer in Cloudy Atmospheres. Springer.

---

附录：Python仿真代码说明

本主题配套的Python仿真程序包含以下四个案例：

案例1：HG相函数特性分析

• 可视化不同不对称因子下的相函数曲线
• 计算和比较Legendre展开系数
• 分析前向/后向散射比例

案例2：大气辐射传输

• 建立一维大气辐射传输模型
• 比较各向同性和各向异性散射
• 计算地表辐射通量

案例3：生物组织光传输

• 蒙特卡洛光子追踪
• 分析穿透深度与不对称因子的关系
• 可视化光子传播路径

案例4：相函数反演

• 基于模拟测量数据的参数反演
• 使用最小二乘优化
• 验证反演精度

运行命令：

python run_simulation.py

生成的文件：

• case1_phase_function.png - 相函数特性分析图
• case1_phase_function.gif - 相函数演化动画
• case2_atmosphere.png - 大气辐射传输结果
• case2_atmosphere.gif - 辐射传输动画
• case3_tissue.png - 生物组织光传输结果
• case3_tissue.gif - 光子传播动画
• case4_inversion.png - 相函数反演结果
• case4_inversion.gif - 反演过程动画

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
主题063：各向异性散射介质辐射仿真程序
==========================================
本程序包含四个仿真案例：
1. Henyey-Greenstein相函数特性分析
2. 大气辐射传输中的各向异性散射
3. 生物组织中的光传输
4. 散射相函数反演

作者: 仿真专家
日期: 2026-03-07
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch
import matplotlib.animation as animation
from scipy import integrate, optimize, special
from scipy.special import legendre, eval_legendre
import warnings
warnings.filterwarnings('ignore')

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

# ============================================================================
# 案例1：Henyey-Greenstein相函数特性分析
# ============================================================================

def case1_phase_function():
    """
    案例1：Henyey-Greenstein相函数特性分析
    分析不同不对称因子g下的相函数特性
    """
    print("=" * 70)
    print("案例1：Henyey-Greenstein相函数特性分析")
    print("=" * 70)
    
    # 散射角
    theta = np.linspace(0, np.pi, 180)
    mu = np.cos(theta)
    
    # 不同不对称因子
    g_values = [-0.8, -0.5, -0.3, 0.0, 0.3, 0.5, 0.8, 0.95]
    colors = plt.cm.RdYlBu(np.linspace(0, 1, len(g_values)))
    
    # Henyey-Greenstein相函数
    def hg_phase_function(mu, g):
        """Henyey-Greenstein相函数"""
        if abs(g) < 1e-10:
            return np.ones_like(mu)
        return (1 - g**2) / ((1 + g**2 - 2*g*mu)**1.5)
    
    # 计算各g值下的相函数
    phase_functions = {}
    for g in g_values:
        phase_functions[g] = hg_phase_function(mu, g)
    
    # 计算Legendre展开系数
    def legendre_coefficients(g, n_terms=10):
        """计算HG相函数的Legendre展开系数"""
        coeffs = []
        for l in range(n_terms):
            # a_l = g^l 对于HG相函数
            coeffs.append(g**l)
        return coeffs
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 1. 极坐标下的相函数
    ax1 = fig.add_subplot(2, 3, 1, projection='polar')
    for i, g in enumerate(g_values):
        # 归一化
        pf = phase_functions[g]
        pf_norm = pf / np.max(pf)
        ax1.plot(theta, pf_norm, color=colors[i], linewidth=2, label=f'g={g}')
    ax1.set_theta_zero_location('N')
    ax1.set_theta_direction(-1)
    ax1.set_title('极坐标相函数分布', fontsize=12, fontweight='bold', pad=20)
    ax1.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1))
    
    # 2. 直角坐标下的相函数
    ax2 = fig.add_subplot(2, 3, 2)
    for i, g in enumerate(g_values):
        ax2.semilogy(theta * 180/np.pi, phase_functions[g], 
                     color=colors[i], linewidth=2, label=f'g={g}')
    ax2.set_xlabel('散射角 (度)', fontsize=11)
    ax2.set_ylabel('相函数值', fontsize=11)
    ax2.set_title('Henyey-Greenstein相函数', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.axvline(x=90, color='k', linestyle='--', alpha=0.5)
    
    # 3. Legendre展开系数
    ax3 = fig.add_subplot(2, 3, 3)
    l_terms = np.arange(10)
    for i, g in enumerate([0.0, 0.3, 0.5, 0.8, 0.95]):
        coeffs = legendre_coefficients(g, 10)
        ax3.semilogy(l_terms, np.abs(coeffs), 'o-', 
                     color=colors[i+3], linewidth=2, label=f'g={g}')
    ax3.set_xlabel('Legendre阶数 l', fontsize=11)
    ax3.set_ylabel('|a_l|', fontsize=11)
    ax3.set_title('Legendre展开系数', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 前向/后向散射比例
    ax4 = fig.add_subplot(2, 3, 4)
    forward_ratio = []
    backward_ratio = []
    g_list = np.linspace(-0.95, 0.95, 20)
    
    for g in g_list:
        pf = hg_phase_function(mu, g)
        # 归一化
        pf = pf / (0.5 * integrate.simpson(pf * np.sin(theta), theta))
        # 前向散射 (0-90度)
        forward_idx = theta <= np.pi/2
        forward = 0.5 * integrate.simpson(pf[forward_idx] * np.sin(theta[forward_idx]), 
                                          theta[forward_idx])
        # 后向散射 (90-180度)
        backward_idx = theta > np.pi/2
        backward = 0.5 * integrate.simpson(pf[backward_idx] * np.sin(theta[backward_idx]), 
                                           theta[backward_idx])
        forward_ratio.append(forward)
        backward_ratio.append(backward)
    
    ax4.plot(g_list, forward_ratio, 'b-', linewidth=2, label='前向散射')
    ax4.plot(g_list, backward_ratio, 'r-', linewidth=2, label='后向散射')
    ax4.axvline(x=0, color='k', linestyle='--', alpha=0.5)
    ax4.axhline(y=0.5, color='k', linestyle='--', alpha=0.5)
    ax4.set_xlabel('不对称因子 g', fontsize=11)
    ax4.set_ylabel('散射能量比例', fontsize=11)
    ax4.set_title('前向/后向散射比例', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. 传输平均自由程
    ax5 = fig.add_subplot(2, 3, 5)
    sigma_s = 1.0  # 散射系数
    kappa = 0.1    # 吸收系数
    g_range = np.linspace(-0.99, 0.99, 100)
    l_s = 1.0 / sigma_s  # 散射平均自由程
    l_tr = l_s / (1 - g_range)  # 传输平均自由程
    
    ax5.semilogy(g_range, l_tr, 'purple', linewidth=2)
    ax5.axvline(x=0, color='k', linestyle='--', alpha=0.5)
    ax5.axhline(y=l_s, color='r', linestyle='--', alpha=0.5, label=f'l_s = {l_s:.2f}')
    ax5.set_xlabel('不对称因子 g', fontsize=11)
    ax5.set_ylabel('传输平均自由程 l_tr', fontsize=11)
    ax5.set_title('传输平均自由程 vs g', fontsize=12, fontweight='bold')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. 相函数3D可视化示意
    ax6 = fig.add_subplot(2, 3, 6, projection='3d')
    theta_3d = np.linspace(0, np.pi, 50)
    phi_3d = np.linspace(0, 2*np.pi, 50)
    THETA, PHI = np.meshgrid(theta_3d, phi_3d)
    
    g_3d = 0.8
    R = hg_phase_function(np.cos(THETA), g_3d)
    R = R / np.max(R) * 0.8
    
    X = R * np.sin(THETA) * np.cos(PHI)
    Y = R * np.sin(THETA) * np.sin(PHI)
    Z = R * np.cos(THETA)
    
    ax6.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
    ax6.set_title(f'3D相函数 (g={g_3d})', fontsize=12, fontweight='bold')
    ax6.set_xlabel('X')
    ax6.set_ylabel('Y')
    ax6.set_zlabel('Z')
    
    plt.tight_layout()
    plt.savefig('case1_phase_function.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("  结果已保存至: case1_phase_function.png")
    
    # 创建动画：相函数随g值演化
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax1 = axes[0]
    line_polar, = ax1.plot([], [], 'b-', linewidth=2)
    ax1.set_xlim(0, np.pi)
    ax1.set_ylim(0, 10)
    ax1.set_xlabel('散射角 (rad)', fontsize=11)
    ax1.set_ylabel('相函数值', fontsize=11)
    ax1.set_title('HG相函数演化', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[1]
    # 极坐标
    ax2_polar = fig.add_axes([0.6, 0.15, 0.35, 0.7], projection='polar')
    line_radial, = ax2_polar.plot([], [], 'r-', linewidth=2)
    ax2_polar.set_theta_zero_location('N')
    ax2_polar.set_theta_direction(-1)
    ax2_polar.set_title('极坐标分布', fontsize=12, fontweight='bold', pad=20)
    ax2_polar.set_ylim(0, 5)
    
    # 隐藏原来的ax2
    ax2.axis('off')
    
    g_anim_values = np.linspace(-0.9, 0.9, 100)
    
    def animate_phase(frame):
        g = g_anim_values[frame]
        pf = hg_phase_function(mu, g)
        
        line_polar.set_data(theta, pf)
        
        # 更新极坐标
        pf_norm = pf / np.max(pf) * 3
        line_radial.set_data(theta, pf_norm)
        
        ax1.set_title(f'HG相函数 (g={g:.2f})', fontsize=12, fontweight='bold')
        
        return line_polar, line_radial
    
    anim = animation.FuncAnimation(fig, animate_phase, frames=100, interval=100, blit=True)
    anim.save('case1_phase_function.gif', writer='pillow', fps=10, dpi=100)
    plt.close()
    
    print("  动画已保存至: case1_phase_function.gif")
    
    return phase_functions, g_values


# ============================================================================
# 案例2：大气辐射传输中的各向异性散射
# ============================================================================

def case2_atmosphere():
    """
    案例2：大气辐射传输中的各向异性散射
    比较各向同性和各向异性散射对太阳辐射传输的影响
    """
    print("\n" + "=" * 70)
    print("案例2：大气辐射传输中的各向异性散射")
    print("=" * 70)
    
    # 大气参数
    tau_total = 5.0  # 总光学厚度
    n_layers = 50    # 分层数
    tau = np.linspace(0, tau_total, n_layers)
    dtau = tau[1] - tau[0]
    
    # 散射反照率
    omega = 0.9
    
    # 太阳入射角度
    mu_0 = np.cos(np.deg2rad(30))  # 天顶角30度
    
    # 不对称因子
    g_iso = 0.0   # 各向同性
    g_aniso = 0.75  # 各向异性（典型气溶胶值）
    
    # 离散坐标法参数
    n_angles = 16
    mu_angles, weights = np.polynomial.legendre.leggauss(n_angles)
    # 转换到[0,1]区间
    mu_angles = 0.5 * (mu_angles + 1)
    weights = 0.5 * weights
    
    # Henyey-Greenstein相函数
    def hg_phase(mu_scatter, g):
        if abs(g) < 1e-10:
            return np.ones_like(mu_scatter)
        return (1 - g**2) / ((1 + g**2 - 2*g*mu_scatter)**1.5)
    
    # 计算散射相函数矩阵
    def compute_phase_matrix(mu_angles, g):
        n = len(mu_angles)
        P = np.zeros((n, n))
        for i in range(n):
            for j in range(n):
                # 散射角余弦
                mu_scatter = mu_angles[i] * mu_angles[j]
                P[i, j] = hg_phase(mu_scatter, g)
        # 归一化
        for i in range(n):
            P[i, :] = P[i, :] / np.sum(P[i, :] * weights)
        return P
    
    # 求解辐射传递方程
    def solve_rte(tau, mu_angles, weights, omega, g, mu_0, max_iter=200):
        n_layers = len(tau)
        n_angles = len(mu_angles)
        
        # 辐射强度 [layer, angle]
        I = np.zeros((n_layers, n_angles))
        
        # 相函数矩阵
        P = compute_phase_matrix(mu_angles, g)
        
        # 源项
        S = np.zeros((n_layers, n_angles))
        
        # 迭代求解
        for iteration in range(max_iter):
            I_old = I.copy()
            
            # 计算源项（包含散射）
            for k in range(n_layers):
                for i in range(n_angles):
                    scatter_sum = 0
                    for j in range(n_angles):
                        scatter_sum += weights[j] * P[i, j] * I[k, j]
                    S[k, i] = omega * scatter_sum
            
            # 正向传播 (mu > 0)
            for i in range(n_angles // 2, n_angles):
                I[0, i] = 1.0  # 入射边界
                for k in range(1, n_layers):
                    I[k, i] = (I[k-1, i] * (mu_angles[i] / dtau) + S[k, i]) / \
                              (mu_angles[i] / dtau + 1)
            
            # 反向传播 (mu < 0)
            for i in range(n_angles // 2):
                I[-1, i] = 0.0  # 出射边界
                for k in range(n_layers - 2, -1, -1):
                    I[k, i] = (I[k+1, i] * (-mu_angles[i] / dtau) + S[k, i]) / \
                              (-mu_angles[i] / dtau + 1)
            
            # 检查收敛
            if np.max(np.abs(I - I_old)) < 1e-6:
                print(f"    g={g}: 收敛于第 {iteration+1} 次迭代")
                break
        
        return I
    
    # 求解各向同性散射
    print("  求解各向同性散射 (g=0)...")
    I_iso = solve_rte(tau, mu_angles, weights, omega, g_iso, mu_0)
    
    # 求解各向异性散射
    print("  求解各向异性散射 (g=0.75)...")
    I_aniso = solve_rte(tau, mu_angles, weights, omega, g_aniso, mu_0)
    
    # 计算辐射通量
    flux_iso = np.zeros(n_layers)
    flux_aniso = np.zeros(n_layers)
    for k in range(n_layers):
        flux_iso[k] = 2 * np.pi * np.sum(weights * mu_angles * I_iso[k, :])
        flux_aniso[k] = 2 * np.pi * np.sum(weights * mu_angles * I_aniso[k, :])
    
    # 创建可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 1. 辐射强度分布对比
    ax = axes[0, 0]
    for i in [0, n_angles//4, n_angles//2, 3*n_angles//4]:
        ax.plot(I_iso[:, i], tau, '--', linewidth=2, label=f'μ={mu_angles[i]:.2f} (iso)')
        ax.plot(I_aniso[:, i], tau, '-', linewidth=2, label=f'μ={mu_angles[i]:.2f} (aniso)')
    ax.invert_yaxis()
    ax.set_xlabel('辐射强度', fontsize=11)
    ax.set_ylabel('光学厚度 τ', fontsize=11)
    ax.set_title('辐射强度分布对比', fontsize=12, fontweight='bold')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)
    
    # 2. 辐射通量对比
    ax = axes[0, 1]
    ax.plot(flux_iso, tau, 'b--', linewidth=2, label='各向同性 (g=0)')
    ax.plot(flux_aniso, tau, 'r-', linewidth=2, label='各向异性 (g=0.75)')
    ax.invert_yaxis()
    ax.set_xlabel('辐射通量', fontsize=11)
    ax.set_ylabel('光学厚度 τ', fontsize=11)
    ax.set_title('辐射通量对比', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 透射率对比
    ax = axes[0, 2]
    T_iso = I_iso[-1, :] / I_iso[0, :]
    T_aniso = I_aniso[-1, :] / I_aniso[0, :]
    ax.plot(mu_angles, T_iso, 'bo--', linewidth=2, markersize=6, label='各向同性')
    ax.plot(mu_angles, T_aniso, 'rs-', linewidth=2, markersize=6, label='各向异性')
    ax.set_xlabel('角度余弦 μ', fontsize=11)
    ax.set_ylabel('透射率', fontsize=11)
    ax.set_title('角度分辨透射率', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 4. 不同g值的透射率
    ax = axes[1, 0]
    g_values = [0.0, 0.3, 0.5, 0.7, 0.9]
    colors = plt.cm.viridis(np.linspace(0, 1, len(g_values)))
    for i, g in enumerate(g_values):
        print(f"  计算 g={g}...")
        I_g = solve_rte(tau, mu_angles, weights, omega, g, mu_0, max_iter=100)
        flux_g = np.array([2 * np.pi * np.sum(weights * mu_angles * I_g[k, :]) 
                          for k in range(n_layers)])
        ax.plot(flux_g, tau, color=colors[i], linewidth=2, label=f'g={g}')
    ax.invert_yaxis()
    ax.set_xlabel('辐射通量', fontsize=11)
    ax.set_ylabel('光学厚度 τ', fontsize=11)
    ax.set_title('不同g值的辐射通量', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 5. 有效穿透深度
    ax = axes[1, 1]
    g_range = np.linspace(0, 0.95, 20)
    penetration_depths = []
    for g in g_range:
        I_g = solve_rte(tau, mu_angles, weights, omega, g, mu_0, max_iter=50)
        flux_g = np.array([2 * np.pi * np.sum(weights * mu_angles * I_g[k, :]) 
                          for k in range(n_layers)])
        # 定义穿透深度为通量衰减到1/e的位置
        if flux_g[0] > 0:
            target_flux = flux_g[0] / np.e
            idx = np.argmin(np.abs(flux_g - target_flux))
            penetration_depths.append(tau[idx])
        else:
            penetration_depths.append(0)
    
    ax.plot(g_range, penetration_depths, 'purple', linewidth=2, marker='o')
    ax.set_xlabel('不对称因子 g', fontsize=11)
    ax.set_ylabel('有效穿透深度 (τ)', fontsize=11)
    ax.set_title('穿透深度 vs 不对称因子', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 6. 能量守恒检查
    ax = axes[1, 2]
    # 计算总能量
    total_energy_iso = np.sum(weights * I_iso[0, :])
    total_energy_aniso = np.sum(weights * I_aniso[0, :])
    
    categories = ['入射', '透射', '反射']
    x_pos = np.arange(len(categories))
    
    # 各向同性
    incident_iso = np.sum(weights * I_iso[0, :])
    transmitted_iso = np.sum(weights * I_iso[-1, :])
    reflected_iso = incident_iso - transmitted_iso
    values_iso = [incident_iso, transmitted_iso, reflected_iso]
    
    # 各向异性
    incident_aniso = np.sum(weights * I_aniso[0, :])
    transmitted_aniso = np.sum(weights * I_aniso[-1, :])
    reflected_aniso = incident_aniso - transmitted_aniso
    values_aniso = [incident_aniso, transmitted_aniso, reflected_aniso]
    
    width = 0.35
    ax.bar(x_pos - width/2, values_iso, width, label='各向同性', alpha=0.8)
    ax.bar(x_pos + width/2, values_aniso, width, label='各向异性', alpha=0.8)
    ax.set_ylabel('能量', fontsize=11)
    ax.set_title('能量守恒检查', fontsize=12, fontweight='bold')
    ax.set_xticks(x_pos)
    ax.set_xticklabels(categories)
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('case2_atmosphere.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("  结果已保存至: case2_atmosphere.png")
    
    # 创建动画：辐射传输过程
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax1 = axes[0]
    line_flux_iso, = ax1.plot([], [], 'b--', linewidth=2, label='各向同性')
    line_flux_aniso, = ax1.plot([], [], 'r-', linewidth=2, label='各向异性')
    ax1.set_xlim(0, 2)
    ax1.set_ylim(0, tau_total)
    ax1.invert_yaxis()
    ax1.set_xlabel('辐射通量', fontsize=11)
    ax1.set_ylabel('光学厚度 τ', fontsize=11)
    ax1.set_title('辐射通量演化', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[1]
    # 角度分布
    line_angle_iso, = ax2.plot([], [], 'b--', linewidth=2, label='各向同性')
    line_angle_aniso, = ax2.plot([], [], 'r-', linewidth=2, label='各向异性')
    ax2.set_xlim(0, 1)
    ax2.set_ylim(0, 2)
    ax2.set_xlabel('角度余弦 μ', fontsize=11)
    ax2.set_ylabel('辐射强度', fontsize=11)
    ax2.set_title('角度分布演化', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    def animate_atmosphere(frame):
        layer = int(frame * n_layers / 100)
        layer = min(layer, n_layers - 1)
        
        # 更新通量
        flux_iso_anim = flux_iso[:layer+1]
        flux_aniso_anim = flux_aniso[:layer+1]
        tau_anim = tau[:layer+1]
        
        line_flux_iso.set_data(flux_iso_anim, tau_anim)
        line_flux_aniso.set_data(flux_aniso_anim, tau_anim)
        
        # 更新角度分布
        line_angle_iso.set_data(mu_angles, I_iso[layer, :])
        line_angle_aniso.set_data(mu_angles, I_aniso[layer, :])
        
        return line_flux_iso, line_flux_aniso, line_angle_iso, line_angle_aniso
    
    anim = animation.FuncAnimation(fig, animate_atmosphere, frames=100, interval=100, blit=True)
    anim.save('case2_atmosphere.gif', writer='pillow', fps=10, dpi=100)
    plt.close()
    
    print("  动画已保存至: case2_atmosphere.gif")
    
    return I_iso, I_aniso, flux_iso, flux_aniso


# ============================================================================
# 案例3：生物组织中的光传输
# ============================================================================

def case3_tissue():
    """
    案例3：生物组织中的光传输
    使用蒙特卡洛方法模拟光子在生物组织中的传输
    """
    print("\n" + "=" * 70)
    print("案例3：生物组织中的光传输")
    print("=" * 70)
    
    # 生物组织光学参数
    tissue_params = {
        'skin': {'mu_a': 0.1, 'mu_s': 100.0, 'g': 0.85},
        'brain': {'mu_a': 0.2, 'mu_s': 200.0, 'g': 0.90},
        'muscle': {'mu_a': 0.15, 'mu_s': 150.0, 'g': 0.88},
        'blood': {'mu_a': 2.0, 'mu_s': 500.0, 'g': 0.95}
    }
    
    # 选择组织类型
    tissue_type = 'brain'
    params = tissue_params[tissue_type]
    mu_a = params['mu_a']  # 吸收系数 (1/cm)
    mu_s = params['mu_s']  # 散射系数 (1/cm)
    g = params['g']        # 不对称因子
    mu_t = mu_a + mu_s     # 总衰减系数
    
    # 计算传输系数
    mu_tr = mu_a + mu_s * (1 - g)  # 传输系数
    
    print(f"  组织类型: {tissue_type}")
    print(f"  吸收系数 μa = {mu_a} cm^-1")
    print(f"  散射系数 μs = {mu_s} cm^-1")
    print(f"  不对称因子 g = {g}")
    print(f"  传输系数 μtr = {mu_tr:.2f} cm^-1")
    
    # 蒙特卡洛模拟参数
    n_photons = 5000  # 光子数
    max_steps = 1000  # 最大步数
    
    # 初始化
    photon_positions = []
    photon_weights = []
    photon_status = []  # 0:存活, 1:吸收, 2:逃逸
    
    print(f"  模拟 {n_photons} 个光子...")
    
    for p in range(n_photons):
        # 光子初始位置和方向
        x, y, z = 0.0, 0.0, 0.0
        ux, uy, uz = 0.0, 0.0, 1.0  # 初始向上传播
        weight = 1.0
        
        path_x = [x]
        path_y = [y]
        path_z = [z]
        status = 0
        
        for step in range(max_steps):
            # 采样自由程
            s = -np.log(np.random.random()) / mu_t
            
            # 更新位置
            x += s * ux
            y += s * uy
            z += s * uz
            
            path_x.append(x)
            path_y.append(y)
            path_z.append(z)
            
            # 检查边界 (模拟1cm厚的组织)
            if z > 1.0 or z < 0:
                status = 2  # 逃逸
                break
            
            # 吸收或散射
            if np.random.random() < mu_a / mu_t:
                # 吸收
                weight *= 0.9
                if weight < 0.01:
                    status = 1
                    break
            else:
                # 散射 - HG相函数采样
                if abs(g) < 1e-6:
                    cos_theta = 2 * np.random.random() - 1
                else:
                    # HG相函数采样
                    xi = np.random.random()
                    cos_theta = (1 + g**2 - ((1 - g**2) / (1 - g + 2*g*xi))**2) / (2*g)
                
                sin_theta = np.sqrt(1 - cos_theta**2)
                phi = 2 * np.pi * np.random.random()
                
                # 更新方向
                if abs(uz) > 0.99999:
                    ux_new = sin_theta * np.cos(phi)
                    uy_new = sin_theta * np.sin(phi)
                    uz_new = cos_theta * np.sign(uz)
                else:
                    ux_new = sin_theta * (ux * uz * np.cos(phi) - uy * np.sin(phi)) / np.sqrt(1 - uz**2) + ux * cos_theta
                    uy_new = sin_theta * (uy * uz * np.cos(phi) + ux * np.sin(phi)) / np.sqrt(1 - uz**2) + uy * cos_theta
                    uz_new = -sin_theta * np.cos(phi) * np.sqrt(1 - uz**2) + uz * cos_theta
                
                # 归一化
                norm = np.sqrt(ux_new**2 + uy_new**2 + uz_new**2)
                ux, uy, uz = ux_new/norm, uy_new/norm, uz_new/norm
        
        photon_positions.append((np.array(path_x), np.array(path_y), np.array(path_z)))
        photon_weights.append(weight)
        photon_status.append(status)
    
    # 统计结果
    n_absorbed = np.sum(np.array(photon_status) == 1)
    n_escaped = np.sum(np.array(photon_status) == 2)
    n_alive = np.sum(np.array(photon_status) == 0)
    
    print(f"  吸收: {n_absorbed} ({n_absorbed/n_photons*100:.1f}%)")
    print(f"  逃逸: {n_escaped} ({n_escaped/n_photons*100:.1f}%)")
    print(f"  存活: {n_alive} ({n_alive/n_photons*100:.1f}%)")
    
    # 计算穿透深度统计
    max_depths = []
    for px, py, pz in photon_positions:
        max_depths.append(np.max(pz))
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 10))
    
    # 1. 光子路径3D投影
    ax1 = fig.add_subplot(2, 3, 1, projection='3d')
    n_show = min(50, n_photons)
    for i in range(n_show):
        px, py, pz = photon_positions[i]
        color = 'red' if photon_status[i] == 1 else 'blue' if photon_status[i] == 2 else 'green'
        ax1.plot(px, py, pz, color=color, linewidth=0.5, alpha=0.6)
    ax1.set_xlabel('X (cm)')
    ax1.set_ylabel('Y (cm)')
    ax1.set_zlabel('Z (cm)')
    ax1.set_title('光子路径3D视图', fontsize=12, fontweight='bold')
    
    # 2. XY平面投影
    ax2 = fig.add_subplot(2, 3, 2)
    for i in range(n_show):
        px, py, pz = photon_positions[i]
        color = 'red' if photon_status[i] == 1 else 'blue' if photon_status[i] == 2 else 'green'
        ax2.plot(px, py, color=color, linewidth=0.5, alpha=0.6)
    ax2.set_xlabel('X (cm)', fontsize=11)
    ax2.set_ylabel('Y (cm)', fontsize=11)
    ax2.set_title('XY平面投影', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    ax2.grid(True, alpha=0.3)
    
    # 3. XZ平面投影
    ax3 = fig.add_subplot(2, 3, 3)
    for i in range(n_show):
        px, py, pz = photon_positions[i]
        color = 'red' if photon_status[i] == 1 else 'blue' if photon_status[i] == 2 else 'green'
        ax3.plot(px, pz, color=color, linewidth=0.5, alpha=0.6)
    ax3.axhline(y=1.0, color='k', linestyle='--', alpha=0.5, label='组织边界')
    ax3.set_xlabel('X (cm)', fontsize=11)
    ax3.set_ylabel('Z (cm)', fontsize=11)
    ax3.set_title('XZ平面投影', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    ax3.legend()
    
    # 4. 穿透深度分布
    ax4 = fig.add_subplot(2, 3, 4)
    ax4.hist(max_depths, bins=30, color='purple', alpha=0.7, edgecolor='black')
    ax4.axvline(x=np.mean(max_depths), color='r', linestyle='--', 
                linewidth=2, label=f'均值: {np.mean(max_depths):.3f} cm')
    ax4.axvline(x=1/mu_tr, color='g', linestyle='--', 
                linewidth=2, label=f'理论: {1/mu_tr:.3f} cm')
    ax4.set_xlabel('最大穿透深度 (cm)', fontsize=11)
    ax4.set_ylabel('光子数', fontsize=11)
    ax4.set_title('穿透深度分布', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. 不同g值的穿透深度比较
    ax5 = fig.add_subplot(2, 3, 5)
    g_values = [0.0, 0.5, 0.7, 0.85, 0.9, 0.95]
    penetration_means = []
    
    for g_test in g_values:
        mu_tr_test = mu_a + mu_s * (1 - g_test)
        # 简化计算：使用指数衰减近似
        mean_depth = 1 / mu_tr_test
        penetration_means.append(mean_depth)
    
    ax5.plot(g_values, penetration_means, 'o-', linewidth=2, markersize=8, color='orange')
    ax5.set_xlabel('不对称因子 g', fontsize=11)
    ax5.set_ylabel('平均穿透深度 (cm)', fontsize=11)
    ax5.set_title('穿透深度 vs 不对称因子', fontsize=12, fontweight='bold')
    ax5.grid(True, alpha=0.3)
    
    # 6. 能量沉积分布
    ax6 = fig.add_subplot(2, 3, 6)
    z_bins = np.linspace(0, 1, 21)
    energy_deposit = np.zeros(len(z_bins)-1)
    
    for i in range(n_photons):
        px, py, pz = photon_positions[i]
        weight = photon_weights[i]
        for j in range(len(pz)-1):
            z_mid = (pz[j] + pz[j+1]) / 2
            if 0 <= z_mid < 1:
                bin_idx = int(z_mid / 0.05)
                if bin_idx < len(energy_deposit):
                    energy_deposit[bin_idx] += weight * mu_a / mu_t
    
    z_centers = (z_bins[:-1] + z_bins[1:]) / 2
    ax6.barh(z_centers, energy_deposit, height=0.04, color='brown', alpha=0.7, edgecolor='black')
    ax6.set_xlabel('能量沉积', fontsize=11)
    ax6.set_ylabel('深度 Z (cm)', fontsize=11)
    ax6.set_title('能量沉积分布', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case3_tissue.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("  结果已保存至: case3_tissue.png")
    
    # 创建动画：光子传播过程
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    ax1 = axes[0]
    lines_xy = [ax1.plot([], [], linewidth=0.8, alpha=0.6)[0] for _ in range(30)]
    ax1.set_xlim(-0.5, 0.5)
    ax1.set_ylim(-0.5, 0.5)
    ax1.set_xlabel('X (cm)', fontsize=11)
    ax1.set_ylabel('Y (cm)', fontsize=11)
    ax1.set_title('XY平面光子传播', fontsize=12, fontweight='bold')
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[1]
    lines_xz = [ax2.plot([], [], linewidth=0.8, alpha=0.6)[0] for _ in range(30)]
    ax2.set_xlim(-0.5, 0.5)
    ax2.set_ylim(0, 1)
    ax2.axhline(y=1.0, color='k', linestyle='--', alpha=0.5)
    ax2.set_xlabel('X (cm)', fontsize=11)
    ax2.set_ylabel('Z (cm)', fontsize=11)
    ax2.set_title('XZ平面光子传播', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    def animate_tissue(frame):
        progress = (frame + 1) / 100
        
        for i, (line_xy, line_xz) in enumerate(zip(lines_xy, lines_xz)):
            if i < len(photon_positions):
                px, py, pz = photon_positions[i]
                n_points = int(len(px) * progress)
                if n_points > 1:
                    line_xy.set_data(px[:n_points], py[:n_points])
                    line_xz.set_data(px[:n_points], pz[:n_points])
                    # 根据状态设置颜色
                    if photon_status[i] == 1:
                        line_xy.set_color('red')
                        line_xz.set_color('red')
                    elif photon_status[i] == 2:
                        line_xy.set_color('blue')
                        line_xz.set_color('blue')
                    else:
                        line_xy.set_color('green')
                        line_xz.set_color('green')
        
        return lines_xy + lines_xz
    
    anim = animation.FuncAnimation(fig, animate_tissue, frames=100, interval=100, blit=True)
    anim.save('case3_tissue.gif', writer='pillow', fps=10, dpi=100)
    plt.close()
    
    print("  动画已保存至: case3_tissue.gif")
    
    return photon_positions, photon_weights, photon_status


# ============================================================================
# 案例4：散射相函数反演
# ============================================================================

def case4_inversion():
    """
    案例4：散射相函数反演
    基于模拟测量数据反演HG相函数参数
    """
    print("\n" + "=" * 70)
    print("案例4：散射相函数反演")
    print("=" * 70)
    
    # 生成模拟测量数据
    theta_meas = np.linspace(0, np.pi, 36)  # 测量角度
    mu_meas = np.cos(theta_meas)
    
    # 真实参数
    g_true = 0.75
    
    # HG相函数
    def hg_phase(mu, g):
        if abs(g) < 1e-10:
            return np.ones_like(mu)
        return (1 - g**2) / ((1 + g**2 - 2*g*mu)**1.5)
    
    # 生成带噪声的测量数据
    np.random.seed(42)
    phase_true = hg_phase(mu_meas, g_true)
    noise_level = 0.05
    phase_meas = phase_true * (1 + noise_level * np.random.randn(len(theta_meas)))
    phase_meas = np.maximum(phase_meas, 0.1)  # 确保正值
    
    print(f"  真实不对称因子: g_true = {g_true}")
    print(f"  测量噪声水平: {noise_level*100:.1f}%")
    
    # 反演方法1：最小二乘法
    def residual(g, mu, phase_meas):
        phase_model = hg_phase(mu, g[0])
        return phase_meas - phase_model
    
    result_ls = optimize.least_squares(residual, x0=[0.5], args=(mu_meas, phase_meas), 
                                        bounds=([-0.99], [0.99]))
    g_inverted_ls = result_ls.x[0]
    
    print(f"  最小二乘反演结果: g = {g_inverted_ls:.4f}")
    print(f"  相对误差: {abs(g_inverted_ls - g_true)/g_true*100:.2f}%")
    
    # 反演方法2：最大似然估计
    def neg_log_likelihood(g, mu, phase_meas):
        phase_model = hg_phase(mu, g[0])
        # 假设高斯噪声
        sigma = noise_level * phase_model
        nll = 0.5 * np.sum(((phase_meas - phase_model) / sigma)**2 + np.log(2*np.pi*sigma**2))
        return nll
    
    result_mle = optimize.minimize(neg_log_likelihood, x0=[0.5], args=(mu_meas, phase_meas),
                                   bounds=[(-0.99, 0.99)])
    g_inverted_mle = result_mle.x[0]
    
    print(f"  最大似然反演结果: g = {g_inverted_mle:.4f}")
    print(f"  相对误差: {abs(g_inverted_mle - g_true)/g_true*100:.2f}%")
    
    # 反演方法3：Legendre展开系数拟合
    def legendre_fit(g, mu, phase_meas, n_terms=5):
        phase_model = hg_phase(mu, g[0])
        # 计算Legendre系数
        coeffs_model = []
        coeffs_meas = []
        for l in range(n_terms):
            P_l = eval_legendre(l, mu)
            coeff_model = np.sum(phase_model * P_l) / len(mu)
            coeff_meas = np.sum(phase_meas * P_l) / len(mu)
            coeffs_model.append(coeff_model)
            coeffs_meas.append(coeff_meas)
        return np.array(coeffs_model) - np.array(coeffs_meas)
    
    result_leg = optimize.least_squares(legendre_fit, x0=[0.5], 
                                        args=(mu_meas, phase_meas, 5),
                                        bounds=([-0.99], [0.99]))
    g_inverted_leg = result_leg.x[0]
    
    print(f"  Legendre拟合反演结果: g = {g_inverted_leg:.4f}")
    print(f"  相对误差: {abs(g_inverted_leg - g_true)/g_true*100:.2f}%")
    
    # 创建可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 1. 测量数据与拟合结果对比
    ax = axes[0, 0]
    theta_deg = theta_meas * 180 / np.pi
    ax.semilogy(theta_deg, phase_meas, 'ko', markersize=8, label='测量数据')
    ax.semilogy(theta_deg, phase_true, 'b-', linewidth=2, label=f'真实 (g={g_true})')
    ax.semilogy(theta_deg, hg_phase(mu_meas, g_inverted_ls), 'r--', linewidth=2, 
                label=f'最小二乘 (g={g_inverted_ls:.3f})')
    ax.set_xlabel('散射角 (度)', fontsize=11)
    ax.set_ylabel('相函数值', fontsize=11)
    ax.set_title('相函数反演结果对比', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 残差分析
    ax = axes[0, 1]
    residual_ls = phase_meas - hg_phase(mu_meas, g_inverted_ls)
    residual_mle = phase_meas - hg_phase(mu_meas, g_inverted_mle)
    residual_leg = phase_meas - hg_phase(mu_meas, g_inverted_leg)
    
    ax.plot(theta_deg, residual_ls, 'o-', label='最小二乘', linewidth=2)
    ax.plot(theta_deg, residual_mle, 's-', label='最大似然', linewidth=2)
    ax.plot(theta_deg, residual_leg, '^-', label='Legendre', linewidth=2)
    ax.axhline(y=0, color='k', linestyle='--', alpha=0.5)
    ax.set_xlabel('散射角 (度)', fontsize=11)
    ax.set_ylabel('残差', fontsize=11)
    ax.set_title('残差分析', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 不同噪声水平下的反演精度
    ax = axes[0, 2]
    noise_levels = np.linspace(0.01, 0.2, 10)
    errors_ls = []
    errors_mle = []
    
    for nl in noise_levels:
        phase_noisy = phase_true * (1 + nl * np.random.randn(len(theta_meas)))
        phase_noisy = np.maximum(phase_noisy, 0.1)
        
        result = optimize.least_squares(residual, x0=[0.5], args=(mu_meas, phase_noisy),
                                        bounds=([-0.99], [0.99]))
        errors_ls.append(abs(result.x[0] - g_true) / g_true * 100)
        
        result_mle = optimize.minimize(neg_log_likelihood, x0=[0.5], 
                                       args=(mu_meas, phase_noisy),
                                       bounds=[(-0.99, 0.99)])
        errors_mle.append(abs(result_mle.x[0] - g_true) / g_true * 100)
    
    ax.plot(noise_levels*100, errors_ls, 'o-', linewidth=2, label='最小二乘')
    ax.plot(noise_levels*100, errors_mle, 's-', linewidth=2, label='最大似然')
    ax.set_xlabel('噪声水平 (%)', fontsize=11)
    ax.set_ylabel('相对误差 (%)', fontsize=11)
    ax.set_title('反演精度 vs 噪声水平', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 4. 不同测量角度数的影响
    ax = axes[1, 0]
    n_angles_list = [10, 18, 36, 72, 180]
    errors_nangles = []
    
    for n_ang in n_angles_list:
        theta_n = np.linspace(0, np.pi, n_ang)
        mu_n = np.cos(theta_n)
        phase_n = hg_phase(mu_n, g_true) * (1 + noise_level * np.random.randn(n_ang))
        phase_n = np.maximum(phase_n, 0.1)
        
        result = optimize.least_squares(residual, x0=[0.5], args=(mu_n, phase_n),
                                        bounds=([-0.99], [0.99]))
        errors_nangles.append(abs(result.x[0] - g_true) / g_true * 100)
    
    ax.semilogy(n_angles_list, errors_nangles, 'o-', linewidth=2, markersize=8, color='purple')
    ax.set_xlabel('测量角度数', fontsize=11)
    ax.set_ylabel('相对误差 (%)', fontsize=11)
    ax.set_title('反演精度 vs 测量角度数', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 5. 参数空间搜索
    ax = axes[1, 1]
    g_search = np.linspace(-0.9, 0.9, 100)
    chi2_values = []
    
    for g_test in g_search:
        phase_test = hg_phase(mu_meas, g_test)
        chi2 = np.sum(((phase_meas - phase_test) / (noise_level * phase_test))**2)
        chi2_values.append(chi2)
    
    ax.plot(g_search, chi2_values, 'b-', linewidth=2)
    ax.axvline(x=g_true, color='g', linestyle='--', linewidth=2, label=f'真实值 g={g_true}')
    ax.axvline(x=g_inverted_ls, color='r', linestyle='--', linewidth=2, 
               label=f'反演值 g={g_inverted_ls:.3f}')
    ax.set_xlabel('不对称因子 g', fontsize=11)
    ax.set_ylabel('χ²', fontsize=11)
    ax.set_title('参数空间搜索', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 6. 反演结果汇总
    ax = axes[1, 2]
    methods = ['真实值', '最小二乘', '最大似然', 'Legendre']
    g_values = [g_true, g_inverted_ls, g_inverted_mle, g_inverted_leg]
    colors = ['green', 'red', 'blue', 'orange']
    
    bars = ax.bar(methods, g_values, color=colors, alpha=0.7, edgecolor='black')
    ax.axhline(y=g_true, color='g', linestyle='--', alpha=0.5)
    ax.set_ylabel('不对称因子 g', fontsize=11)
    ax.set_title('反演结果汇总', fontsize=12, fontweight='bold')
    ax.set_ylim(0, 1)
    ax.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar, val in zip(bars, g_values):
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.4f}', ha='center', va='bottom', fontsize=10)
    
    plt.tight_layout()
    plt.savefig('case4_inversion.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("  结果已保存至: case4_inversion.png")
    
    # 创建动画：反演收敛过程
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax1 = axes[0]
    line_fit, = ax1.plot([], [], 'r-', linewidth=2, label='当前拟合')
    ax1.semilogy(theta_deg, phase_meas, 'ko', markersize=6, label='测量数据')
    ax1.semilogy(theta_deg, phase_true, 'b--', linewidth=2, label='真实值')
    ax1.set_xlabel('散射角 (度)', fontsize=11)
    ax1.set_ylabel('相函数值', fontsize=11)
    ax1.set_title('反演收敛过程', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[1]
    line_error, = ax2.plot([], [], 'b-', linewidth=2)
    ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
    ax2.set_xlabel('迭代次数', fontsize=11)
    ax2.set_ylabel('相对误差 (%)', fontsize=11)
    ax2.set_title('误差收敛', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 模拟迭代过程
    g_iterations = np.linspace(0.1, g_inverted_ls, 50)
    error_history = []
    
    def animate_inversion(frame):
        g_current = g_iterations[frame]
        phase_current = hg_phase(mu_meas, g_current)
        
        line_fit.set_data(theta_deg, phase_current)
        
        error = abs(g_current - g_true) / g_true * 100
        error_history.append(error)
        
        line_error.set_data(range(len(error_history)), error_history)
        ax2.set_xlim(0, 50)
        ax2.set_ylim(0, max(error_history) * 1.1)
        
        ax1.set_title(f'反演收敛过程 (g={g_current:.3f})', fontsize=12, fontweight='bold')
        
        return line_fit, line_error
    
    anim = animation.FuncAnimation(fig, animate_inversion, frames=50, interval=100, blit=True)
    anim.save('case4_inversion.gif', writer='pillow', fps=10, dpi=100)
    plt.close()
    
    print("  动画已保存至: case4_inversion.gif")
    
    return g_inverted_ls, g_inverted_mle, g_inverted_leg


# ============================================================================
# 主程序
# ============================================================================

def main():
    """
    主程序：运行所有仿真案例
    """
    print("\n" + "=" * 70)
    print("主题063：各向异性散射介质辐射仿真")
    print("=" * 70)
    print("本程序包含四个仿真案例：")
    print("  1. Henyey-Greenstein相函数特性分析")
    print("  2. 大气辐射传输中的各向异性散射")
    print("  3. 生物组织中的光传输")
    print("  4. 散射相函数反演")
    print("=" * 70 + "\n")
    
    # 运行案例1
    try:
        result1 = case1_phase_function()
        print("\n案例1完成！\n")
    except Exception as e:
        print(f"案例1运行出错: {e}")
        import traceback
        traceback.print_exc()
    
    # 运行案例2
    try:
        result2 = case2_atmosphere()
        print("\n案例2完成！\n")
    except Exception as e:
        print(f"案例2运行出错: {e}")
        import traceback
        traceback.print_exc()
    
    # 运行案例3
    try:
        result3 = case3_tissue()
        print("\n案例3完成！\n")
    except Exception as e:
        print(f"案例3运行出错: {e}")
        import traceback
        traceback.print_exc()
    
    # 运行案例4
    try:
        result4 = case4_inversion()
        print("\n案例4完成！\n")
    except Exception as e:
        print(f"案例4运行出错: {e}")
        import traceback
        traceback.print_exc()
    
    print("\n" + "=" * 70)
    print("所有仿真案例已完成！")
    print("生成的文件：")
    print("  - case1_phase_function.png")
    print("  - case1_phase_function.gif")
    print("  - case2_atmosphere.png")
    print("  - case2_atmosphere.gif")
    print("  - case3_tissue.png")
    print("  - case3_tissue.gif")
    print("  - case4_inversion.png")
    print("  - case4_inversion.gif")
    print("=" * 70)


if __name__ == "__main__":
    main()