模拟分析与结构a

# 导入相关的库
import numpy as np
import scipy.constants as const
import scipy.stats as stats
import matplotlib.pyplot as plt

# 定义物理参数
m = const.electron_mass # 电子质量
e = const.e # 电子电荷量
c = const.c # 光速
hbar = const.hbar # 约化普朗克常量
alpha = const.fine_structure # 精细结构常数
omega0 = 1e-6 * const.e # 激光光子能量 (1 MeV)
a0 = 200 # 激光场归一化强度
gamma = 1000 # 电子洛伦兹因子
beta = np.sqrt(1 - 1 / gamma**2) # 电子速度 / c（光速归一化速度）

# Define inverse transform sampling method
def inverse_transform_sampling(pdf, cdf_inv, n_samples):
    # pdf: probability density function
    # cdf_inv: inverse cumulative distribution function
    # n_samples: number of samples to generate
    u = np.random.uniform(size=n_samples) # generate uniform random numbers in [0, 1]
    x = cdf_inv(u) # apply inverse cdf to get samples from target distribution
    return x

# Define Compton scattering functions
def photon_energy_after_scattering(omega0, gamma, beta, theta):
    # omega0: initial photon energy
    # gamma: electron Lorentz factor
    # beta: electron velocity / c
    # theta: scattering angle in electron rest frame
    omega = omega0 / (gamma * (1 - beta * np.cos(theta))) # final photon energy
    return omega

def photon_angle_after_scattering(omega0, gamma, beta, theta):
    # omega0: initial photon energy
    # gamma: electron Lorentz factor
    # beta: electron velocity / c
    # theta: scattering angle in electron rest frame
    omega = photon_energy_after_scattering(omega0, gamma, beta, theta) # final photon energy
    phi = np.arcsin(np.sin(theta) / (gamma * (1 + omega / (gamma * omega0)))) # scattering angle in lab frame
    return phi

def photon_polarization_after_scattering(omega0, gamma, beta, theta):
    # omega0: initial photon energy
    # gamma: electron Lorentz factor
    # beta: electron velocity / c
    # theta: scattering angle in electron rest frame
    omega = photon_energy_after_scattering(omega0, gamma, beta, theta) # final photon energy
    phi = photon_angle_after_scattering(omega0, gamma, beta, theta) # scattering angle in lab frame
    chi = np.arctan2(np.sin(theta), gamma * (np.cos(theta) - beta)) # azimuthal angle of polarization vector
    epsilon = np.sqrt((omega0 / omega)**2 * np.sin(phi)**2 + np.cos(phi)**2) # degree of linear polarization
    return chi, epsilon

def electron_spin_after_scattering(omega0, gamma, beta, theta):
    # omega0: initial photon energy
    # gamma: electron Lorentz factor
    # beta: electron velocity / c
    # theta: scattering angle in electron rest frame
    omega = photon_energy_after_scattering(omega0, gamma, beta, theta) # final photon energy
    phi = photon_angle_after_scattering(omega0, gamma, beta, theta) # scattering angle in lab frame
    s = np.random.choice([-1, 1]) # spin projection along local quantization axis (random choice)
    psi = np.arctan2(-beta * np.sin(theta), 1 - beta * np.cos(theta)) # angle between local quantization axis and z-axis
    sz = s * np.cos(psi) # spin projection along z-axis
    return s, sz

# Define Klein-Nishina formula for differential cross section
def klein_nishina(omega0, gamma, beta, theta):
    # omega0: initial photon energy
    # gamma: electron Lorentz factor
    # beta: electron velocity / c
    # theta: scattering angle in electron rest frame
    omega = photon_energy_after_scattering(omega0, gamma, beta, theta) # final photon energy
    r0 = e**2 / (m * c**2) # classical electron radius
    f = omega / omega0 # energy fraction
    g = 1 / (1 + f * (1 - np.cos(theta))) # kinematic factor
    dsigma = 2 * np.pi * r0**2 * (f / g)**2 * (f + 1 / f - np.sin(theta)**2) # differential cross section
    return dsigma

# Define inverse cumulative distribution function for Klein-Nishina formula
def klein_nishina_cdf_inv(omega0, gamma, beta, u):
    # omega0: initial photon energy
    # gamma: electron Lorentz factor
    # beta: electron velocity / c
    # u: uniform random number in [0, 1]
    r0 = e**2 / (m * c**2) # classical electron radius
    sigma0 = 2 * np.pi * r0**2 / (gamma**2 * (1 - beta)**2) # total cross section
    fmax = 4 * gamma**2 / (1 + 2 * gamma)**2 # maximum energy fraction
    gmin = 1 / (1 + fmax) # minimum kinematic factor
    gmax = gamma * (1 + beta) # maximum kinematic factor
    def F(g): # cumulative distribution function as a function of g
        f = 1 / g # energy fraction
        x = 1 - f * (1 - np.cos(np.pi)) # auxiliary variable
        y = 1 - f * (1 - np.cos(0)) # auxiliary variable
        return sigma0 * (np.log(gmax / gmin) + (fmax + 1 / fmax - 2) * np.log(y / x) + (y - x) / 2) / sigma0 - u
    g = stats.brentq(F, gmin, gmax) # find the root of F(g) using Brent's method
    theta = np.arccos(1 - 1 / (f * g) + 1 / f) # scattering angle in electron rest frame
    return theta

# Define number of samples and arrays for storing data
n_samples = 10000 # number of scattering events to simulate
omega_array = np.zeros(n_samples) # array for final photon energy
phi_array = np.zeros(n_samples) # array for scattering angle in lab frame
chi_array = np.zeros(n_samples) # array for azimuthal angle of polarization vector
epsilon_array = np.zeros(n_samples) # array for degree of linear polarization
s_array = np.zeros(n_samples) # array for spin projection along local quantization axis
sz_array = np.zeros(n_samples) # array for spin projection along z-axis

# Define loop for simulating Compton scattering events using inverse transform sampling method and Compton scattering functions
for i in range(n_samples):
    theta = inverse_transform_sampling(klein_nishina, klein_nishina_cdf_inv, 1)[0] # generate random scattering angle from Klein-Nishina formula using inverse transform sampling method
    omega_array[i] = photon_energy_after_scattering(omega0, gamma, beta, theta) # calculate final photon energy
    phi_array[i] = photon_angle_after_scattering(omega0, gamma, beta, theta) # calculate scattering angle in lab frame
    chi_array[i], epsilon_array[i] = photon_polarization_after_scattering(omega0, gamma, beta, theta) # calculate azimuthal angle and degree of linear polarization
    s_array[i], sz_array[i] = electron_spin_after_scattering(omega0, gamma, beta, theta) # calculate spin projection along local quantization axis and z-axis

# Analyze and visualize