use std::f64::consts::PI;
pub const NUM_OPERATORS: usize = 15;
#[derive(Clone, Debug, PartialEq)]
pub struct EftCoefficients {
pub c_p: [f64; NUM_OPERATORS],
pub c_n: [f64; NUM_OPERATORS],
}
impl EftCoefficients {
pub fn new() -> Self {
Self {
c_p: [0.0; NUM_OPERATORS],
c_n: [0.0; NUM_OPERATORS],
}
}
pub fn isoscalar(&self, op: usize) -> f64 {
if op == 0 || op > NUM_OPERATORS {
return 0.0;
}
(self.c_p[op - 1] + self.c_n[op - 1]) / 2.0
}
pub fn isovector(&self, op: usize) -> f64 {
if op == 0 || op > NUM_OPERATORS {
return 0.0;
}
(self.c_p[op - 1] - self.c_n[op - 1]) / 2.0
}
pub fn set_si(mut self, c: f64) -> Self {
self.c_p[0] = c;
self.c_n[0] = c;
self
}
pub fn set_sd(mut self, c_p: f64, c_n: f64) -> Self {
self.c_p[3] = c_p;
self.c_n[3] = c_n;
self
}
}
impl Default for EftCoefficients {
fn default() -> Self {
Self::new()
}
}
pub fn standard_si_cross_section_cm2(c1: f64, mu_gev: f64) -> f64 {
mu_gev * mu_gev * c1 * c1 / PI * 1.782_661_92e-27 * 1.782_661_92e-27 * 1e4
}
pub fn standard_sd_cross_section_cm2(c4: f64, mu_gev: f64, j_nucleus: f64) -> f64 {
3.0 * mu_gev * mu_gev * c4 * c4 * (j_nucleus + 1.0) / (PI * j_nucleus)
* 1.782_661_92e-27
* 1.782_661_92e-27
* 1e4
}
pub fn nuclear_response_m(q_fm_inv: f64, a: u32, z: u32) -> f64 {
let q_r = q_fm_inv * nuclear_radius_fm(a);
let f = helm_form_factor_sq(q_r, a);
let n = a - z;
(z as f64 + n as f64).powi(2) * f
}
pub fn nuclear_response_sigma_prime(
q_fm_inv: f64,
a: u32,
spin_proton: f64,
spin_neutron: f64,
j_nucleus: f64,
) -> f64 {
let r = nuclear_radius_fm(a);
let q_r = q_fm_inv * r;
let suppression = (-q_r * q_r / 4.0).exp();
let structure = (j_nucleus + 1.0) / (3.0 * j_nucleus) * (spin_proton + spin_neutron).powi(2);
structure * suppression
}
pub fn nuclear_response_sigma_double_prime(
q_fm_inv: f64,
a: u32,
spin_proton: f64,
spin_neutron: f64,
j_nucleus: f64,
) -> f64 {
let r = nuclear_radius_fm(a);
let q_r = q_fm_inv * r;
let suppression = (-q_r * q_r / 4.0).exp();
let q2_over_mn2 = (q_fm_inv * 0.197).powi(2) / (0.938 * 0.938);
let structure = (j_nucleus + 1.0) / (12.0 * j_nucleus) * (spin_proton - spin_neutron).powi(2);
structure * suppression * q2_over_mn2
}
pub fn nuclear_response_delta(q_fm_inv: f64, a: u32, z: u32) -> f64 {
let q_r = q_fm_inv * nuclear_radius_fm(a);
let f = helm_form_factor_sq(q_r, a);
let q2_over_mn2 = (q_fm_inv * 0.197).powi(2) / (0.938 * 0.938);
(z as f64).powi(2) * f * q2_over_mn2
}
pub fn nuclear_response_phi_pp(q_fm_inv: f64, a: u32, z: u32) -> f64 {
let q_r = q_fm_inv * nuclear_radius_fm(a);
let f = helm_form_factor_sq(q_r, a);
let n = a - z;
((z * (z - 1)) as f64 + (n * (n - 1)) as f64) * f / 4.0
}
#[derive(Clone, Debug, PartialEq)]
pub struct NuclearTarget {
pub a: u32,
pub z: u32,
pub j_nucleus: f64,
pub sp: f64,
pub sn: f64,
}
#[derive(Clone, Debug, PartialEq)]
pub struct DmProperties {
pub m_dm_gev: f64,
pub j_chi: f64,
pub rho_dm_gev_cm3: f64,
pub v_mean_km_s: f64,
}
fn helm_form_factor_sq(q_r: f64, a: u32) -> f64 {
if q_r.abs() < 1e-10 {
return 1.0;
}
let s = 0.9;
let r_n = nuclear_radius_fm(a);
let j1 = (q_r.sin() - q_r * q_r.cos()) / (q_r * q_r * q_r);
let f = 3.0 * j1 * (-0.5 * (q_r * s / r_n).powi(2)).exp();
f * f
}
fn nuclear_radius_fm(a: u32) -> f64 {
let c = 1.23 * (a as f64).powf(1.0 / 3.0) - 0.6;
let a_skin = 0.52;
(c * c + 7.0 / 3.0 * PI * PI * a_skin * a_skin - 5.0 * 0.9 * 0.9).sqrt()
}
pub fn operator_matrix_element_sq(
op: usize,
q_fm_inv: f64,
v_perp_sq: f64,
j_chi: f64,
target: &NuclearTarget,
) -> f64 {
let m_response = nuclear_response_m(q_fm_inv, target.a, target.z);
let sigma_p =
nuclear_response_sigma_prime(q_fm_inv, target.a, target.sp, target.sn, target.j_nucleus);
let sigma_pp = nuclear_response_sigma_double_prime(
q_fm_inv,
target.a,
target.sp,
target.sn,
target.j_nucleus,
);
let delta = nuclear_response_delta(q_fm_inv, target.a, target.z);
let phi_pp = nuclear_response_phi_pp(q_fm_inv, target.a, target.z);
let s_chi = j_chi * (j_chi + 1.0) / 3.0;
match op {
1 => m_response,
3 => sigma_pp * v_perp_sq,
4 => sigma_p * s_chi,
5 => m_response * v_perp_sq + delta,
6 => sigma_pp * s_chi,
7 => sigma_p,
8 => m_response * v_perp_sq,
9 => sigma_p * s_chi,
10 => sigma_pp,
11 => m_response,
12 => phi_pp * s_chi,
13 => sigma_pp * s_chi * v_perp_sq,
14 => sigma_p,
15 => phi_pp,
_ => 0.0,
}
}
pub fn eft_differential_rate(
coeffs: &EftCoefficients,
e_r_kev: f64,
m_nucleus_gev: f64,
target: &NuclearTarget,
dm: &DmProperties,
) -> f64 {
let mu = dm.m_dm_gev * m_nucleus_gev / (dm.m_dm_gev + m_nucleus_gev);
let e_r_gev = e_r_kev * 1e-6;
let q_gev = (2.0 * m_nucleus_gev * e_r_gev).sqrt();
let q_fm_inv = q_gev / 0.197;
let v_min = (m_nucleus_gev * e_r_gev / (2.0 * mu * mu)).sqrt() * 3e8 / 1e3;
let v_mean = dm.v_mean_km_s;
if v_min > v_mean * 3.0 {
return 0.0;
}
let eta = (-v_min * v_min / (v_mean * v_mean)).exp() / (v_mean * 1e3);
let v_perp_sq = (v_mean * 1e3).powi(2) - (q_gev / (2.0 * mu * 1.782_661_92e-27)).powi(2);
let v_perp_sq = v_perp_sq.max(0.0);
let mut rate = 0.0;
for op in 1..=NUM_OPERATORS {
let c_eff = coeffs.c_p[op - 1] * target.z as f64
+ coeffs.c_n[op - 1] * (target.a - target.z) as f64;
if c_eff.abs() < 1e-50 {
continue;
}
let me_sq = operator_matrix_element_sq(op, q_fm_inv, v_perp_sq, dm.j_chi, target);
rate += c_eff * c_eff * me_sq;
}
let n_dm = dm.rho_dm_gev_cm3 / dm.m_dm_gev;
rate * n_dm * eta / (2.0 * dm.m_dm_gev * 1.782_661_92e-27)
}