/
plasma
/
compton_ff_yield
Обзор
Документация
Войти
/
plasma
/
compton_ff_yield
Код
Запросы
0
Задачи
Вики
Пакеты
0
Релизы
0
CI/CD
Аналитика
Безопасность
master
model.py
116 строк
4 KB
Sergey Rykovanov
Initial commit
13 июл 2026, 19:42
13 июл 2026, 19:42
65c7569
Код
Авторство
О чём код?
from dataclasses import dataclass, replace from pathlib import Path import numpy as np import yaml from scipy.constants import c, elementary_charge, epsilon_0, hbar, m_e, physical_constants from scipy.integrate import quad from units import gamma_from_energy SIGMA_T = physical_constants["Thomson cross section"][0] @dataclass class Case: electron_energy_MeV: float charge_pC: float epsn_m_rad: float beta_star_m: float sigma_le_m: float laser_energy_J: float wavelength_m: float a0_peak: float sigma_p0_m: float beta_ff: float = 0.0 sigma_J_perp_m: float = 0.0 sigma_JZ_m: float = 0.0 sigma_Jt_s: float = 0.0 scattering_model: str = "kn" def load_case(path=None): path = Path(__file__).with_name("case_paper.yaml") if path is None else Path(path) with Path(path).open() as stream: return Case(**yaml.safe_load(stream)) def sigma_kn(kappa): if kappa < 1e-3: return SIGMA_T * (1 - 2*kappa + 26*kappa**2/5 - 133*kappa**3/10) log_term = np.log1p(2*kappa) return SIGMA_T * 0.75 * ( (1 + kappa) / kappa**3 * (2*kappa*(1 + kappa)/(1 + 2*kappa) - log_term) + log_term/(2*kappa) - (1 + 3*kappa)/(1 + 2*kappa)**2 ) def parameters(case): if not -0.999 < case.beta_ff <= 1: raise ValueError("beta_ff must be in (-1, 1].") if case.scattering_model not in {"thomson", "kn"}: raise ValueError("scattering_model must be 'thomson' or 'kn'.") gamma = gamma_from_energy(case.electron_energy_MeV) beta = np.sqrt(1 - gamma**-2) eps = case.epsn_m_rad / (gamma * beta) sigma_e2 = eps * case.beta_star_m omega = 2*np.pi*c / case.wavelength_m photon_energy = hbar * omega photons = case.laser_energy_J / photon_energy electrons = case.charge_pC * 1e-12 / elementary_charge field = case.a0_peak * m_e * c * omega / elementary_charge intensity = epsilon_0 * c * field**2 / 2 photon_density = intensity / (c * photon_energy) sigma_lp = photons / ((2*np.pi)**1.5 * photon_density * case.sigma_p0_m**2) rayleigh = 4*np.pi * case.sigma_p0_m**2 / case.wavelength_m kappa_kn = gamma * (1 + beta) * photon_energy / (m_e*c**2) sigma = SIGMA_T if case.scattering_model == "thomson" else sigma_kn(kappa_kn) return { "gamma": gamma, "electrons": electrons, "photons": photons, "sigma_e2": sigma_e2, "sigma_lp": sigma_lp, "rayleigh": rayleigh, "kappa_e": sigma_e2 / case.beta_star_m**2, "kappa_p": case.sigma_p0_m**2 / rayleigh**2, "h": (1 - case.beta_ff) / (2*(1 + case.beta_ff)), "g": case.beta_ff / (1 + case.beta_ff), "sigma_t": c * case.sigma_Jt_s, "overlap2": sigma_e2 + case.sigma_p0_m**2 + case.sigma_J_perp_m**2, "cross_section": sigma, "pulse_duration_fs": sigma_lp / c * 1e15, } def yield_mean(case, rtol=2e-9, kappa_e=None, kappa_p=None): p = parameters(case) ke = p["kappa_e"] if kappa_e is None else kappa_e kp = p["kappa_p"] if kappa_p is None else kappa_p def integrand(u): s = u / p["overlap2"] K_e, K_p = s*ke, s*kp longitudinal = 1 + 2*K_p*case.sigma_JZ_m**2 lambda_p = K_p / longitudinal A_pp = 1/p["sigma_lp"]**2 + K_e/2 + 2*p["h"]**2*lambda_p A_pe = K_e/2 + p["h"]*lambda_p A_ee = 1/case.sigma_le_m**2 + K_e/2 + lambda_p/2 determinant = A_pp*A_ee - A_pe**2 if p["sigma_t"]: vectors = np.array([[1, 0, -1], [0, 1, 0], [1, 1, 0], [p["h"], .5, p["g"]], [0, 0, 1]]) weights = [1/p["sigma_lp"]**2, 1/case.sigma_le_m**2, K_e/2, 2*lambda_p, 1/p["sigma_t"]**2] matrix = sum(weight*np.outer(vector, vector) for weight, vector in zip(weights, vectors)) determinant = p["sigma_t"]**2 * np.linalg.det(matrix) return np.exp(-u) / (p["overlap2"] * np.sqrt(longitudinal * determinant)) integral, _ = quad(integrand, 0, np.inf, epsabs=0, epsrel=rtol, limit=250) return p["cross_section"] * p["electrons"] * p["photons"] * integral / (2*np.pi * case.sigma_le_m * p["sigma_lp"]) def ideal_case(case): return replace(case, sigma_J_perp_m=0, sigma_JZ_m=0, sigma_Jt_s=0)