diff --git a/Cargo.toml b/Cargo.toml index f4738f1..a335ee2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "RustBCA" -version = "2.8.4" +version = "2.9.0" default-run = "RustBCA" authors = ["Jon Drobny ", "Jon Drobny "] edition = "2021" @@ -20,7 +20,7 @@ rand_distr = "0.4.3" toml = "0.7.4" anyhow = "1.0.71" itertools = "0.10.5" -rayon = "1.7.0" +rayon = "1.10.0" geo = {version = "0.25", optional = false} indicatif = {version = "0.15.0", features=["rayon"]} serde = { version = "1.0.163", features = ["derive"] } diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml index f6f64f1..106f8fd 100644 --- a/examples/boron_nitride_0D.toml +++ b/examples/boron_nitride_0D.toml @@ -1,3 +1,4 @@ +#Use with 0D geometry option only [options] name = "boron_nitride_" weak_collision_order = 0 diff --git a/examples/boron_nitride_sphere.toml b/examples/boron_nitride_sphere.toml index fd81322..d780ce2 100644 --- a/examples/boron_nitride_sphere.toml +++ b/examples/boron_nitride_sphere.toml @@ -1,3 +1,4 @@ +#Use with SPHERE geometry option only [options] name = "boron_nitride_" track_trajectories = true diff --git a/examples/boron_nitride_wire.toml b/examples/boron_nitride_wire.toml index 8f7aa84..d172890 100644 --- a/examples/boron_nitride_wire.toml +++ b/examples/boron_nitride_wire.toml @@ -1,3 +1,4 @@ +#Use with 2D geometry option only [options] name = "boron_nitride_" weak_collision_order = 0 diff --git a/examples/boron_nitride_wire_homogeneous.toml b/examples/boron_nitride_wire_homogeneous.toml index c7f40a1..57156b7 100644 --- a/examples/boron_nitride_wire_homogeneous.toml +++ b/examples/boron_nitride_wire_homogeneous.toml @@ -1,3 +1,4 @@ +#Use with 0D geometry option only [options] name = "boron_nitride_h_" weak_collision_order = 0 diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index e4f3cd8..235cc5a 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -1,3 +1,4 @@ +#Use with 2D geometry option only [options] name = "2000.0eV_0.0001deg_He_TiO2_Al_Si" track_trajectories = false diff --git a/examples/layered_geometry_1D.toml b/examples/layered_geometry_1D.toml index 27364ea..ae0a3e0 100644 --- a/examples/layered_geometry_1D.toml +++ b/examples/layered_geometry_1D.toml @@ -1,3 +1,4 @@ +#Use with 2D geometry option only [options] name = "2000.0eV_0.0001deg_He_TiO2_Al_Si" track_trajectories = false diff --git a/examples/lithium_vapor_shield.toml b/examples/lithium_vapor_shield.toml index e96f994..d458937 100644 --- a/examples/lithium_vapor_shield.toml +++ b/examples/lithium_vapor_shield.toml @@ -1,3 +1,4 @@ +#Use with 1D geometry option only [options] name = "lithium_vapor_shield_" track_trajectories = true diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index e5fee07..f2194ab 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -1,3 +1,4 @@ +#Use with 2D geometry input only [options] name = "2000.0eV_0.0001deg_He_H_TiO2_Al_Si" track_trajectories = true diff --git a/examples/test_morse.py b/examples/test_morse.py index 6de9119..227203c 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -180,8 +180,8 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): R_E_2 = np.zeros(num_energies) for index, energy in enumerate(energies): - R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=False) - R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=False) + R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=True) + R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=True) plt.semilogx(energies, R_N, label='R_N Morse-Kr-C H-Ni, Es=1.5eV', color='purple') plt.semilogx(energies, R_N_2, label='R_N Morse H-Ni, Es=1.5eV', color='green') diff --git a/examples/test_rustbca.py b/examples/test_rustbca.py index 3a8ca74..dd3d2b5 100644 --- a/examples/test_rustbca.py +++ b/examples/test_rustbca.py @@ -12,6 +12,7 @@ def main(): + #test rotation to and from RustBCA coordinates #nx, ny, nz is the normal vector (out of surface) @@ -35,6 +36,31 @@ def main(): assert(abs(uy) < 1e-6) assert(abs(uz) < 1e-6) + #test vectorized rotation to and from RustBCA coordinates + + num_rot_test = 1000000 + + #nx, ny, nz is the normal vector (out of surface) + nx = [-np.sqrt(2)/2]*num_rot_test + ny = [-np.sqrt(2)/2]*num_rot_test + nz = [0.0]*num_rot_test + + #ux, uy, uz is the particle direction (simulation coordinates) + ux = [1.0]*num_rot_test + uy = [0.0]*num_rot_test + uz = [0.0]*num_rot_test + + start = time.time() + ux, uy, uz = rotate_given_surface_normal_vec_py(nx, ny, nz, ux, uy, uz) + ux, uy, uz = rotate_back_vec_py(nx, ny, nz, ux, uy, uz) + stop = time.time() + print(f'Time to rotate: {(stop - start)/num_rot_test} sec/vector') + + #After rotating and rotating back, effect should be where you started (minus fp error) + assert(abs(ux[0] - 1.0) < 1e-6) + assert(abs(uy[0]) < 1e-6) + assert(abs(uz[0]) < 1e-6) + #scripts/materials.py has a number of potential ions and targets ion = helium ion['Eb'] = 0.0 @@ -128,7 +154,7 @@ def main(): print('Data processing complete.') print(f'RustBCA Y: {len(sputtered[:, 0])/number_ions} Yamamura Y: {yamamura_yield}') print(f'RustBCA R: {len(reflected[:, 0])/number_ions} Thomas R: {thomas}') - print(f'Time per ion: {delta_time/number_ions*1e3} us/{ion["symbol"]}') + print(f'Time per ion: {delta_time/number_ions} s/{ion["symbol"]}') #Next up is the layered target version. I'll add a 50 Angstrom layer of W-H to the top of the target. @@ -145,7 +171,7 @@ def main(): print(f'Running RustBCA for {number_ions} {ion["symbol"]} ions on {target["symbol"]} with hydrogenated layer at {energies_eV[0]/1000.} keV...') print(f'This may take several minutes.') - #Not the different argument order; when a breaking change is due, this will + #Note the different argument order; when a breaking change is due, this will #be back-ported to the other bindings as well for consistency. output, incident, stopped = compound_bca_list_1D_py( ux, uy, uz, energies_eV, [ion['Z']]*number_ions, @@ -153,6 +179,7 @@ def main(): [target['Ec'], 1.0], [target['Es'], 1.5], [target['Eb'], 0.0], [[target['n']/10**30, target['n']/10**30], [target['n']/10**30, 0.0]], [50.0, 1e6] ) + output = np.array(output) Z = output[:, 0] @@ -173,6 +200,80 @@ def main(): plt.plot([50.0, 50.0], [0.0, np.max(heights)*1.1]) plt.gca().set_ylim([0.0, np.max(heights)*1.1]) + number_ions = 10000 + + #1 keV is above the He on W sputtering threshold of ~150 eV + energies_eV = 1000.0*np.ones(number_ions) + + #Working with angles of exactly 0 is problematic due to gimbal lock + angle = 0.0001 + + #In RustBCA's 0D geometry, +x -> into the surface + ux = np.cos(angle*np.pi/180.)*np.ones(number_ions) + uy = np.sin(angle*np.pi/180.)*np.ones(number_ions) + uz = np.zeros(number_ions) + + Z1 = np.array([ion['Z']]*number_ions) + m1 = np.array([ion['m']]*number_ions) + Ec1 = np.array([ion['Ec']]*number_ions) + Es1 = np.array([ion['Es']]*number_ions) + + print(f'Running tracked RustBCA for {number_ions} {ion["symbol"]} ions on {target["symbol"]} with 1% {ion["symbol"]} at {energies_eV[0]/1000.} keV...') + print(f'This may take several minutes.') + + start = time.time() + # Note that simple_bca_list_py expects number densities in 1/Angstrom^3 + output, incident, incident_index = compound_bca_list_tracked_py( + energies_eV, ux, uy, uz, Z1, m1, Ec1, Es1, + [target['Z'], ion['Z']], [target['m'], ion['m']], + [target['Ec'], ion['Ec']], [target['Es'], ion['Es']], [target['n']/10**30, 0.01*target['n']/10**30], + [target['Eb'], 0.0]) + stop = time.time() + delta_time = stop - start + + output = np.array(output) + Z = output[:, 0] + m = output[:, 1] + E = output[:, 2] + x = output[:, 3] + y = output[:, 4] + z = output[:, 5] + ux = output[:, 6] + uy = output[:, 7] + uz = output[:, 8] + + #For the python bindings, these conditionals can be used to distinguish + #between sputtered, reflected, and implanted particles in the output list + sputtered = output[np.logical_and(Z == target['Z'], E > 0), :] + reflected = output[np.logical_and(Z == ion['Z'], x < 0), :] + implanted = output[np.logical_and(Z == ion['Z'], x > 0), :] + + plt.figure(1) + plt.title(f'Sputtered {target["symbol"]} Energy Distribution') + plt.xlabel('E [eV]') + plt.ylabel('f(E) [A.U.]') + plt.hist(sputtered[:, 2], bins=100, density=True, histtype='step') + + plt.figure(2) + plt.title(f'Implanted {ion["symbol"]} Depth Distribution') + plt.xlabel('x [A]') + plt.ylabel('f(x) [A.U.]') + plt.hist(implanted[:, 3], bins=100, density=True, histtype='step') + + thomas = thomas_reflection(ion, target, energies_eV[0]) + yamamura_yield = yamamura(ion, target, energies_eV[0]) + + print('Data processing complete.') + print(f'RustBCA Y: {len(sputtered[:, 0])/number_ions} Yamamura Y: {yamamura_yield}') + print(f'RustBCA R: {len(reflected[:, 0])/number_ions} Thomas R: {thomas}') + print(f'Time per ion: {delta_time/number_ions} s/{ion["symbol"]}') + + plt.figure() + plt.plot(incident_index) + plt.xlabel('Particle number') + plt.ylabel('Particle index') + plt.legend(['Incident', 'Indicies']) + plt.show() if __name__ == '__main__': diff --git a/examples/titanium_dioxide_0D.toml b/examples/titanium_dioxide_0D.toml index 70491c6..178f855 100644 --- a/examples/titanium_dioxide_0D.toml +++ b/examples/titanium_dioxide_0D.toml @@ -1,3 +1,4 @@ +#Use with 0D geometry option only [options] name = "10.0" track_trajectories = false diff --git a/examples/tungsten_tiles_trimesh.toml b/examples/tungsten_tiles_trimesh.toml index 43e5644..793590c 100644 --- a/examples/tungsten_tiles_trimesh.toml +++ b/examples/tungsten_tiles_trimesh.toml @@ -1,3 +1,4 @@ +#Use with TRIMESH geometry option only [options] name = "tungsten_tiles_" track_trajectories = true diff --git a/examples/tungsten_twist_trimesh.toml b/examples/tungsten_twist_trimesh.toml index 01051d9..0664f3e 100644 --- a/examples/tungsten_twist_trimesh.toml +++ b/examples/tungsten_twist_trimesh.toml @@ -1,3 +1,4 @@ +#Use with TRIMESH geometry option only [options] name = "tungsten_twist_" track_trajectories = true diff --git a/src/geometry.rs b/src/geometry.rs index 8e3d31a..a915629 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -66,6 +66,7 @@ impl Geometry for Mesh0D { let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); + assert!(densities.len() > 0, "Input Error: density list empty."); let total_density: f64 = densities.iter().sum(); @@ -135,6 +136,7 @@ impl Geometry for Mesh1D { let layer_thicknesses = geometry_input.layer_thicknesses.clone(); let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone(); + assert!(electronic_stopping_correction_factors.len() > 0, "Input Error: Electronic stopping correction factor list empty."); let n = layer_thicknesses.len(); let mut layers: Vec = Vec::with_capacity(n); @@ -433,7 +435,7 @@ impl Geometry for Mesh2D { let simulation_boundary_points = geometry_input.simulation_boundary_points.clone(); let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone(); - + assert!(electronic_stopping_correction_factors.len() > 0, "Input Error: Electronic stopping correction factor list empty."); let n = triangles.len(); let mut cells: Vec = Vec::with_capacity(n); diff --git a/src/lib.rs b/src/lib.rs index 05817c3..3a9f53e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -39,7 +39,7 @@ use std::slice; use std::sync::Mutex; //itertools -use itertools::izip; +use itertools::{izip}; //Math use std::f64::consts::FRAC_2_SQRT_PI; @@ -89,6 +89,7 @@ pub fn libRustBCA(py: Python, m: &PyModule) -> PyResult<()> { m.add_function(wrap_pyfunction!(simple_bca_py, m)?)?; m.add_function(wrap_pyfunction!(simple_bca_list_py, m)?)?; m.add_function(wrap_pyfunction!(compound_bca_list_py, m)?)?; + m.add_function(wrap_pyfunction!(compound_bca_list_tracked_py, m)?)?; m.add_function(wrap_pyfunction!(compound_bca_list_1D_py, m)?)?; m.add_function(wrap_pyfunction!(sputtering_yield, m)?)?; m.add_function(wrap_pyfunction!(reflection_coefficient, m)?)?; @@ -97,7 +98,11 @@ pub fn libRustBCA(py: Python, m: &PyModule) -> PyResult<()> { #[cfg(feature = "parry3d")] m.add_function(wrap_pyfunction!(rotate_given_surface_normal_py, m)?)?; #[cfg(feature = "parry3d")] + m.add_function(wrap_pyfunction!(rotate_given_surface_normal_vec_py, m)?)?; + #[cfg(feature = "parry3d")] m.add_function(wrap_pyfunction!(rotate_back_py, m)?)?; + #[cfg(feature = "parry3d")] + m.add_function(wrap_pyfunction!(rotate_back_vec_py, m)?)?; Ok(()) } @@ -794,7 +799,7 @@ pub extern "C" fn simple_bca_c(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64 } #[cfg(feature = "python")] -///compound_tagged_bca_list_py(ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2) +///compound_\\\\\\_bca_list_py(ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2) /// runs a BCA simulation for a list of particles and outputs a list of sputtered, reflected, and implanted particles. /// Args: /// energies (list(f64)): initial ion energies in eV. @@ -914,6 +919,136 @@ pub fn compound_bca_list_py(energies: Vec, ux: Vec, uy: Vec, uz: (total_output, incident) } +#[cfg(feature = "python")] +///compound_bca_list_tracked_py(ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2) +/// runs a BCA simulation for a list of particles and outputs a list of sputtered, reflected, and implanted particles. +/// Args: +/// energies (list(f64)): initial ion energies in eV. +/// ux (list(f64)): initial ion directions x. ux != 0.0 to avoid gimbal lock +/// uy (list(f64)): initial ion directions y. +/// uz (list(f64)): initial ion directions z. +/// Z1 (list(f64)): initial ion atomic numbers. +/// m1 (list(f64)): initial ion masses in amu. +/// Ec1 (list(f64)): ion cutoff energies in eV. If ion energy < Ec1, it stops in the material. +/// Es1 (list(f64)): ion surface binding energies. Assumed planar. +/// Z2 (list(f64)): target material species atomic numbers. +/// m2 (list(f64)): target material species masses in amu. +/// Ec2 (list(f64)): target material species cutoff energies in eV. If recoil energy < Ec2, it stops in the material. +/// Es2 (list(f64)): target species surface binding energies. Assumed planar. +/// n2 (list(f64)): target material species atomic number densities in inverse cubic Angstroms. +/// Eb2 (list(f64)): target material species bulk binding energies in eV. +/// Returns: +/// output (NX9 list of f64): each row in the list represents an output particle (implanted, +/// sputtered, or reflected). Each row consists of: +/// [Z, m (amu), E (eV), x, y, z, (angstrom), ux, uy, uz] +/// incident (list(bool)): whether each row of output was an incident ion or originated in the target +/// incident_index (list(usize)): index of incident particle that caused this particle to be emitted +#[pyfunction] +pub fn compound_bca_list_tracked_py(energies: Vec, ux: Vec, uy: Vec, uz: Vec, Z1: Vec, m1: Vec, Ec1: Vec, Es1: Vec, Z2: Vec, m2: Vec, Ec2: Vec, Es2: Vec, n2: Vec, Eb2: Vec) -> (Vec<[f64; 9]>, Vec, Vec) { + let mut total_output = vec![]; + let mut incident = vec![]; + let mut incident_index = vec![]; + let num_species_target = Z2.len(); + let num_incident_ions = energies.len(); + + assert_eq!(ux.len(), num_incident_ions, "Input error: list of x-directions is not the same length as list of incident energies."); + assert_eq!(uy.len(), num_incident_ions, "Input error: list of y-directions is not the same length as list of incident energies."); + assert_eq!(uz.len(), num_incident_ions, "Input error: list of z-directions is not the same length as list of incident energies."); + assert_eq!(Z1.len(), num_incident_ions, "Input error: list of incident atomic numbers is not the same length as list of incident energies."); + assert_eq!(m1.len(), num_incident_ions, "Input error: list of incident atomic masses is not the same length as list of incident energies."); + assert_eq!(Es1.len(), num_incident_ions, "Input error: list of incident surface binding energies is not the same length as list of incident energies."); + assert_eq!(Ec1.len(), num_incident_ions, "Input error: list of incident cutoff energies is not the same length as list of incident energies."); + + assert_eq!(m2.len(), num_species_target, "Input error: list of target atomic masses is not the same length as atomic numbers."); + assert_eq!(Ec2.len(), num_species_target, "Input error: list of target cutoff energies is not the same length as atomic numbers."); + assert_eq!(Es2.len(), num_species_target, "Input error: list of target surface binding energies is not the same length as atomic numbers."); + assert_eq!(Eb2.len(), num_species_target, "Input error: list of target bulk binding energies is not the same length as atomic numbers."); + assert_eq!(n2.len(), num_species_target, "Input error: list of target number densities is not the same length as atomic numbers."); + + let options = Options::default_options(true); + //options.high_energy_free_flight_paths = true; + + let x = -2.*(n2.iter().sum::()*10E30).powf(-1./3.); + let y = 0.0; + let z = 0.0; + + let material_parameters = material::MaterialParameters { + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: Eb2, + Es: Es2, + Ec: Ec2, + Ed: vec![0.0; num_species_target], + Z: Z2, + m: m2, + interaction_index: vec![0; num_species_target], + surface_binding_model: SurfaceBindingModel::INDIVIDUAL, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let geometry_input = geometry::Mesh0DInput { + length_unit: "ANGSTROM".to_string(), + densities: n2, + electronic_stopping_correction_factor: 1.0 + }; + + let m = material::Material::::new(&material_parameters, &geometry_input); + + let mut finished_particles: Vec = Vec::new(); + + let incident_particles: Vec = izip!(energies, ux, uy, uz, Z1, Ec1, Es1, m1) + .enumerate() + .map(|(index, (energy, ux_, uy_, uz_, Z1_, Ec1_, Es1_, m1_))| { + let mut p = particle::Particle::default_incident( + m1_, + Z1_, + energy, + Ec1_, + Es1_, + x, + ux_, + uy_, + uz_ + ); + p.tag = index as i32; + p + }).collect(); + + finished_particles.par_extend( + incident_particles.into_par_iter() + .map(|particle| bca::single_ion_bca(particle, &m, &options)) + .flatten() + ); + + for particle in finished_particles { + if (particle.left) | (particle.incident) { + incident.push(particle.incident); + incident_index.push(particle.tag as usize); + let mut energy_out; + if particle.stopped { + energy_out = 0.; + } else { + energy_out = particle.E/EV; + } + total_output.push( + [ + particle.Z, + particle.m/AMU, + energy_out, + particle.pos.x/ANGSTROM, + particle.pos.y/ANGSTROM, + particle.pos.z/ANGSTROM, + particle.dir.x, + particle.dir.y, + particle.dir.z, + ] + ) + } + } + + (total_output, incident, incident_index) +} + #[cfg(feature = "python")] ///reflect_single_ion_py(ion, target, vx, vy, vz) ///Performs a single BCA ion trajectory in target material with specified incident velocity. @@ -1382,7 +1517,8 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); let c = into_surface.dot(&RUSTBCA_DIRECTION); let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let rotation_matrix = if c != -1.0 { + + let rotation_matrix = if (c + 1.0).abs() > 1e-6 { Matrix3::identity() + vx + vx*vx/(1. + c) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation @@ -1394,12 +1530,12 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu // ux must not be exactly 1.0 to avoid gimbal lock in RustBCA // simple_bca does not normalize direction before proceeding, must be done manually - *ux = incident.x + DELTA; assert!( - *ux > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = ({}, {}, {}) u = ({}, {}, {})", - nx, ny, nz, ux, uy, uz + incident.x + DELTA > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. c={} n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", + (c + 1.0).abs(), nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z ); + *ux = incident.x + DELTA; *uy = incident.y - DELTA; *uz = incident.z; let mag = (ux.powf(2.) + uy.powf(2.) + uz.powf(2.)).sqrt(); @@ -1409,6 +1545,8 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu *uz /= mag; } + + #[cfg(all(feature = "python", feature = "parry3d"))] #[pyfunction] /// rotate_given_surface_normal_py(nx, ny, nz, ux, uy, uz) @@ -1432,6 +1570,46 @@ pub fn rotate_given_surface_normal_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f6 (ux, uy, uz) } +#[cfg(all(feature = "python", feature = "parry3d"))] +#[pyfunction] +/// rotate_given_surface_normal_vec_py(nx, ny, nz, ux, uy, uz) +/// -- +/// +/// This function takes a particle direction and a normal vector and rotates from simulation to RustBCA coordinates. +/// Args: +/// nx (list(f64)): surface normal in global frame x-component. +/// ny (list(f64)): surface normal in global frame y-component. +/// nz (list(f64)): surface normal in global frame z-component. +/// ux (list(f64)): particle direction in global frame x-component. +/// uy (list(f64)): particle direction in global frame normal y-component. +/// uz (list(f64)): particle direction in global frame normal z-component. +/// Returns: +/// direction (list(f64), list(f64), list(f64)): direction vector of particle in RustBCA coordinates. +/// Note: non-incident particles will be returned with ux, uy, uz = (0, 0, 0) +pub fn rotate_given_surface_normal_vec_py(nx: Vec, ny: Vec, nz: Vec, ux: Vec, uy: Vec, uz: Vec) -> (Vec, Vec, Vec) { + + let length = nx.len(); + + let mut ux_new = Vec::with_capacity(length); + let mut uy_new = Vec::with_capacity(length); + let mut uz_new = Vec::with_capacity(length); + + (0..length).into_iter().for_each(|index| { + + let mut ux_ = ux[index]; + let mut uy_ = uy[index]; + let mut uz_ = uz[index]; + + rotate_given_surface_normal(nx[index], ny[index], nz[index], &mut ux_, &mut uy_, &mut uz_); + ux_new.push(ux_); + uy_new.push(uy_); + uz_new.push(uz_); + + }); + + (ux_new, uy_new, uz_new) +} + #[cfg(feature = "parry3d")] #[no_mangle] pub extern "C" fn rotate_back(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut f64, uz: &mut f64) { @@ -1487,6 +1665,36 @@ pub fn rotate_back_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f64, uz: f64) -> ( (ux, uy, uz) } +#[cfg(all(feature = "python", feature = "parry3d"))] +#[pyfunction] +/// rotate_back_vec_py(nx, ny, nz, ux, uy, uz) +/// -- +/// +/// This function takes a RustBCA particle direction and a normal vector and rotates back from RustBCA to simulation coordinates. +/// Args: +/// nx (list(f64)): surface normal in global frame x-component. +/// ny (list(f64)): surface normal in global frame y-component. +/// nz (list(f64)): surface normal in global frame z-component. +/// ux (list(f64)): particle direction in global frame x-component. +/// uy (list(f64)): particle direction in global frame normal y-component. +/// uz (list(f64)): particle direction in global frame normal z-component. +/// Returns: +/// direction (list(f64), list(f64), list(f64)): direction vector of particle in simulation coordinates. +pub fn rotate_back_vec_py(nx: Vec, ny: Vec, nz: Vec, ux: Vec, uy: Vec, uz: Vec) -> (Vec, Vec, Vec) { + + let (ux_new, (uy_new, uz_new)) = (nx, ny, nz, ux, uy, uz).into_par_iter().map(|(nx_, ny_, nz_, ux_, uy_, uz_)| { + + let mut ux_mut = ux_; + let mut uy_mut = uy_; + let mut uz_mut = uz_; + rotate_back(nx_, ny_, nz_, &mut ux_mut, &mut uy_mut, &mut uz_mut); + + (ux_mut, (uy_mut, uz_mut)) + }).unzip(); + + (ux_new, uy_new, uz_new) +} + #[cfg(feature = "python")] /// A helper function to unpack a python float from a python any. fn unpack(python_float: &PyAny) -> f64 {