Added tagged BCA output. Now particles carry tags and weights through…

.github/workflows/rustbca_compile_check.yml

@@ -15,6 +15,9 @@ jobs:
 
     steps:
     - uses: actions/checkout@v2
+    - name: udpate
+      run: |
+        sudo apt-get update
     - name: Install curl
       run: |
         sudo apt-get install curl

src/bca.rs

@@ -377,15 +377,19 @@ pub fn choose_collision_partner<T: Geometry>(particle_1: &particle::Particle, ma
 
     //Choose recoil Z, M
     let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil);
-
-    return (species_index,
-        particle::Particle::new(
-            M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil,
-            x_recoil, y_recoil, z_recoil,
-            cosx, cosy, cosz,
-            false, options.track_recoil_trajectories, interaction_index
-        )
-    )
+    let mut new_particle = particle::Particle::new(
+        M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil,
+        x_recoil, y_recoil, z_recoil,
+        cosx, cosy, cosz,
+        false, options.track_recoil_trajectories, interaction_index
+    );
+
+    // Carry through tag, weight, and tracked vector from originating particle
+    new_particle.weight = particle_1.weight;
+    new_particle.tag = particle_1.tag;
+    new_particle.tracked_vector = particle_1.tracked_vector;
+
+    return (species_index, new_particle)
 }
 
 /// Calculate the distance of closest approach of two particles given a particular binary collision geometry.

src/input.rs

@@ -351,7 +351,9 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
                             ux: cosx,
                             uy: cosy,
                             uz: cosz,
-                            interaction_index: interaction_index
+                            interaction_index: interaction_index,
+                            tag: 0,
+                            weight: 1.0,
                         }
                     );
                 }
@@ -395,7 +397,9 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
                             ux: cosx,
                             uy: cosy,
                             uz: cosz,
-                            interaction_index
+                            interaction_index,
+                            weight: 1.0,
+                            tag: 0,
                         }
                     );
                 }

src/lib.rs

@@ -123,6 +123,226 @@ pub struct InputCompoundBCA {
     pub Eb2: *mut f64,
 }
 
+#[derive(Debug)]
+#[repr(C)]
+pub struct InputTaggedBCA {
+    pub len: usize,
+    /// x y z
+    pub positions: *mut [f64; 3],
+    /// vx, vy, vz
+    pub velocities: *mut [f64; 3],
+    pub Z1: f64,
+    pub m1: f64,
+    pub Ec1: f64,
+    pub Es1: f64,
+    pub num_species_target: usize,
+    pub Z2: *mut f64,
+    pub m2: *mut f64,
+    pub n2: *mut f64,
+    pub Ec2: *mut f64,
+    pub Es2: *mut f64,
+    pub Eb2: *mut f64,
+    pub tags: *mut i32,
+    pub weights: *mut f64,
+}
+
+#[derive(Debug)]
+#[repr(C)]
+pub struct OutputTaggedBCA {
+    pub len: usize,
+    pub particles: *mut [f64; 9],
+    pub weights: *mut f64,
+    pub tags: *mut i32
+}
+
+#[no_mangle]
+pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTaggedBCA {
+
+    let mut total_output = vec![];
+    let mut output_tags = vec![];
+    let mut output_weights = vec![];
+
+    #[cfg(feature = "distributions")]
+    let options = Options {
+        name: "test".to_string(),
+        track_trajectories: false,
+        track_recoils: true,
+        track_recoil_trajectories: false,
+        write_buffer_size: 8000,
+        weak_collision_order: 3,
+        suppress_deep_recoils: false,
+        high_energy_free_flight_paths: false,
+        electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED,
+        mean_free_path_model: MeanFreePathModel::LIQUID,
+        interaction_potential: vec![vec![InteractionPotential::KR_C]],
+        scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]],
+        num_threads: 1,
+        num_chunks: 1,
+        use_hdf5: false,
+        root_finder: vec![vec![Rootfinder::DEFAULTNEWTON]],
+        track_displacements: false,
+        track_energy_losses: false,
+        energy_min: 0.0,
+        energy_max: 10.0,
+        energy_num: 11,
+        angle_min: 0.0,
+        angle_max: 90.0,
+        angle_num: 11,
+        x_min: 0.0,
+        y_min: -10.0,
+        z_min: -10.0,
+        x_max: 10.0,
+        y_max: 10.0,
+        z_max: 10.0,
+        x_num: 11,
+        y_num: 11,
+        z_num: 11,
+    };
+
+    #[cfg(not(feature = "distributions"))]
+    let options = Options {
+        name: "test".to_string(),
+        track_trajectories: false,
+        track_recoils: true,
+        track_recoil_trajectories: false,
+        write_buffer_size: 8000,
+        weak_collision_order: 3,
+        suppress_deep_recoils: false,
+        high_energy_free_flight_paths: false,
+        electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED,
+        mean_free_path_model: MeanFreePathModel::LIQUID,
+        interaction_potential: vec![vec![InteractionPotential::KR_C]],
+        scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]],
+        num_threads: 1,
+        num_chunks: 1,
+        use_hdf5: false,
+        root_finder: vec![vec![Rootfinder::DEFAULTNEWTON]],
+        track_displacements: false,
+        track_energy_losses: false,
+    };
+
+    let Z2 = unsafe { slice::from_raw_parts(input.Z2, input.num_species_target).to_vec() };
+    let m2 = unsafe { slice::from_raw_parts(input.m2, input.num_species_target).to_vec() };
+    let n2 = unsafe { slice::from_raw_parts(input.n2, input.num_species_target).to_vec() };
+    let Ec2 = unsafe { slice::from_raw_parts(input.Ec2, input.num_species_target).to_vec() };
+    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 positions = unsafe { slice::from_raw_parts(input.positions, input.num_species_target).to_vec() };
+    let tags = unsafe { slice::from_raw_parts(input.tags, input.num_species_target).to_vec() };
+    let weights = unsafe { slice::from_raw_parts(input.weights, input.num_species_target).to_vec() };
+
+    let x = -2.*(n2.iter().sum::<f64>()*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,
+        Z: Z2,
+        m: m2,
+        interaction_index: vec![0; input.num_species_target],
+        surface_binding_model: SurfaceBindingModel::AVERAGE,
+        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::<Mesh0D>::new(&material_parameters, &geometry_input);
+
+    let velocities = unsafe { slice::from_raw_parts(input.velocities, input.len) };
+
+    let mut index: usize = 0;
+    for velocity in velocities {
+
+        let vx = velocity[0];
+        let vy = velocity[1];
+        let vz = velocity[2];
+
+        let v = (vx*vx + vy*vy + vz*vz).sqrt();
+
+        let E1 = 0.5*input.m1*AMU*v*v;
+
+        let ux = vx/v;
+        let uy = vy/v;
+        let uz = vz/v;
+
+        let p = particle::Particle {
+            m: input.m1*AMU,
+            Z: input.Z1,
+            E: E1,
+            Ec: input.Ec1*EV,
+            Es: input.Es1*EV,
+            pos: Vector::new(x, y, z),
+            dir: Vector::new(ux, uy, uz),
+            pos_origin: Vector::new(x, y, z),
+            pos_old: Vector::new(x, y, z),
+            dir_old: Vector::new(ux, uy, uz),
+            energy_origin: E1,
+            asymptotic_deflection: 0.0,
+            stopped: false,
+            left: false,
+            incident: true,
+            first_step: true,
+            trajectory: vec![],
+            energies: vec![],
+            track_trajectories: false,
+            number_collision_events: 0,
+            backreflected: false,
+            interaction_index : 0,
+            weight: weights[index],
+            tag: tags[index],
+            tracked_vector: Vector::new(positions[index][0], positions[index][1], positions[index][2]),
+        };
+
+        let output = bca::single_ion_bca(p, &m, &options);
+
+        for particle in output {
+
+            if (particle.left) | (particle.incident) {
+                total_output.push(
+                    [
+                        particle.Z,
+                        particle.m/AMU,
+                        particle.E/EV,
+
+                        particle.tracked_vector.x/ANGSTROM,
+                        particle.tracked_vector.y/ANGSTROM,
+                        particle.tracked_vector.z/ANGSTROM,
+
+                        particle.dir.x,
+                        particle.dir.y,
+                        particle.dir.z
+                    ]
+                );
+                output_tags.push(particle.tag);
+                output_weights.push(particle.weight);
+            }
+        }
+        index += 1;
+    }
+
+    let len = total_output.len();
+    let particles = total_output.as_mut_ptr();
+    let tags_ptr = output_tags.as_mut_ptr();
+    let weights_ptr = output_weights.as_mut_ptr();
+
+    std::mem::forget(total_output);
+    OutputTaggedBCA {
+        len,
+        particles,
+        tags: tags_ptr,
+        weights: weights_ptr
+    }
+}
+
+
 #[no_mangle]
 pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA {
 
@@ -250,7 +470,10 @@ pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA {
             track_trajectories: false,
             number_collision_events: 0,
             backreflected: false,
-            interaction_index : 0
+            interaction_index : 0,
+            weight: 1.0,
+            tag: 0,
+            tracked_vector: Vector::new(0.0, 0.0, 0.0),
         };
 
 
@@ -420,7 +643,10 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA {
             track_trajectories: false,
             number_collision_events: 0,
             backreflected: false,
-            interaction_index : 0
+            interaction_index : 0,
+            weight: 1.0,
+            tag: 0,
+            tracked_vector: Vector::new(0.0, 0.0, 0.0),
         };
 
 
@@ -456,6 +682,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA {
     }
 }
 
+
 #[no_mangle]
 pub extern "C" fn simple_bca_c(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, Z1: f64, m1: f64, Ec1: f64, Es1: f64, Z2: f64, m2: f64, Ec2: f64, Es2: f64, n2: f64, Eb2: f64) -> OutputBCA {
     let mut output = simple_bca(x, y, z, ux, uy, uz, E1, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2);
@@ -585,7 +812,10 @@ pub fn simple_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, Z1
         track_trajectories: false,
         number_collision_events: 0,
         backreflected: false,
-        interaction_index : 0
+        interaction_index : 0,
+        weight: 1.0,
+        tag: 0,
+        tracked_vector: Vector::new(0.0, 0.0, 0.0),
     };
 
     let material_parameters = material::MaterialParameters {

src/particle.rs

@@ -18,7 +18,7 @@ pub struct ParticleParameters {
     pub Es: Vec<f64>,
     pub pos: Vec<(f64, f64, f64)>,
     pub dir: Vec<(f64, f64, f64)>,
-    pub interaction_index: Vec<usize>
+    pub interaction_index: Vec<usize>,
 }
 
 #[cfg(not(feature = "hdf5_input"))]
@@ -35,7 +35,7 @@ pub struct ParticleParameters {
     pub Es: Vec<f64>,
     pub pos: Vec<(f64, f64, f64)>,
     pub dir: Vec<(f64, f64, f64)>,
-    pub interaction_index: Vec<usize>
+    pub interaction_index: Vec<usize>,
 }
 
 /// HDF5 version of particle input.
@@ -55,6 +55,8 @@ pub struct ParticleInput {
     pub uy: f64,
     pub uz: f64,
     pub interaction_index: usize,
+    pub weight: f64,
+    pub tag: i32,
 }
 
 /// Particle object. Particles in rustbca include incident ions and material atoms.
@@ -82,6 +84,9 @@ pub struct Particle {
     pub number_collision_events: usize,
     pub backreflected: bool,
     pub interaction_index: usize,
+    pub weight: f64,
+    pub tag: i32,
+    pub tracked_vector: Vector,
 }
 impl Particle {
     /// Construct a particle object from input.
@@ -118,6 +123,9 @@ impl Particle {
             number_collision_events: 0,
             backreflected: false,
             interaction_index: input.interaction_index,
+            weight: 1.0,
+            tag: 0,
+            tracked_vector: Vector::new(0.0, 0.0, 0.0),
         }
     }
 
@@ -148,6 +156,9 @@ impl Particle {
             number_collision_events: 0,
             backreflected: false,
             interaction_index,
+            weight: 1.0,
+            tag: 0,
+            tracked_vector: Vector::new(0.0, 0.0, 0.0),
         }
     }
 

src/structs.rs

@@ -1,5 +1,5 @@
 /// 3D vector.
-#[derive(Clone)]
+#[derive(Clone, Copy)]
 pub struct Vector {
     pub x: f64,
     pub y: f64,