Skip to content

Powf to powi #307

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
drobnyjt wants to merge 5 commits into
base: main
Choose a base branch
from
Open

Powf to powi #307

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c17d50d
Fix Biersack-Varelas stopping power.
drobnyjt Jun 23, 2026
8c41b67
Changed all invocations of powf(<integer>) to appropriate powi, sqrt,…
Jul 21, 2026
895d510
Fixed int->float
Jul 21, 2026
83acc2a
Fix to beta
Jul 21, 2026
97e32ce
Update rustbca_compile_check.yml
drobnyjt Jul 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .github/workflows/rustbca_compile_check.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
name: RustBCA Compile check

on:
workflow_dispatch:
push:
branches: [ main, dev ]
pull_request:
Expand Down
4 changes: 2 additions & 2 deletions src/bca.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -538,7 +538,7 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: &
/// Mendenhall-Weller scattering integrand.
fn scattering_integral_mw(x: f64, beta: f64, reduced_energy: f64, interaction_potential: InteractionPotential) -> f64 {
//Function for scattering integral - see Mendenhall and Weller, 1991 & 2005
return (1. - interactions::phi(x, interaction_potential)/x/reduced_energy - beta*beta/x/x).powf(-0.5);
return 1./(1. - interactions::phi(x, interaction_potential)/x/reduced_energy - beta*beta/x/x).sqrt();
}

/// Gauss-Legendre scattering integrand.
Expand Down Expand Up @@ -776,7 +776,7 @@ pub fn mendenhall_weller(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par
let beta: f64 = impact_parameter/a;

//Scattering integral quadrature from Mendenhall and Weller 2005
let lambda0 = (0.5 + beta*beta/x0/x0/2. - interactions::dphi(x0, interaction_potential)/2./reduced_energy).powf(-1./2.);
let lambda0 = 1./(0.5 + beta*beta/x0/x0/2. - interactions::dphi(x0, interaction_potential)/2./reduced_energy).sqrt();
let alpha = 1./12.*(1. + lambda0 + 5.*(0.4206*scattering_integral_mw(x0/0.9072, beta, reduced_energy, interaction_potential) + 0.9072*scattering_integral_mw(x0/0.4206, beta, reduced_energy, interaction_potential)));
PI*(1. - beta*alpha/x0)
}
Expand Down
10 changes: 5 additions & 5 deletions src/geometry.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ impl Geometry for Mesh0D {

let total_density: f64 = densities.iter().sum();

let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;

let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();

Expand Down Expand Up @@ -180,8 +180,8 @@ impl Geometry for Mesh1D {
layer_top = layer_bottom;
}

let top_energy_barrier_thickness = layers[0].densities.iter().sum::<f64>().powf(-1./3.)/SQRTPI*2.;
let bottom_energy_barrier_thickness = layers[layers.len() - 1].densities.iter().sum::<f64>().powf(-1./3.)/SQRTPI*2.;
let top_energy_barrier_thickness = 1./layers[0].densities.iter().sum::<f64>().cbrt()/SQRTPI*2.;
let bottom_energy_barrier_thickness = 1./layers[layers.len() - 1].densities.iter().sum::<f64>().cbrt()/SQRTPI*2.;

Mesh1D {
layers,
Expand Down Expand Up @@ -315,7 +315,7 @@ impl Geometry for HomogeneousMesh2D {

let total_density: f64 = densities.iter().sum();

let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;

let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();

Expand Down Expand Up @@ -361,7 +361,7 @@ impl Geometry for HomogeneousMesh2D {
true
} else {
if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
let distance = ((x - p.x()).powf(2.) + (y - p.y()).powf(2.)).sqrt();
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
Expand Down
26 changes: 11 additions & 15 deletions src/interactions.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,15 +166,12 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent
1./(1. + (r/a).powi(2))
},
InteractionPotential::LENNARD_JONES_12_6{sigma, ..} => {
let n = 11.;
1./(1. + (r/sigma).powf(n))
1./(1. + (r/sigma).powi(11))
},
InteractionPotential::LENNARD_JONES_65_6{sigma, ..} => {
let n = 6.;
1./(1. + (r/sigma).powf(n))
1./(1. + (r/sigma).powi(6))
},
InteractionPotential::FOUR_EIGHT{alpha, beta} => {
let n = 8.;
1./(1. + r.powi(8)/beta)
}
InteractionPotential::MORSE{D, alpha, r0} => {
Expand Down Expand Up @@ -294,15 +291,14 @@ pub fn polynomial_coefficients(relative_energy: f64, impact_parameter: f64, inte
let epsilon_ev = epsilon/EV;
let sigma_angstroms = sigma/ANGSTROM;
let relative_energy_ev = relative_energy/EV;
//vec![1., 0., -impact_parameter.powi(2), 0., 0., 0., 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, 0., 0., 0., 0., 0., -4.*epsilon_ev*sigma.powf(12.)/relative_energy_ev]
vec![1.0, -impact_parameter.powi(2), 0.0, 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, 0.0, 0.0, -4.*epsilon_ev*sigma.powf(12.)/relative_energy_ev]
vec![1.0, -impact_parameter.powi(2), 0.0, 4.*epsilon_ev*sigma.powi(6)/relative_energy_ev, 0.0, 0.0, -4.*epsilon_ev*sigma.powi(12)/relative_energy_ev]
},
InteractionPotential::LENNARD_JONES_65_6{sigma, epsilon} => {
let impact_parameter_angstroms = impact_parameter/ANGSTROM;
let epsilon_ev = epsilon/EV;
let sigma_angstroms = sigma/ANGSTROM;
let relative_energy_ev = relative_energy/EV;
vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, -4.*epsilon_ev*sigma.powf(6.5)/relative_energy_ev]
vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon_ev*sigma.powi(6)/relative_energy_ev, -4.*epsilon_ev*sigma.powf(6.5)/relative_energy_ev]
},
InteractionPotential::FOUR_EIGHT{alpha, beta} => {
//Note: I've transformed to angstroms here to help the rootfinder with numerical issues.
Expand Down Expand Up @@ -343,12 +339,12 @@ pub fn four_eight(r: f64, alpha: f64, beta: f64) -> f64 {

/// Lennard-Jones 12-6
pub fn lennard_jones(r: f64, sigma: f64, epsilon: f64) -> f64 {
4.*epsilon*((sigma/r).powf(12.) - (sigma/r).powf(6.))
4.*epsilon*((sigma/r).powi(12) - (sigma/r).powi(6))
}

/// Lennard-Jones 6.5-6
pub fn lennard_jones_65_6(r: f64, sigma: f64, epsilon: f64) -> f64 {
4.*epsilon*((sigma/r).powf(6.5) - (sigma/r).powf(6.))
4.*epsilon*((sigma/r).powf(6.5) - (sigma/r).powi(6))
}

/// Morse potential
Expand All @@ -366,7 +362,7 @@ pub fn doca_four_eight(r: f64, impact_parameter: f64, relative_energy: f64, alph
let a = alpha.powf(1./4.);
let b = beta.powf(1./8.);
let b4 = beta.sqrt();
(r/b).powf(8.) - (-(a*r/b/b) + 1.)/relative_energy - (impact_parameter*r.powf(3.)/b4)
(r/b).powi(8) - (-(a*r/b/b) + 1.)/relative_energy - (impact_parameter*r.powi(3)/b4)
}

/// Distance of closest approach function for Morse potential.
Expand All @@ -391,17 +387,17 @@ pub fn doca_lennard_jones_65_6(r: f64, p: f64, relative_energy: f64, sigma: f64,

/// Distance of closest approach function for LJ 12-6 potential.
pub fn doca_lennard_jones(r: f64, p: f64, relative_energy: f64, sigma: f64, epsilon: f64) -> f64 {
(r/sigma).powf(12.) - 4.*epsilon/relative_energy*(1. - (r/sigma).powf(6.)) - p.powi(2)*r.powf(10.)/sigma.powf(12.)
(r/sigma).powi(12) - 4.*epsilon/relative_energy*(1. - (r/sigma).powi(6)) - p.powi(2)*r.powi(10)/sigma.powi(12)
}

/// First derivative w.r.t. `r` of the distance of closest approach function for LJ 12-6 potential.
pub fn diff_doca_lennard_jones(r: f64, p: f64, relative_energy: f64, sigma: f64, epsilon: f64) -> f64 {
12.*(r/sigma).powf(11.)/sigma + 4.*epsilon/relative_energy*6.*(r/sigma).powf(5.)/sigma - 10.*p.powi(2)*r.powf(9.)/sigma.powf(12.)
12.*(r/sigma).powi(11)/sigma + 4.*epsilon/relative_energy*6.*(r/sigma).powi(5)/sigma - 10.*p.powi(2)*r.powi(9)/sigma.powi(12)
}

/// First derivative w.r.t. `r` of the distance of closest approach function for LJ 6.5-6 potential.
pub fn diff_doca_lennard_jones_65_6(r: f64, p: f64, relative_energy: f64, sigma: f64, epsilon: f64) -> f64 {
6.5*(r/sigma).powf(5.5)/sigma + 4.*epsilon/relative_energy*0.5*(sigma*r).powf(-0.5) - (p/sigma).powi(2)*4.5*(r/sigma).powf(3.5)/sigma
6.5*(r/sigma).powf(5.5)/sigma + 4.*epsilon/relative_energy*0.5/(sigma*r).sqrt() - (p/sigma).powi(2)*4.5*(r/sigma).powf(3.5)/sigma
}

/// W-W cublic spline potential from Bjorkas et al.
Expand All @@ -427,7 +423,7 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 {
-0.050264585985867E4
];

(a[0] + a[1]*x + a[2]*x.powi(2) + a[3]*x.powi(3) + a[4]*x.powf(4.) + a[5]*x.powf(5.))*EV
(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

} else {

Expand Down
18 changes: 9 additions & 9 deletions src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -235,7 +235,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag
let tags = unsafe { slice::from_raw_parts(input.tags, input.len).to_vec() };
let weights = unsafe { slice::from_raw_parts(input.weights, input.len).to_vec() };

let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
let x = -2.*(n2.iter().sum::<f64>()*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -368,7 +368,7 @@ pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut
let Es2 = unsafe { slice::from_raw_parts(Es2, *num_species_target as usize).to_vec() };
let Eb2 = unsafe { slice::from_raw_parts(Eb2, *num_species_target as usize).to_vec() };

let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
let x = -2.*(n2.iter().sum::<f64>()*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -438,7 +438,7 @@ pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut
#[no_mangle]
pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA {

let x = -2.*(input.n2*10E30).powf(-1./3.);
let x = -2.*(input.n2*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -560,7 +560,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA {
let Es2 = unsafe { slice::from_raw_parts(input.Es2, input.num_species_target).to_vec() };
let Eb2 = unsafe { slice::from_raw_parts(input.Eb2, input.num_species_target).to_vec() };

let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
let x = -2.*(n2.iter().sum::<f64>()*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -699,7 +699,7 @@ pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track

//println!("Z2: {} m2: {} n2: {} Ec2: {} Es2: {} Eb2: {}", Z2[0], m2[0], n2[0], Ec2[0], Es2[0], Eb2[0]);

let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
let x = -2.*(n2.iter().sum::<f64>()*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -844,7 +844,7 @@ pub fn compound_bca_list_py(energies: Vec<f64>, ux: Vec<f64>, uy: Vec<f64>, uz:

let options = Options::default_options(true);

let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
let x = -2.*(n2.iter().sum::<f64>()*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -968,7 +968,7 @@ pub fn compound_bca_list_tracked_py(energies: Vec<f64>, ux: Vec<f64>, uy: Vec<f6
let options = Options::default_options(true);
//options.high_energy_free_flight_paths = true;

let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
let x = -2.*(n2.iter().sum::<f64>()*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -1333,7 +1333,7 @@ pub fn simple_bca_list_py(energies: Vec<f64>, usx: Vec<f64>, usy: Vec<f64>, usz:
assert_eq!(energies.len(), usy.len());
assert_eq!(energies.len(), usz.len());

let x = -2.*(n2*10E30).powf(-1./3.);
let x = -2.*(n2*1E30).powf(-1./3.);
let y = 0.0;
let z = 0.0;

Expand Down Expand Up @@ -1525,7 +1525,7 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu
*ux = incident.x;
*uy = incident.y;
*uz = incident.z;
let mag = (ux.powf(2.) + uy.powf(2.) + uz.powf(2.)).sqrt();
let mag = (ux*ux + uy*uy + uz*uz).sqrt();

*ux /= mag;
*uy /= mag;
Expand Down
9 changes: 5 additions & 4 deletions src/material.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ impl <T: Geometry> Material<T> {

/// Determines the local mean free path from the formula sum(n(x, y))^(-1/3)
pub fn mfp(&self, x: f64, y: f64, z: f64) -> f64 {
return self.total_number_density(x, y, z).powf(-1./3.);
return 1./self.total_number_density(x, y, z).cbrt();
}

/// Total number density of triangle that contains or is nearest to (x, y).
Expand Down Expand Up @@ -274,7 +274,7 @@ impl <T: Geometry> Material<T> {

for (n, Zb) in self.number_densities(x, y, z).iter().zip(&self.Z) {

let beta = (1. - (1. + E/Ma/C.powi(2)).powf(-2.)).sqrt();
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
Expand All @@ -296,11 +296,12 @@ impl <T: Geometry> Material<T> {
let S_high = prefactor*(eb + 1. + B/eb).ln();

//Lindhard-Scharff electronic stopping
let S_low = LINDHARD_SCHARFF_PREFACTOR*(Za.powf(7./6.)*Zb)/(Za.powf(2./3.) + Zb.powf(2./3.)).powf(3./2.)*(E/Q/Ma*AMU).sqrt();
//let S_low = LINDHARD_SCHARFF_PREFACTOR*(Za.powf(7./6.)*Zb)/(Za.powf(2./3.) + Zb.powf(2./3.)).powf(3./2.)*(E/Q/Ma*AMU).sqrt();
let S_low = LINDHARD_SCHARFF_PREFACTOR*(Za*Za.cbrt().sqrt()*Zb)/(Za.cbrt().powi(2) + Zb.cbrt().powi(2)).powi(3).sqrt()*(E/Q/Ma*AMU).sqrt();

let stopping_power = match electronic_stopping_mode {
//Biersack-Varelas Interpolation
ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low*ck),
ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low)*ck,
//Oen-Robinson
ElectronicStoppingMode::LOW_ENERGY_LOCAL => S_low*ck,
//Lindhard-Scharff
Expand Down
Loading