Move particle advance and rotate into Particle.

Author: drobnyjt

src/bca.rs

@@ -134,10 +134,10 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
                     total_asymptotic_deflection += binary_collision_result.asymptotic_deflection;
 
                     //Rotate particle 1, 2 by lab frame scattering angles
-                    particle::rotate_particle(&mut particle_1, binary_collision_result.psi,
+                    particle_1.rotate(binary_collision_result.psi,
                         binary_collision_geometry.phi_azimuthal);
 
-                    particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil,
+                    particle_2.rotate(-binary_collision_result.psi_recoil,
                         binary_collision_geometry.phi_azimuthal);
 
                     particle_2.dir_old.x = particle_2.dir.x;
@@ -184,11 +184,11 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
 
             //Advance particle in space and track total distance traveled
             #[cfg(not(feature = "accelerated_ions"))]
-            let distance_traveled = particle::particle_advance(&mut particle_1,
+            let distance_traveled = particle_1.advance(
                 binary_collision_geometries[0].mfp, total_asymptotic_deflection);
 
             #[cfg(feature = "accelerated_ions")]
-            let distance_traveled = particle::particle_advance(&mut particle_1,
+            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

src/lib.rs

@@ -984,9 +984,46 @@ 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)
 /// 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
 #[pyfunction]
-pub fn compound_bca_list_py(num_species_target: usize, energies: Vec<f64>, ux: Vec<f64>, uy: Vec<f64>, uz: Vec<f64>, Z1: Vec<f64>, m1: Vec<f64>, Ec1: Vec<f64>, Es1: Vec<f64>, Z2: Vec<f64>, m2: Vec<f64>, Ec2: Vec<f64>, Es2: Vec<f64>, n2: Vec<f64>, Eb2: Vec<f64>) -> Vec<[f64; 6]> {
+pub fn compound_bca_list_py(energies: Vec<f64>, ux: Vec<f64>, uy: Vec<f64>, uz: Vec<f64>, Z1: Vec<f64>, m1: Vec<f64>, Ec1: Vec<f64>, Es1: Vec<f64>, Z2: Vec<f64>, m2: Vec<f64>, Ec2: Vec<f64>, Es2: Vec<f64>, n2: Vec<f64>, Eb2: Vec<f64>) -> (Vec<[f64; 9]>, Vec<bool>) {
     let mut total_output = vec![];
+    let mut incident = 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.");
 
     #[cfg(feature = "distributions")]
     let options = Options {
@@ -1107,8 +1144,10 @@ pub fn compound_bca_list_py(num_species_target: usize, energies: Vec<f64>, ux: V
         let output = bca::single_ion_bca(p, &m, &options);
 
         for particle in output {
-
             if (particle.left) | (particle.incident) {
+
+                incident.push(particle.incident);
+
                 if particle.stopped {
                     energy_out = 0.
                 } else {
@@ -1119,6 +1158,9 @@ pub fn compound_bca_list_py(num_species_target: usize, energies: Vec<f64>, ux: V
                         particle.Z,
                         particle.m/AMU,
                         energy_out,
+                        particle.pos.x,
+                        particle.pos.y,
+                        particle.pos.z,
                         particle.dir.x,
                         particle.dir.y,
                         particle.dir.z,
@@ -1128,7 +1170,7 @@ pub fn compound_bca_list_py(num_species_target: usize, energies: Vec<f64>, ux: V
         }
         index += 1;
     }
-    total_output
+    (total_output, incident)
 }
 
 #[cfg(feature = "python")]

src/particle.rs

@@ -83,7 +83,7 @@ pub struct Particle {
     pub left: bool,
     pub incident: bool,
     pub first_step: bool,
-    pub trajectory: Vec<Vector4>,
+    pub trajectory: Vec<TrajectoryElement>,
     pub energies: Vec<EnergyLoss>,
     pub track_trajectories: bool,
     pub number_collision_events: usize,
@@ -122,7 +122,7 @@ impl Particle {
             left: false,
             incident: true,
             first_step: true,
-            trajectory: vec![Vector4::new(input.E, input.x, input.y, input.z)],
+            trajectory: vec![TrajectoryElement::new(input.E, input.x, input.y, input.z)],
             energies: vec![],
             track_trajectories: options.track_trajectories,
             number_collision_events: 0,
@@ -170,7 +170,7 @@ impl Particle {
     /// If `track_trajectories`, add the current (E, x, y, z) to the trajectory.
     pub fn add_trajectory(&mut self) {
         if self.track_trajectories {
-            self.trajectory.push(Vector4 {E: self.E, x: self.pos.x, y: self.pos.y, z: self.pos.z});
+            self.trajectory.push(TrajectoryElement {E: self.E, x: self.pos.x, y: self.pos.y, z: self.pos.z});
         }
     }
 
@@ -190,60 +190,62 @@ impl Particle {
             self.m*speed*self.dir.z,
         )
     }
-}
 
-/// Rotate a particle by a deflection psi at an azimuthal angle phi.
-pub fn rotate_particle(particle_1: &mut particle::Particle, psi: f64, phi: f64) {
-    let cosx: f64 = particle_1.dir.x;
-    let cosy: f64 = particle_1.dir.y;
-    let cosz: f64 = particle_1.dir.z;
-    let cphi: f64 = phi.cos();
-    let sphi: f64 = phi.sin();
-    let sa = (1. - cosx*cosx).sqrt();
-
-    //Particle direction update formulas from original TRIDYN paper, see Moeller and Eckstein 1988
-    let cpsi: f64 = psi.cos();
-    let spsi: f64 = psi.sin();
-    let cosx_new: f64 = cpsi*cosx + spsi*cphi*sa;
-    let cosy_new: f64 = cpsi*cosy - spsi/sa*(cphi*cosx*cosy - sphi*cosz);
-    let cosz_new: f64 = cpsi*cosz - spsi/sa*(cphi*cosx*cosz + sphi*cosy);
-
-    let dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new};
-
-    particle_1.dir.assign(&dir_new);
-    particle_1.dir.normalize();
-}
+    /// Rotate a particle by a deflection psi at an azimuthal angle phi.
+    pub fn rotate(&mut self, psi: f64, phi: f64) {
+        let cosx: f64 = self.dir.x;
+        let cosy: f64 = self.dir.y;
+        let cosz: f64 = self.dir.z;
+        let cphi: f64 = phi.cos();
+        let sphi: f64 = phi.sin();
+        let sa = (1. - cosx*cosx).sqrt();
+
+        //Particle direction update formulas from original TRIDYN paper, see Moeller and Eckstein 1988
+        let cpsi: f64 = psi.cos();
+        let spsi: f64 = psi.sin();
+        let cosx_new: f64 = cpsi*cosx + spsi*cphi*sa;
+        let cosy_new: f64 = cpsi*cosy - spsi/sa*(cphi*cosx*cosy - sphi*cosz);
+        let cosz_new: f64 = cpsi*cosz - spsi/sa*(cphi*cosx*cosz + sphi*cosy);
 
-/// Push particle in space according to previous direction and return the distance traveled.
-pub fn particle_advance(particle_1: &mut particle::Particle, mfp: f64, asymptotic_deflection: f64) -> f64 {
+        let dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new};
 
-    if particle_1.E > particle_1.Ec {
-        particle_1.add_trajectory();
+        self.dir.assign(&dir_new);
+        self.dir.normalize();
     }
 
-    //Update previous position
-    particle_1.pos_old.x = particle_1.pos.x;
-    particle_1.pos_old.y = particle_1.pos.y;
-    particle_1.pos_old.z = particle_1.pos.z;
+    /// Push particle in space according to previous direction and return the distance traveled.
+    pub fn advance(&mut self, mfp: f64, asymptotic_deflection: f64) -> f64 {
+
+        if self.E > self.Ec {
+            self.add_trajectory();
+        }
 
-    //In order to keep average denisty constant, must add back previous asymptotic deflection
-    let distance_traveled = mfp + particle_1.asymptotic_deflection - asymptotic_deflection;
-    //let distance_traveled = mfp - asymptotic_deflection;
+        //Update previous position
+        self.pos_old.x = self.pos.x;
+        self.pos_old.y = self.pos.y;
+        self.pos_old.z = self.pos.z;
 
-    //dir has been updated, so use previous direction to advance in space
-    particle_1.pos.x += particle_1.dir_old.x*distance_traveled;
-    particle_1.pos.y += particle_1.dir_old.y*distance_traveled;
-    particle_1.pos.z += particle_1.dir_old.z*distance_traveled;
-    particle_1.asymptotic_deflection = asymptotic_deflection;
+        //In order to keep average denisty constant, must add back previous asymptotic deflection
+        let distance_traveled = mfp + self.asymptotic_deflection - asymptotic_deflection;
+        //let distance_traveled = mfp - asymptotic_deflection;
 
-    //Update previous direction
-    particle_1.dir_old.x = particle_1.dir.x;
-    particle_1.dir_old.y = particle_1.dir.y;
-    particle_1.dir_old.z = particle_1.dir.z;
+        //dir has been updated, so use previous direction to advance in space
+        self.pos.x += self.dir_old.x*distance_traveled;
+        self.pos.y += self.dir_old.y*distance_traveled;
+        self.pos.z += self.dir_old.z*distance_traveled;
+        self.asymptotic_deflection = asymptotic_deflection;
 
-    return distance_traveled;
+        //Update previous direction
+        self.dir_old.x = self.dir.x;
+        self.dir_old.y = self.dir.y;
+        self.dir_old.z = self.dir.z;
+
+        return distance_traveled;
+    }
 }
 
+
+
 pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) {
     let E = particle.E;
 

src/structs.rs

@@ -44,18 +44,18 @@ impl Vector {
     }
 }
 
-/// Vector4 is a trajectory-tracking object that includes x, y, z, and the current energy.
+/// TrajectoryElement is a trajectory-tracking object that includes x, y, z, and the current energy.
 #[derive(Clone)]
-pub struct Vector4 {
+pub struct TrajectoryElement {
     pub E: f64,
     pub x: f64,
     pub y: f64,
     pub z: f64,
 }
 
-impl Vector4 {
-    pub fn new(E: f64, x: f64, y: f64, z: f64) -> Vector4 {
-        Vector4 {
+impl TrajectoryElement {
+    pub fn new(E: f64, x: f64, y: f64, z: f64) -> TrajectoryElement {
+        TrajectoryElement {
             E,
             x,
             y,

src/tests.rs

@@ -660,9 +660,6 @@ fn test_surface_refraction() {
     let cosy_new = cosy_new/dir_mag;
     let cosz_new = cosz_new/dir_mag;
 
-    //let delta_theta = particle::refraction_angle(cosx, E, E + Es);
-    //particle::rotate_particle(&mut particle_1, delta_theta, 0.);
-
     let normal = Vector::new(1.0, 0.0, 0.0);
     particle::surface_refraction(&mut particle_1, normal, Es);
 
@@ -690,9 +687,6 @@ fn test_surface_refraction() {
         println!("{} {} {}", particle_1.dir.x, particle_1.dir.y, particle_1.dir.z);
     }
 
-    //let delta_theta = particle::refraction_angle(particle_1.dir.x, particle_1.E, particle_1.E - Es);
-    //particle::rotate_particle(&mut particle_1, delta_theta, 0.);
-
     let normal = Vector::new(1.0, 0.0, 0.0);
     particle::surface_refraction(&mut particle_1, normal, -Es);
 
@@ -856,10 +850,10 @@ fn test_momentum_conservation() {
                         particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z);
 
                         //Rotate particle 1, 2 by lab frame scattering angles
-                        particle::rotate_particle(&mut particle_1, binary_collision_result.psi,
+                        particle_1.rotate(binary_collision_result.psi,
                             binary_collision_geometries[0].phi_azimuthal);
 
-                        particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil,
+                        particle_2.rotate(-binary_collision_result.psi_recoil,
                             binary_collision_geometries[0].phi_azimuthal);
 
                         //Subtract total energy from all simultaneous collisions and electronic stopping
@@ -894,7 +888,7 @@ fn test_momentum_conservation() {
 }
 
 #[test]
-fn test_rotate_particle() {
+fn test_rotate() {
     let mass = 1.;
     let Z = 1.;
     let E = 1.;
@@ -912,18 +906,18 @@ fn test_rotate_particle() {
     let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0);
 
     //Check that rotation in 2D works
-    particle::rotate_particle(&mut particle, psi, phi);
+    particle.rotate(psi, phi);
     assert!(approx_eq!(f64, particle.dir.x, 0., epsilon = 1E-12), "particle.dir.x: {} Should be ~0.", particle.dir.x);
     assert!(approx_eq!(f64, particle.dir.y, 1., epsilon = 1E-12), "particle.dir.y: {} Should be ~1.", particle.dir.y);
 
     //Check that rotating back by negative psi returns to the previous values
-    particle::rotate_particle(&mut particle, -psi, phi);
+    particle.rotate(-psi, phi);
     assert!(approx_eq!(f64, particle.dir.x, cosx, epsilon = 1E-12), "particle.dir.x: {} Should be ~{}", particle.dir.x, cosx);
     assert!(approx_eq!(f64, particle.dir.y, cosy, epsilon = 1E-12), "particle.dir.y: {} Should be ~{}", particle.dir.y, cosy);
 
     //Check that azimuthal rotation by 180 degrees works correctly
     let phi = PI;
-    particle::rotate_particle(&mut particle, psi, phi);
+    particle.rotate(psi, phi);
     assert!(approx_eq!(f64, particle.dir.x, 1., epsilon = 1E-12), "particle.dir.x: {} Should be ~1.", particle.dir.x);
     assert!(approx_eq!(f64, particle.dir.y, 0., epsilon = 1E-12), "particle.dir.y: {} Should be ~0.", particle.dir.y);
 
@@ -950,7 +944,7 @@ fn test_particle_advance() {
 
     let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0);
 
-    let distance_traveled = particle::particle_advance(&mut particle, mfp, asymptotic_deflection);
+    let distance_traveled = particle.advance(mfp, asymptotic_deflection);
 
     assert_eq!(particle.pos.x, (1. - 0.5)*cosx);
     assert_eq!(particle.pos.y, (1. - 0.5)*cosy);