diff --git a/Cargo.toml b/Cargo.toml index 8c6213a..1bd4fa3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,6 +28,7 @@ hdf5 = {version = "0.7.1", optional = true} rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} ndarray = {version = "0.14.0", features = ["serde"], optional = true} parry3d-f64 = {optional = true, version="0.2.0"} +cached = {features=["proc_macro"], version="2.0.2"} [dependencies.pyo3] version = "0.19.0" diff --git a/src/bca.rs b/src/bca.rs index e1e300e..a277eb6 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -82,14 +82,6 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate //Choose impact parameters and azimuthal angles for all collisions, and determine mean free path let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, material, options, rng); - #[cfg(feature = "accelerated_ions")] - let distance_to_target = if !material.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z) { - let (x, y, z) = material.geometry.closest_point(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); - ((x - particle_1.pos.x).powi(2) + (y - particle_1.pos.y).powi(2) + (z - particle_1.pos.z).powi(2)).sqrt() - } else { - 0. - }; - let mut total_energy_lost_to_recoils = 0.; let mut total_asymptotic_deflection = 0.; let mut normalized_distance_of_closest_approach = 0.; @@ -159,7 +151,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate let n = material.total_number_density(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); //We just need the lindhard screening length here, so the particular potential is not important - let a: f64 = interactions::screening_length(Za, Zb, InteractionPotential::MOLIERE); + let a: f64 = interactions::lindhard_screening_length_lookup(Za as u64, Zb as u64); let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E; let estimated_range_of_recoils = (reduced_energy.powf(0.3) + 0.1).powi(3)/n/a/a; @@ -181,14 +173,9 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate } //Advance particle in space and track total distance traveled - #[cfg(not(feature = "accelerated_ions"))] let distance_traveled = particle_1.advance( binary_collision_geometries[0].mfp, total_asymptotic_deflection); - #[cfg(feature = "accelerated_ions")] - let distance_traveled = particle_1.advance( - binary_collision_geometries[0].mfp + distance_to_target - material.geometry.get_energy_barrier_thickness(), total_asymptotic_deflection); - //Subtract total energy from all simultaneous collisions and electronic stopping let energy_lost_to_electronic_stopping = bca::subtract_electronic_stopping_energy(&mut particle_1, material, distance_traveled, normalized_distance_of_closest_approach, strong_collision_Z, @@ -235,7 +222,8 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle let E: f64 = particle_1.E; let Ec: f64 = particle_1.Ec; //We just need the Lindhard screening length here, so the particular potential is not important - let a: f64 = interactions::screening_length(Za, Zb, InteractionPotential::MOLIERE); + let a: f64 = interactions::lindhard_screening_length_lookup(Za as u64, Zb as u64); + let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E; //Minimum energy transfer for generating scattering event set to cutoff energy @@ -510,8 +498,8 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: & _ => x0*a*(theta/2.).sin() }; - let psi = theta.sin().atan2(Ma/Mb + theta.cos());//.abs(); - let psi_recoil = theta.sin().atan2(1. - theta.cos());//.abs(); + let psi = theta.sin().atan2(Ma/Mb + theta.cos()); + let psi_recoil = theta.sin().atan2(1. - theta.cos()); let recoil_energy = 4.*(Ma*Mb)/(Ma + Mb).powi(2)*E0*(theta/2.).sin().powi(2); Ok(BinaryCollisionResult::new(theta, psi, psi_recoil, recoil_energy, asymptotic_deflection, x0)) @@ -525,7 +513,9 @@ fn scattering_integral_mw(x: f64, beta: f64, reduced_energy: f64, interaction_po } /// Gauss-Legendre scattering integrand. -fn scattering_function_gl(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: &dyn Fn(f64) -> f64) -> Result { +fn scattering_function_gl(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: F) -> Result + where F: Fn(f64) -> f64 +{ let result = 4.*impact_parameter*u/(r0*(1. - interaction_potential(r0/(1. - u*u))/relative_energy - impact_parameter*impact_parameter*(1. - u*u).powi(2)/r0/r0).sqrt()); if result.is_nan() { @@ -537,7 +527,9 @@ fn scattering_function_gl(u: f64, impact_parameter: f64, r0: f64, relative_energ } /// Gauss-Mehler scattering integrand. -fn scattering_function_gm(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: &dyn Fn(f64) -> f64) -> Result { +fn scattering_function_gm(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: F) -> Result + where F: Fn(f64) -> f64 +{ let result = impact_parameter/r0/(1. - interaction_potential(r0/u)/relative_energy - (impact_parameter*u/r0).powi(2)).sqrt(); if result.is_nan() { @@ -549,23 +541,27 @@ fn scattering_function_gm(u: f64, impact_parameter: f64, r0: f64, relative_energ } /// Compute the scattering integral for a given relative energy, distance of closest approach `r0`, and interaction potential using a Gauss-Mehler, n-point quadrature. -fn scattering_integral_gauss_mehler(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: &dyn Fn(f64) -> f64, n_points: usize) -> f64 { +fn scattering_integral_gauss_mehler(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: F, n_points: usize) -> f64 + where F: Fn(f64) -> f64 + Clone +{ let x: Vec = (1..=n_points).map(|i| ((2.*i as f64 - 1.)/4./n_points as f64*PI).cos()).collect(); let w: Vec = (1..=n_points).map(|i| PI/n_points as f64*((2.*i as f64 - 1.)/4./n_points as f64*PI).sin()).collect(); PI - x.iter().zip(w) - .map(|(&x, w)| w*scattering_function_gm(x, impact_parameter, r0, relative_energy, interaction_potential) + .map(|(&x, w)| w*scattering_function_gm(x, impact_parameter, r0, relative_energy, interaction_potential.clone()) .with_context(|| format!("Numerical error: NaN in Gauss-Mehler scattering integral at x = {} with Er = {} eV and p = {} A.", x, relative_energy/EV, impact_parameter/ANGSTROM)) .unwrap()).sum::() } /// Compute the scattering integral for a given relative energy, distance of closest approach `r0`, and interaction potential using a Gauss-Legendre, 5-point quadrature. -fn scattering_integral_gauss_legendre(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: &dyn Fn(f64) -> f64) -> f64 { +fn scattering_integral_gauss_legendre(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: F) -> f64 + where F: Fn(f64) -> f64 + Clone +{ let x: Vec = [0., -0.538469, 0.538469, -0.90618, 0.90618].iter().map(|x| x/2. + 1./2.).collect(); let w: Vec = [0.568889, 0.478629, 0.478629, 0.236927, 0.236927].iter().map(|w| w/2.).collect(); PI - x.iter().zip(w) - .map(|(&x, w)| w*scattering_function_gl(x, impact_parameter, r0, relative_energy, interaction_potential) + .map(|(&x, w)| w*scattering_function_gl(x, impact_parameter, r0, relative_energy, interaction_potential.clone()) .with_context(|| format!("Numerical error: NaN in Gauss-Legendre scattering integral at x = {} with Er = {} eV and p = {} A.", x, relative_energy/EV, impact_parameter/ANGSTROM)) .unwrap()).sum::() } @@ -583,18 +579,6 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact let coefficients = interactions::polynomial_coefficients(relative_energy, impact_parameter, interaction_potential); let roots = real_polynomial_roots(coefficients.clone(), polynom_complex_threshold).unwrap(); - /* - println!("p={}", impact_parameter/ANGSTROM); - - for coefficient in &coefficients { - println!("{}", {coefficient}); - } - - for root in &roots { - println!("{} A", {root}); - } - */ - let max_root = roots.iter().cloned().fold(f64::NAN, f64::max); let inverse_transformed_root = interactions::inverse_transform(max_root, interaction_potential); @@ -642,8 +626,6 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame let g = |r: f64| -> f64 {interactions::distance_of_closest_approach_function_singularity_free(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)* interactions::scaling_function(r, impact_parameter, interaction_potential)}; - //Using upper bound const, ~10, construct upper bound as a plateau near 0 and linear increase away from that - //let upper_bound = f64::max(upper_bound_const*p, upper_bound_const*a); let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential); let roots = match derivative_free { @@ -683,14 +665,16 @@ pub fn newton_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par let f = |r: f64| -> f64 {interactions::distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)}; let df = |r: f64| -> f64 {interactions::diff_distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)}; - //Guess for large reduced energy from Mendenhall and Weller 1991 - //For small energies, use pure Newton-Raphson with arbitrary guess of 1 - let mut x0 = beta; let mut xn: f64; - if reduced_energy > 5. { + + //Guess for large reduced energy from Mendenhall and Weller 1991 + //For small energies, use pure Newton-Raphson with arbitrary guess of beta + let mut x0 = if reduced_energy > 5. { let inv_er_2 = 0.5/reduced_energy; - x0 = inv_er_2 + (inv_er_2*inv_er_2 + beta*beta).sqrt(); - } + inv_er_2 + (inv_er_2*inv_er_2 + beta*beta).sqrt() + } else { + beta + }; //Newton-Raphson to determine distance of closest approach let mut err: f64 = tolerance + 1.; diff --git a/src/geometry.rs b/src/geometry.rs index 616d9ad..624d36e 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -57,7 +57,7 @@ impl Geometry for Mesh0D { "NM" => NM, "M" => 1., _ => input.length_unit.parse() - .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", &input.length_unit.as_str())), }; @@ -295,7 +295,7 @@ impl Geometry for HomogeneousMesh2D { "NM" => NM, "M" => 1., _ => input.length_unit.parse() - .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", &input.length_unit.as_str())), }; @@ -353,16 +353,14 @@ impl Geometry for HomogeneousMesh2D { fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { if self.inside(x, y, z) { true - } else { - if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { - let distance = ((x - p.x()).powi(2) + (y - p.y()).powi(2)).sqrt(); - distance < self.energy_barrier_thickness - } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { - true - } else { - panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); - } - } + } else if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + let distance = ((x - p.x()).powi(2) + (y - p.y()).powi(2)).sqrt(); + distance < self.energy_barrier_thickness + } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + true + } else { + panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); + } } fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { @@ -401,7 +399,7 @@ impl Mesh2D { /// Finds the cell that is nearest to (x, y). fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { - let mut min_distance: f64 = std::f64::MAX; + let mut min_distance: f64 = f64::MAX; let mut index: usize = 0; for (cell_index, cell) in self.mesh.iter().enumerate() { @@ -760,7 +758,7 @@ impl Triangle2D { /// Calculates the shortest distance from this triangle to the point (x, y). pub fn distance_to(&self, x: f64, y: f64) -> f64 { - let mut distance_to = std::f64::MAX; + let mut distance_to = f64::MAX; for segment in &self.segments { let length_2 = (segment.2 - segment.0)*(segment.2 - segment.0) + (segment.3 - segment.1)*(segment.3 - segment.1); diff --git a/src/interactions.rs b/src/interactions.rs index 4687cf2..593e4e5 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -1,4 +1,5 @@ use super::*; +use std::sync::LazyLock; /// Analytic solutions to outermost root of the interaction potential. pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { @@ -9,7 +10,6 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { InteractionPotential::WW => 50.*ANGSTROM, _ => 10.*ANGSTROM, } - } fn smootherstep(x: f64, k: f64, x0: f64) -> f64 { @@ -53,7 +53,6 @@ pub fn interaction_potential(r: f64, a: f64, Za: f64, Zb: f64, interaction_poten InteractionPotential::FOUR_EIGHT{alpha, beta} => { four_eight(r, alpha, beta) } - } } @@ -271,17 +270,61 @@ pub fn dphi(xi: f64, interaction_potential: InteractionPotential) -> f64 { pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 { match interaction_potential { //ZBL screening length, Eckstein (4.1.8) - InteractionPotential::ZBL => 0.88534*A0/(Za.powf(0.23) + Zb.powf(0.23)), + InteractionPotential::ZBL => zbl_screening_length_lookup(Za as u64, Zb as u64), //Lindhard/Firsov screening length, Eckstein (4.1.5) - InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN | InteractionPotential::WW => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), - InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), + InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN | InteractionPotential::WW => lindhard_screening_length_lookup(Za as u64, Zb as u64), + InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64), InteractionPotential::MORSE{D, alpha, r0} => alpha, - InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => 0.88534*A0/(Z1.powf(0.23) + Z2.powf(0.23)), - InteractionPotential::KRC_MORSE{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), - InteractionPotential::FOUR_EIGHT{..} =>0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), + InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => zbl_screening_length_lookup(Za as u64, Zb as u64), + InteractionPotential::KRC_MORSE{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64), + InteractionPotential::FOUR_EIGHT{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64), } } +// It turns out it's faster (~10% speedup) to just generate every possible screening length as a lookup table +// LazyLock is a thread-safe value that is initialized whenever it is first accessed +// It will block other threads while it runs, but it should run extremely quickly and only once +const Z_MAX: usize = 120; +static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( + || + std::array::from_fn( + |i| { + // standard 2D to 1D array indexing + // I always have to look this up; copied from here + // (e.g. https://stackoverflow.com/questions/5494974/convert-1d-array-index-to-2d-array-index) + let Za = i / Z_MAX; + let Zb = i % Z_MAX; + lindhard_screening_length(Za as f64, Zb as f64) + } + ) +); +static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( + || + std::array::from_fn( + |i| { + let Za = i / Z_MAX; + let Zb = i % Z_MAX; + zbl_screening_length(Za as f64, Zb as f64) + } + ) +); + +pub fn zbl_screening_length(Za: f64, Zb: f64) -> f64{ + 0.88534*A0/(Za.powf(0.23) + Zb.powf(0.23)) +} + +pub fn lindhard_screening_length(Za: f64, Zb: f64) -> f64 { + 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.) +} + +pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { + LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] +} + +pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ + ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] +} + /// Coefficients of inverse-polynomial interaction potentials. /// Note: these go in order of decreasing orders of r, e.g., [r^n, r^n-1, r^n-2, ... r, 1]. pub fn polynomial_coefficients(relative_energy: f64, impact_parameter: f64, interaction_potential: InteractionPotential) -> Vec { @@ -414,11 +457,11 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 { } else if x <= x2 { - let a = [1.389653276380862E4, - -3.596912431628216E4, - 3.739206756369099E4, - -1.933748081656593E4, - 0.495516793802426E4, + let a = [1.389653276380862E4, + -3.596912431628216E4, + 3.739206756369099E4, + -1.933748081656593E4, + 0.495516793802426E4, -0.050264585985867E4]; (a[0] + a[1]*x + a[2]*x.powi(2) + a[3]*x.powi(3) + a[4]*x.powi(4) + a[5]*x.powi(5))*EV diff --git a/src/lib.rs b/src/lib.rs index 2c9c9d9..df3955f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -364,7 +364,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag } #[no_mangle] -pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut f64, uy: &mut f64, uz: &mut f64, E1: &mut f64, Z1: &mut f64, m1: &mut f64, Ec1: &mut f64, Es1: &mut f64, Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64) { +pub unsafe extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut f64, uy: &mut f64, uz: &mut f64, E1: &mut f64, Z1: &mut f64, m1: &mut f64, Ec1: &mut f64, Es1: &mut f64, Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64) { assert!(E1 > &mut 0.0); @@ -689,7 +689,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { } #[no_mangle] -pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track_recoils: &mut bool, +pub unsafe extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track_recoils: &mut bool, ux: *mut f64, uy: *mut f64, uz: *mut f64, E1: *mut f64, Z1: *mut f64, m1: *mut f64, Ec1: *mut f64, Es1: *mut f64, num_species_target: &mut c_int, @@ -1152,7 +1152,7 @@ pub fn reflect_single_ion_py(ion: &PyDict, target: &PyDict, vx: f64, vy: f64, vz let options = Options::default_options(false); - let velocity2 = vx.powf(2.) + vy.powf(2.) + vz.powf(2.); //m/s + let velocity2 = vx*vx + vy*vy + vz*vz; //m^2/s^2 let energy_eV = 0.5*m1*AMU*velocity2/EV; //EV let ux = vx/velocity2.sqrt(); diff --git a/src/main.rs b/src/main.rs index a3b0301..8b02d4b 100644 --- a/src/main.rs +++ b/src/main.rs @@ -41,7 +41,6 @@ use rand::{SeedableRng, rngs::ChaCha8Rng}; //Load internal modules pub mod material; pub mod particle; -pub mod tests; pub mod interactions; pub mod bca; pub mod geometry; diff --git a/src/material.rs b/src/material.rs index 7843e87..8904d4a 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,5 +1,6 @@ use super::*; use rand::RngExt; +use std::sync::LazyLock; ///This helper function is a workaround to issue #368 in serde fn default_surface_binding_model() -> SurfaceBindingModel { @@ -278,7 +279,7 @@ impl Material { let S_high = bethe_bloch_stopping_power_cross_section(Za, *Zb, E, Ma); - // correction applied only to LS component + // ck correction applied only to LS component 1./(1./S_high + 1./(S_low*ck)) }, //Lindhard-Scharff, Oen-Robinson, Lindhard Equipartition @@ -288,7 +289,7 @@ impl Material { ElectronicStoppingMode::INTERPOLATEDPLUS{ci} => { let S_high = bethe_bloch_stopping_power_cross_section(Za, *Zb, E, Ma); - // correction applied only to LS component + // ck correction applied only to LS component (S_high.powf(-ci) + (S_low*ck).powf(-ci)).powf(-1./ci) }, @@ -299,30 +300,56 @@ impl Material { } } +const Z_MAX: usize = 120; +// Generating lookup tables for all possibilities turns out to be faster than calculating on the fly +static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( + || + std::array::from_fn( + |i| { + let Za = i / Z_MAX; + let Zb = i % Z_MAX; + lindhard_scharff_stopping_power_constant(Za as f64, Zb as f64) + } + ) +); + +fn lindhard_scharff_stopping_power_constant(Za: f64, Zb: f64) -> f64 { + LINDHARD_SCHARFF_PREFACTOR*(Za*Za.cbrt().sqrt()*Zb)/(Za.cbrt().powi(2) + Zb.cbrt().powi(2)).powi(3).sqrt()*(AMU/Q).sqrt() +} + pub fn lindhard_scharff_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { - LINDHARD_SCHARFF_PREFACTOR*(Za*Za.cbrt().sqrt()*Zb)/(Za.cbrt().powi(2) + Zb.cbrt().powi(2)).powi(3).sqrt()*(E/Q/Ma*AMU).sqrt() + LS_STOPPING_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]*(E/Ma).sqrt() } +static BV_EMPIRICAL_MEAN_IONIZATON_POT: LazyLock<[f64; Z_MAX]> = LazyLock::new( + || + std::array::from_fn( + |Zb| { + let I0 = if (Zb as f64) < 13. { + 12. + 7./ Zb as f64 + } else { + 9.76 + 58.5*(Zb as f64).powf(-1.19) + }; + (Zb as f64)*I0*Q + } + ) +); + pub fn bethe_bloch_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { let beta = (1. - 1./(1. + E/Ma/C.powi(2)).powi(2)).sqrt(); let v = beta*C; - // This term is an empirical fit to the mean ionization potential - let I0 = match Zb < 13. { - true => 12. + 7./Zb, - false => 9.76 + 58.5*Zb.powf(-1.19), - }; - let I = Zb*I0*Q; + let I = BV_EMPIRICAL_MEAN_IONIZATON_POT[Zb as usize]; //See Biersack and Haggmark - this looks like an empirical shell correction - let B = match Zb < 3. { - true => 100.*Za/Zb, - false => 5. + let B = if Zb < 3. { + 100.*Za/Zb + } else { + 5. }; let prefactor = BETHE_BLOCH_PREFACTOR*Zb*Za*Za/beta/beta; let eb = 2.*ME*v*v/I; - prefactor*(eb + 1. + B/eb).ln() } diff --git a/src/math.rs b/src/math.rs index d8c5a59..3afc54a 100644 --- a/src/math.rs +++ b/src/math.rs @@ -4,17 +4,10 @@ use super::*; /// Duff et al., JCGT 2017 /// http://jcgt.org/published/0006/01/01/ pub fn duff_orthonormal_basis(n: Vector) -> (Vector, Vector) { - if n.z < 0. { - let a = 1.0 / (1.0 - n.z); - let b = n.x * n.y * a; - let b1 = Vector::new(1.0 - n.x*n.x*a, -b, n.x); - let b2 = Vector::new(b, n.y * n.y*a- 1.0, -n.y); - (b1, b2) - } else { - let a = 1.0 / (1.0 + n.z); - let b = -n.x * n.y * a; - let b1 = Vector::new(1.0 - n.x*n.x*a, b, -n.x); - let b2 = Vector::new(b, 1.0 - n.y*n.y*a, -n.y); - (b1, b2) - } + let sign = (1.0_f64).copysign(n.z); + let a = -1.0_f64 / (sign + n.z); + let b = n.x*n.y*a; + let b1 = Vector::new(1.0 + sign*n.x*n.x*a, sign*b, -sign*n.x); + let b2 = Vector::new(b, sign + n.y*n.y*a, -n.y); + (b1, b2) } \ No newline at end of file diff --git a/src/particle.rs b/src/particle.rs index cc1cf00..89039fa 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -1,7 +1,5 @@ use super::*; -/// Rustbca's internal representation of the particle_parameters input. - fn default_vec_zero() -> Vec { vec![0] } diff --git a/src/physics.rs b/src/physics.rs index c8d59a7..304fce7 100644 --- a/src/physics.rs +++ b/src/physics.rs @@ -75,5 +75,4 @@ pub fn physics_loop(particle_input_array: Vec