MC模型

一个 3D 蒙特卡罗散射代码，用于跟踪平面平行同质板中的单个地块

Petrich 等人已经描述了模型背后的理论、它的许多参数和模型验证。(2012), https://doi.org/10.1016/j.coldregions.2011.12.004

Monte Carlo 代码本身以 C 语言实现以提高速度，并与 Python 接口集成以方便使用。

许可证：Apache 2.0。

安装

要求

• Python 2.7、3.5、3.6、3.7。应该可以在 Linux、MacOS 和 Windows 上运行。
• 麻木的
• 从源代码编译时，一个 C 编译器，通常是 setuptools 和 wheel。

来自 PyPI

运行pip install mcmodel。完毕。稍后，使用 升级pip install -U mcmodel或使用 卸载 pip uninstall mcmodel。

无需安装即可使用

从https://github.com/cpetrich/MCmodel获取最新的源代码 运行python setup.py build并将目录内容复制built/lib.[...] 到与自己的建模代码相同的目录中。

使用 pip 安装

从https://github.com/cpetrich/MCmodel获取最新的源代码 将在本地构建一个轮子并使用 pip 安装。构建模块并使用python setup.py bdist_wheel, cd dist, . 创建一个 wheel 文件pip install mcmodel-[...].whl。稍后要卸载该软件包，请运行 pip uninstall mcmodel.

测试

要测试构建过程的结果，可以使用python setup.py test.

限制

该代码目前不可重入，因为它使用全局状态来生成随机数生成器、相位函数和保存有关最新粒子轨迹的详细信息的缓冲区。即，每个进程只有一个模型实例。

特殊情况——在散射模型中——“无散射”( w0 == 0) 仅用于k == 0（即无吸收）。

用法

通常的操作模式是：

0. 种子随机数生成器
1. 定义散射相位函数
2. 调用具有平板光学特性的散射码，以及一包光的位置和方向
3. 存储包裹的入口和出口位置和角度以及介质内平面交叉点的信息
4. 如果需要，存储轻包裹的完整轨迹
5. 继续（2）几千或几百万次。
6. 从磁盘读取模拟结果，并确定由于吸收和菲涅耳反射引起的包裹强度
7. 根据需要分析结果

播种 Monte Carlo 模型的随机数生成器的示例是

import os
import mcmodel
mcmmodel.set_seed(int.from_bytes(os.urandom(4),'little'))

存在一个辅助函数来生成常见的散射相位函数，例如

import mcmodel_util
pf = mcmodel_util.make_phase_function('Henyey-Greenstein', (0.98, 0))

相位函数被转移到模型代码中，如

mcmodel.define_phase_function(pf['lookup_cdf'],pf['lookup_phi'])

随机数生成器通常用于生成传入分布（该随机数生成器与蒙特卡罗模型使用的随机数生成器是分开的，也应该是种子）。

import random
random.seed()

来自 Lambert 发射器角分布的样本（即，积雪覆盖的地面上的阴天）是

import math

elevation_angle = math.arcsin(random.uniform(-1., 0))
azimuth = random.uniform(-math.pi, math.pi)

请注意，蒙特卡洛代码适用于仰角。平板上表面的垂直入射光束的仰角为 -pi/2。

域参数和入射角在字典中传递给蒙特卡洛代码，例如

in = {'angle': elevation_angle,
      'azimuth': azimuth,
      'thickness': 1, # slab thickness
      'k': 1}         # densiy of microscopic interactions used by the model

带有如下所述的附加可选输入参数。许多输出数据是通过更新传递给散射模拟的字典获得的。计算单个包裹调用的轨迹

out = {}
mcmodel.simulate(in, out)

其中出口位置和角度是

exit_position = (out['x'], out['y'], out['z'])
exit_elevation_angle = out['angle']
exit_azimuth = out['azimuth']
path_length = out['path_length']

介质内包裹的衰减是通过减小其幅度来解释的。精确的方法取决于传递到上述模型中的厚度和微观消光系数的定义和单位。然而，新兴地块的强度将类似于

exit_intensity = math.exp(-kappa * path_length)

其中 kappa 是吸收系数。

由于菲涅耳反射（全反射除外）导致的强度降低将根据入射角和出射角以及相对折射率来计算。事后看来，这样做并非易事，并且需要在相对折射率为 1 的情况下运行模拟，并适当缩放输入角分布以考虑非平凡的折射率。这超出了本概述的范围。 [TODO：在未来的版本中，添加一个选项，让 Monte Carlo 模型通过随机选择路径来处理这个问题。]

API

C 代码当前导出四个函数（请参阅上面的用法）：

• mcmodel.set_seed(...)
• mcmodel.define_phase_function(..., ...)
• mcmodel.simulate(..., ...)
• mcmodel.get_last_particle_track()

此外，Python 中还定义了一个辅助函数

• mcmodel_util.make_phase_function(...)

参数

模型域是一个厚度为 的水平板thickness，从 z=0上表面延伸z=-thickness到底部。Parcles 在深度z=z0( -thickness <= z0 <= 0) 和x=0和处注入y=0。分别在下界面和上界面的板坯外侧注入 z0 = -thickness和（即仍会发生折射）。z0 = 0

如果使用各向异性散射系数，即sigma_anisotropy != 0，则垂直和水平散射系数分别为 sigma_v = k * w0和sigma_h = k * w0 * (sigma_anisotropy+1)。

输入字典mcmodel.simulate()

输出字典mcmodel.simulate()

每个平面交叉点的六列条目是：

1. 仰角
2. 方位角
3. 交叉点的 x 坐标
4. 交叉点的 y 坐标
5. 交叉点的 z 坐标
6. 到目前为止的包裹路径长度