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

蒙特卡洛光学应用

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

计算物质中电子的非弹性平均自由程（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

PreviousPython密码线性方程组和频率分析 NextPython电网静态和准静态分析及自动化优化

Last updated 2 years ago

Was this helpful?

# -*- 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()