🥥Python蒙特卡洛算法光学数据模型非弹性散射模拟

关键词

Python | Matplotlib | Numpy | 蒙特卡洛方法 | 光学数据 | 电子散射 | 动能分布

梗概

蒙特卡洛光学应用

使用蒙特卡洛的驱动原因：

• 计算物质中电子的非弹性平均自由程（IMFP）

• 物理过程结果的预测

蒙特卡罗模型的可靠性通常通过与实验结果的比较来检验； 不幸的是，有些量只能间接测量或存在一些实验问题。 基本问题在于

• 在选择合适的数量时，可以直接测量，并且可以将计算结果与它们进行比较；

• 使检查有效的这些量的范围最大化（电子能量、原子序数、表面膜的厚度等）

最初使用这种方法时，由于计算机速度较慢，通常采用多次散射计算，路径较长的部分同时利用散射效应的平均进行模拟。 如今，采用了所谓的单散射模型，其中每个散射事件都是单独计算的。 在 Monte-Carlo 代码中，可以使用计算所需的公式和值表。 由于需要在表格给出的值之间进行插值，因此首选公式。

在用于电子显微镜、光谱学和微量分析的能量范围内（即通常为 0.1 至 300 keV），蒙特卡洛代码中使用各种模型来描述电子与原子的基本相互作用——弹性和非弹性碰撞以及电子与材料的其他相互作用 。

电子与材料的物理相互作用公式

粒子 $A_1$ 与原子核 $A_2$ 的相互作用一般用公式来描述

$A_{1}+A_{2} \Rightarrow A_{3}+A_{4}+Q$

其中 $A_3$ 是出射粒子，$A_4$ 是产生的原子核，$Q$ 是发射能量。

非弹性散射

对于简单的解决方案，可以使用连续减速（CSD）的简单模型。 在更精确的单次散射计算中，我们需要知道非弹性截面。 由于带电粒子与具有各种介电结构的物质相互作用的复杂性，非弹性微分截面的计算通常稍微复杂一些。 此外，非弹性截面是双重微分的——通过动量变化（取决于散射角）和电子能量损失。 对于蒙特卡洛模拟，我们需要计算给定电子能量和元素的非弹性平均自由程 (IMFP) 和阻止能力。

模型描述

以下文字描述了用于模拟光电子散射的模型。 该模型已在算法中实现，并用于确定散射统计数据（例如沿其路径非弹性散射 n 次的电子数量）、路径长度分布以及模拟光电子穿过散射介质的角度分布。 该方法已应用于气相散射介质，以模拟在近环境压力 XPS 测量过程中遇到的情况； 但是，它们也可以应用于模拟固体中的电子散射，并且可以在线获得示例。 在此，我们提供了算法中使用的所有方程和方法的推导。

拟合散射截面函数

表面能量损失

计算细节和精度估计

蒙特卡罗计算与实验的比较

光学电子类

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

class Electron():
  
    
    def __init__(self, source, **kwargs):
     
        self.mass = 9.109383E-31 # the electron mass in kg
        self.J_eV = 1.602176E-19 # conversion of eV to Joules. has units J/eV
        self.vector = np.zeros((9))
        self.pathlength = 0 # This stores to total length the electron has
        # travelled
        self.source = source
        #self.kinetic_energy = self.source.getKE()
        self.theta = self.source.theta
        self.phi = self.source.phi
        self.initCoords()
        self.initVelocity()

    def _convertKinEnToSpeed(self, kinetic_energy: float) -> float:
      
        speed = 1E+9 * np.sqrt(2 * kinetic_energy * 
                                     self.J_eV / self.mass)
        return speed
        
    def initCoords(self, *args):
      
        self.vector[0:3] = self.source.getRandPosition()
        #self.vector[0:3] = self.source.getSlicePosition(*args)
        
    def initVelocity(self):
      
        KE = self.source.getKE()
        self.kinetic_energy = KE
        initial_speed = self._convertKinEnToSpeed(KE)
        theta = self.theta.getAngle()
        phi = self.phi.getAngle()
        vx = initial_speed * np.sin(theta) * np.cos(phi)
        vy = initial_speed * np.sin(theta) * np.sin(phi)
        vz = initial_speed * np.cos(theta)
        self.vector[3:6] = np.array([vx,vy,vz])
        
    def updateKineticEnergy(self, delta_KE):
       
        new_KE = self.kineticEnergy() - delta_KE
        if new_KE < 0:
            self.kinetic_energy = 0
        else:
            self.kinetic_energy = new_KE
        self.changeSpeed()
        
    def changeSpeed(self):
      
        old_v = self.vector[3:6]
  
        speed = 1E+9 * np.sqrt(2 * self.kinetic_energy * 
                                     self.J_eV / self.mass)
        new_v = (old_v / np.sqrt(old_v[0]**2 + old_v[1]**2 + old_v[2]**2)
            * speed)
     
        self.vector[3:6] = new_v
        
    def kineticEnergy(self) -> float:
   
        speed = np.linalg.norm(self.vector[3:6]) / 1E+9 
      
        KE = (1/2 * self.mass * speed**2) / self.J_eV
        return KE
       

if __name__ == '__main__':
    source = Disc(100,1)
    vals = []
    e = Electron(source)
    for i in range(2000):
        e.initVelocity()
        vals+=[list(e.vector)]
    vals = np.array(vals)

    
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111, projection='3d')
    xs = vals[:,3]
    ys = vals[:,4]
    zs = vals[:,5]
    ax.scatter(xs, ys, zs, s=20, marker = '.')
  
    ax.set_xlabel('X Label')
    ax.set_ylabel('Y Label')
    ax.set_zlabel('Z Label')
    ax.view_init(10, 45)  
    
    plt.show()        

if __name__ == '__main__':
    source = Disc(100,1) 
    e = Electron(source) 
    e.initVelocity() 
    v1 = e.kineticEnergy() 
    e.changeSpeed(10) 
   
    v2 = e.kineticEnergy() 

源代码

🏈指点迷津 | Brief