From c7dae154661b562771fe85b5de73bcaf2334e6a3 Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 08:45:21 -0700
Subject: [PATCH 01/12] Draft of generalized RustBCA python wrapper(s)

---
 Cargo.toml                          |  3 +-
 examples/make_input_file_and_run.py | 14 +++++
 src/input.rs                        | 18 ++++--
 src/lib.rs                          | 54 +++++++++++++++++-
 src/main.rs                         | 15 +++++
 src/output.rs                       |  1 +
 src/physics.rs                      | 86 +++++++++++++++++++++++++++--
 7 files changed, 179 insertions(+), 12 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 6b2c3744..6933a135 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -29,6 +29,7 @@ rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
 ndarray = {version = "0.17.2", features = ["serde"], optional = true}
 parry3d-f64 = {optional = true, version="0.2.0"}
 pyo3 = {version = "0.29.0", optional=true}
+pythonize = {version = "0.29.0", optional=true}
 
 [dev-dependencies]
 float-cmp = "0.10.0"
@@ -45,5 +46,5 @@ cpr_rootfinder = ["rcpr"]
 distributions = ["ndarray"]
 no_list_output = []
 parry3d = ["parry3d-f64"]
-python = ["pyo3"]
+python = ["pyo3", "pythonize"]
 extended_max_z = []
\ No newline at end of file
diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py
index b48d376c..37e3a161 100644
--- a/examples/make_input_file_and_run.py
+++ b/examples/make_input_file_and_run.py
@@ -255,6 +255,17 @@
     'geometry_input': geometry_0D
 }
 
+input_data['options']['name'] = 'rustbca_input_file'
+rustbca_py(input_data, mode)
+s = np.genfromtxt('rustbca_input_filesputtered.output', delimiter=',')
+
+arrays = rustbca_local_py(input_data, mode)
+sputtered = arrays['sputtered']
+
+np.testing.assert_approx_equal(s[0, 2], np.array(arrays['energy'])[sputtered][0])
+
+input_data['options']['name'] = 'input_file'
+
 # Attempt to cleanup line endings
 input_string = dumps(input_data).replace('\r', '')
 with  open('examples/input_file.toml', 'w') as input_file:
@@ -265,6 +276,9 @@
 
 # Read output files - ensure arrays are at least 2D for indexing
 sputtered = np.atleast_2d(np.genfromtxt('input_filesputtered.output', delimiter=','))
+
+np.testing.assert_approx_equal(sputtered[0, 2], s[0, 2])
+
 reflected = np.atleast_2d(np.genfromtxt('input_filereflected.output', delimiter=','))
 implanted = np.atleast_2d(np.genfromtxt('input_filedeposited.output', delimiter=','))
 
diff --git a/src/input.rs b/src/input.rs
index 02665b17..e6c607f1 100644
--- a/src/input.rs
+++ b/src/input.rs
@@ -381,9 +381,7 @@ impl Options {
     }
 }
 
-pub fn input<T: Geometry>(input_file: String) -> (Vec<particle::ParticleInput>, material::Material<T>, Options, OutputUnits)
-where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
-
+pub fn read_input_file<T: Geometry>(input_file: String) -> <T as Geometry>::InputFileFormat {
     //Read input file, convert to string, and open with toml
     let mut input_toml = String::new();
     let mut file = OpenOptions::new()
@@ -394,7 +392,19 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
         .unwrap_or_else(|_| panic!("Input errror: could not open input file {}.", &input_file));
     file.read_to_string(&mut input_toml).context("Could not convert TOML file to string.").unwrap();
 
-    let input: <T as Geometry>::InputFileFormat = InputFile::new(&input_toml);
+    InputFile::new(&input_toml)
+}
+
+pub fn input<T: Geometry>(input_file: String) -> (Vec<particle::ParticleInput>, material::Material<T>, Options, OutputUnits)
+where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
+
+    let input: <T as Geometry>::InputFileFormat = read_input_file::<T>(input_file);
+
+    process_input_file(input)
+
+}
+
+pub fn process_input_file<T: Geometry>(input: <T as Geometry>::InputFileFormat) -> (Vec<particle::ParticleInput>, material::Material<T>, Options, OutputUnits) {
 
     //Unpack toml information into structs
     let options = (*input.get_options()).clone();
diff --git a/src/lib.rs b/src/lib.rs
index 9d24cae9..003689d2 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -53,6 +53,8 @@ use std::f64::consts::SQRT_2;
 use pyo3::prelude::*;
 #[cfg(feature = "python")]
 use pyo3::types::*;
+#[cfg(feature = "python")]
+use pythonize::*;
 
 //Load internal modules
 pub mod material;
@@ -68,9 +70,7 @@ pub mod consts;
 pub mod structs;
 pub mod sphere;
 pub mod math;
-
-#[cfg(feature = "parry3d")]
-pub mod parry;
+pub mod physics;
 
 pub use crate::enums::*;
 pub use crate::consts::*;
@@ -81,6 +81,10 @@ pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D};
 pub use crate::sphere::{Sphere, SphereInput, InputSphere};
 pub use crate::math::*;
 pub use crate::material::*;
+pub use crate::physics::*;
+
+#[cfg(feature = "parry3d")]
+pub mod parry;
 
 #[cfg(feature = "parry3d")]
 pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, ParryTriMeshInput, InputParryTriMesh};
@@ -139,6 +143,12 @@ mod libRustBCA {
 
     #[pymodule_export]
     use super::scattering_integrals;
+
+    #[pymodule_export]
+    use super::rustbca_py;
+
+    #[pymodule_export]
+    use super::rustbca_local_py;
 }
 
 #[derive(Debug)]
@@ -2177,4 +2187,42 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g
     let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C);
 
     (theta_gm, theta_gl, theta_mw, theta_magic)
+}
+
+#[cfg(feature = "python")]
+#[pyfunction]
+#[pyo3(signature=(input, geometry_type="1D", output_path=""))]
+fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_type: &str, output_path: &str) {
+
+    match geometry_type {
+        "1D" => {
+            let input: <Mesh1D as geometry::Geometry>::InputFileFormat = depythonize(&input).unwrap();
+            let (particle_input_array, material, options, output_units) = input::process_input_file(input);
+            let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build().unwrap();
+            pool.install( ||
+                physics::physics_loop::<Mesh1D>(particle_input_array, material, options, output_units)
+            );
+        }
+        _ => panic!("Input Error: unimplemented geometry mode for rustbca_py. Try `1D`")
+    }
+}
+
+#[cfg(feature = "python")]
+#[pyfunction]
+#[pyo3(signature=(input, geometry_type="1D"))]
+fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_type: &str) -> PyResult<Bound<'py, PyAny>> {
+
+    match geometry_type {
+        "1D" => {
+            let input: <Mesh1D as geometry::Geometry>::InputFileFormat = depythonize(&input).unwrap();
+            let (particle_input_array, material, options, output_units) = input::process_input_file(input);
+            let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build().unwrap();
+            let finished_particles = pool.install(||
+                physics::silent_physics_loop::<Mesh1D>(particle_input_array, material, options, output_units.clone())
+            );
+            let finished_particles_container = physics::process_finished_particles_to_arrays(finished_particles, output_units);
+            Ok(pythonize(python, &finished_particles_container)?)
+        }
+        _ => panic!("Input Error: unimplemented geometry mode for rustbca_py. Try `1D`")
+    }
 }
\ No newline at end of file
diff --git a/src/main.rs b/src/main.rs
index 8b02d4be..67a145dc 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -94,32 +94,47 @@ fn main() {
      match geometry_type {
         GeometryType::MESH0D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh0D>(input_file);
+            //Initialize threads with rayon
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<Mesh0D>(particle_input_array, material, options, output_units);
         },
         GeometryType::MESH1D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh1D>(input_file);
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<Mesh1D>(particle_input_array, material, options, output_units);
         },
         GeometryType::MESH2D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh2D>(input_file);
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<Mesh2D>(particle_input_array, material, options, output_units);
         },
         GeometryType::SPHERE => {
             let (particle_input_array, material, options, output_units) = input::input::<Sphere>(input_file);
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<Sphere>(particle_input_array, material, options, output_units);
         },
         #[cfg(feature = "parry3d")]
         GeometryType::BALL => {
             let (particle_input_array, material, options, output_units) = input::input::<ParryBall>(input_file);
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<ParryBall>(particle_input_array, material, options, output_units);
         }
         #[cfg(feature = "parry3d")]
         GeometryType::TRIMESH => {
             let (particle_input_array, material, options, output_units) = input::input::<ParryTriMesh>(input_file);
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<ParryTriMesh>(particle_input_array, material, options, output_units);
         }
         GeometryType::HOMOGENEOUS2D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::HomogeneousMesh2D>(input_file);
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
             physics_loop::<HomogeneousMesh2D>(particle_input_array, material, options, output_units);
         }
     }
diff --git a/src/output.rs b/src/output.rs
index ad3d7a2e..463c0c50 100644
--- a/src/output.rs
+++ b/src/output.rs
@@ -1,6 +1,7 @@
 use super::*;
 use std::fs::File;
 
+#[derive(Clone, Debug)]
 pub struct OutputUnits {
     pub length_unit: f64,
     pub energy_unit: f64,
diff --git a/src/physics.rs b/src/physics.rs
index 933d358e..cfd37882 100644
--- a/src/physics.rs
+++ b/src/physics.rs
@@ -1,4 +1,86 @@
 use super::*;
+use rayon::iter::*;
+use indicatif::{ProgressBar, ProgressStyle};
+
+pub fn silent_physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) -> Vec<particle::Particle> {
+    
+    let mut finished_particles: Vec<particle::Particle> = Vec::new();
+
+    finished_particles.par_extend(
+        particle_input_array.into_par_iter()
+        .enumerate()
+        .map_init(
+            || if options.seed < 0 { 
+                ChaCha8Rng::seed_from_u64(rand::random())
+            } else {
+                ChaCha8Rng::seed_from_u64(u64::try_from(options.seed).expect("Value Error: seed not u64."))
+            },
+            | rng, (particle_index, particle_input)| {
+                rng.set_stream((particle_index) as u64);
+                bca::single_ion_bca(particle::Particle::from_input(particle_input, &options), &material, &options, rng)
+        }).flatten()
+    );
+    finished_particles
+}
+
+#[derive(Clone, Debug, Serialize, Deserialize)]
+pub struct FinishedParticlesContainer {
+    sputtered: Vec<bool>,
+    implanted: Vec<bool>,
+    atomic_number: Vec<usize>,
+    mass: Vec<f64>,
+    energy: Vec<f64>,
+    x: Vec<f64>,
+    y: Vec<f64>,
+    z: Vec<f64>,
+    ux: Vec<f64>,
+    uy: Vec<f64>,
+    uz: Vec<f64>,
+}
+impl FinishedParticlesContainer {
+    pub fn new() -> FinishedParticlesContainer {
+        FinishedParticlesContainer {
+            sputtered: Vec::new(),
+            implanted: Vec::new(),
+            atomic_number: Vec::new(),
+            mass: Vec::new(),
+            energy: Vec::new(),
+            x: Vec::new(),
+            y: Vec::new(),
+            z: Vec::new(),
+            ux: Vec::new(),
+            uy: Vec::new(),
+            uz: Vec::new(),
+        }
+    }
+}
+
+
+pub fn process_finished_particles_to_arrays(finished_particles: Vec<particle::Particle>, output_units: OutputUnits) -> FinishedParticlesContainer {
+    let mut finished_particles_container = FinishedParticlesContainer::new();
+    for particle in finished_particles {
+
+        let implanted = particle.incident & !particle.left;
+        let sputtered = !particle.incident & particle.left;
+        let reflected = particle.incident & particle.left;
+
+        if implanted | sputtered | reflected {
+            finished_particles_container.sputtered.push(sputtered);
+            finished_particles_container.implanted.push(implanted);
+
+            finished_particles_container.atomic_number.push(particle.Z as usize);
+            finished_particles_container.mass.push(particle.m/output_units.mass_unit);
+            finished_particles_container.energy.push(particle.E/output_units.energy_unit);
+            finished_particles_container.x.push(particle.pos.x/output_units.length_unit);
+            finished_particles_container.y.push(particle.pos.y/output_units.length_unit);
+            finished_particles_container.z.push(particle.pos.z/output_units.length_unit);
+            finished_particles_container.ux.push(particle.dir.x);
+            finished_particles_container.uy.push(particle.dir.y);
+            finished_particles_container.uz.push(particle.dir.z);
+        }
+    }
+    finished_particles_container
+}
 
 pub fn physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) {
 
@@ -15,10 +97,6 @@ pub fn physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::Part
         #[cfg(feature = "distributions")]
         let mut distributions = output::Distributions::new(&options);
 
-        //Initialize threads with rayon
-        println!("Initializing with {} threads...", options.num_threads);
-        rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();
-
         //Create and configure progress bar
         let bar: ProgressBar = ProgressBar::new(total_count);
         bar.set_style(ProgressStyle::default_bar()

From 7256a65cd118cdb9e2bac0d14322b3be586c84ce Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 08:47:37 -0700
Subject: [PATCH 02/12]  remove output_dir arg for now

---
 src/lib.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 003689d2..84c3d5dc 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -2191,8 +2191,8 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g
 
 #[cfg(feature = "python")]
 #[pyfunction]
-#[pyo3(signature=(input, geometry_type="1D", output_path=""))]
-fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_type: &str, output_path: &str) {
+#[pyo3(signature=(input, geometry_type="1D"))]
+fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_type: &str) {
 
     match geometry_type {
         "1D" => {

From 1c104aa36b2a95d50ee0ba41be268238ab4d50a5 Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 08:54:36 -0700
Subject: [PATCH 03/12] Shuffle around terminal print statements to match
 previous behavior

---
 src/main.rs    | 6 ++++++
 src/physics.rs | 4 +---
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/main.rs b/src/main.rs
index 67a145dc..8b9895af 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -103,18 +103,21 @@ fn main() {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh1D>(input_file);
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Mesh1D>(particle_input_array, material, options, output_units);
         },
         GeometryType::MESH2D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh2D>(input_file);
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Mesh2D>(particle_input_array, material, options, output_units);
         },
         GeometryType::SPHERE => {
             let (particle_input_array, material, options, output_units) = input::input::<Sphere>(input_file);
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Sphere>(particle_input_array, material, options, output_units);
         },
         #[cfg(feature = "parry3d")]
@@ -122,6 +125,7 @@ fn main() {
             let (particle_input_array, material, options, output_units) = input::input::<ParryBall>(input_file);
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<ParryBall>(particle_input_array, material, options, output_units);
         }
         #[cfg(feature = "parry3d")]
@@ -129,12 +133,14 @@ fn main() {
             let (particle_input_array, material, options, output_units) = input::input::<ParryTriMesh>(input_file);
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<ParryTriMesh>(particle_input_array, material, options, output_units);
         }
         GeometryType::HOMOGENEOUS2D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::HomogeneousMesh2D>(input_file);
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<HomogeneousMesh2D>(particle_input_array, material, options, output_units);
         }
     }
diff --git a/src/physics.rs b/src/physics.rs
index cfd37882..443444c2 100644
--- a/src/physics.rs
+++ b/src/physics.rs
@@ -3,7 +3,7 @@ use rayon::iter::*;
 use indicatif::{ProgressBar, ProgressStyle};
 
 pub fn silent_physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) -> Vec<particle::Particle> {
-    
+
     let mut finished_particles: Vec<particle::Particle> = Vec::new();
 
     finished_particles.par_extend(
@@ -84,8 +84,6 @@ pub fn process_finished_particles_to_arrays(finished_particles: Vec<particle::Pa
 
 pub fn physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) {
 
-        println!("Processing {} ions...", particle_input_array.len());
-
         let total_count: u64 = particle_input_array.len() as u64;
         assert!(total_count/options.num_chunks > 0, "Input error: chunk size == 0 - reduce num_chunks or increase particle count.");
 

From d661320bdb2c582f8c47031bd4297df43e9c05a7 Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 08:56:31 -0700
Subject: [PATCH 04/12] add missing println

---
 src/main.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/main.rs b/src/main.rs
index 8b9895af..3eeb3837 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -97,6 +97,7 @@ fn main() {
             //Initialize threads with rayon
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Mesh0D>(particle_input_array, material, options, output_units);
         },
         GeometryType::MESH1D => {

From 22d9e2c681209d27523384970c2ad1c5ae1219f9 Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 09:02:28 -0700
Subject: [PATCH 05/12] Wrong order for print statements

---
 src/main.rs | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/main.rs b/src/main.rs
index 3eeb3837..8a3944a0 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -95,53 +95,53 @@ fn main() {
         GeometryType::MESH0D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh0D>(input_file);
             //Initialize threads with rayon
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Mesh0D>(particle_input_array, material, options, output_units);
         },
         GeometryType::MESH1D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh1D>(input_file);
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Mesh1D>(particle_input_array, material, options, output_units);
         },
         GeometryType::MESH2D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh2D>(input_file);
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Mesh2D>(particle_input_array, material, options, output_units);
         },
         GeometryType::SPHERE => {
             let (particle_input_array, material, options, output_units) = input::input::<Sphere>(input_file);
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<Sphere>(particle_input_array, material, options, output_units);
         },
         #[cfg(feature = "parry3d")]
         GeometryType::BALL => {
             let (particle_input_array, material, options, output_units) = input::input::<ParryBall>(input_file);
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<ParryBall>(particle_input_array, material, options, output_units);
         }
         #[cfg(feature = "parry3d")]
         GeometryType::TRIMESH => {
             let (particle_input_array, material, options, output_units) = input::input::<ParryTriMesh>(input_file);
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<ParryTriMesh>(particle_input_array, material, options, output_units);
         }
         GeometryType::HOMOGENEOUS2D => {
             let (particle_input_array, material, options, output_units) = input::input::<geometry::HomogeneousMesh2D>(input_file);
+            println!("Processing {} ions...", particle_input_array.len());
             println!("Initializing with {} threads...", options.num_threads);
             let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            println!("Processing {} ions...", particle_input_array.len());
             physics_loop::<HomogeneousMesh2D>(particle_input_array, material, options, output_units);
         }
     }

From 7811a286ead62a00cf9bdcfe94ff9ec996d0c1e7 Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 09:29:42 -0700
Subject: [PATCH 06/12] Added tests of rustbca_py to test_different_options.py

---
 examples/test_different_options.py | 15 +++++++++++++++
 src/lib.rs                         | 12 ++++++------
 2 files changed, 21 insertions(+), 6 deletions(-)

diff --git a/examples/test_different_options.py b/examples/test_different_options.py
index 587d6ce1..f3ca6290 100644
--- a/examples/test_different_options.py
+++ b/examples/test_different_options.py
@@ -51,6 +51,8 @@
 a TOML file using tomlkit.
 
 It runs the input file with cargo run --release and reads the output files.
+
+It also runs with rustbca_py and ensures the output files are identical.
 '''
 
 def run_test(
@@ -275,6 +277,19 @@ def run_test(
     reflected = np.atleast_2d(np.genfromtxt(f'input_file_{index}reflected.output', delimiter=','))
     implanted = np.atleast_2d(np.genfromtxt(f'input_file_{index}deposited.output', delimiter=','))
 
+    input_data['options']['name'] = f'python_{index}'
+    if run_sim:
+        rustbca_py(input_data, geometry_mode=mode)
+
+    # Read output files - ensure arrays are at least 2D for indexing
+    sputtered_py = np.atleast_2d(np.genfromtxt(f'python_{index}sputtered.output', delimiter=','))
+    reflected_py = np.atleast_2d(np.genfromtxt(f'python_{index}reflected.output', delimiter=','))
+    implanted_py = np.atleast_2d(np.genfromtxt(f'python_{index}deposited.output', delimiter=','))
+
+    np.testing.assert_allclose(sputtered, sputtered_py)
+    np.testing.assert_allclose(reflected, reflected_py)
+    np.testing.assert_allclose(implanted, implanted_py)
+
     return sputtered, reflected, implanted
 
 num_bins = 75
diff --git a/src/lib.rs b/src/lib.rs
index 84c3d5dc..d5b6b135 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -2191,10 +2191,10 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g
 
 #[cfg(feature = "python")]
 #[pyfunction]
-#[pyo3(signature=(input, geometry_type="1D"))]
-fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_type: &str) {
+#[pyo3(signature=(input, geometry_mode="1D"))]
+fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_mode: &str) {
 
-    match geometry_type {
+    match geometry_mode {
         "1D" => {
             let input: <Mesh1D as geometry::Geometry>::InputFileFormat = depythonize(&input).unwrap();
             let (particle_input_array, material, options, output_units) = input::process_input_file(input);
@@ -2209,10 +2209,10 @@ fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_type: &str) {
 
 #[cfg(feature = "python")]
 #[pyfunction]
-#[pyo3(signature=(input, geometry_type="1D"))]
-fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_type: &str) -> PyResult<Bound<'py, PyAny>> {
+#[pyo3(signature=(input, geometry_mode="1D"))]
+fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<Bound<'py, PyAny>> {
 
-    match geometry_type {
+    match geometry_mode {
         "1D" => {
             let input: <Mesh1D as geometry::Geometry>::InputFileFormat = depythonize(&input).unwrap();
             let (particle_input_array, material, options, output_units) = input::process_input_file(input);

From 2f34ef6179a0c266c51e91bb2f0b64216bd94c31 Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 11:32:48 -0700
Subject: [PATCH 07/12] Implemented better python error handling.

---
 src/lib.rs | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index d5b6b135..31db22b5 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -55,6 +55,8 @@ use pyo3::prelude::*;
 use pyo3::types::*;
 #[cfg(feature = "python")]
 use pythonize::*;
+#[cfg(feature = "python")]
+use pyo3::exceptions::PyValueError;
 
 //Load internal modules
 pub mod material;
@@ -2192,7 +2194,7 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g
 #[cfg(feature = "python")]
 #[pyfunction]
 #[pyo3(signature=(input, geometry_mode="1D"))]
-fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_mode: &str) {
+fn rustbca_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<()> {
 
     match geometry_mode {
         "1D" => {
@@ -2202,8 +2204,9 @@ fn rustbca_py<'py>(input: &Bound<'py, PyDict>, geometry_mode: &str) {
             pool.install( ||
                 physics::physics_loop::<Mesh1D>(particle_input_array, material, options, output_units)
             );
-        }
-        _ => panic!("Input Error: unimplemented geometry mode for rustbca_py. Try `1D`")
+            Ok(())
+        },
+       _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode)))
     }
 }
 
@@ -2223,6 +2226,6 @@ fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geomet
             let finished_particles_container = physics::process_finished_particles_to_arrays(finished_particles, output_units);
             Ok(pythonize(python, &finished_particles_container)?)
         }
-        _ => panic!("Input Error: unimplemented geometry mode for rustbca_py. Try `1D`")
+        _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode)))
     }
 }
\ No newline at end of file

From e9216119fb8487446d74ac1a59f1097633171e9b Mon Sep 17 00:00:00 2001
From: Jon Drobny <contact@jondrobny.com>
Date: Wed, 5 Aug 2026 11:49:40 -0700
Subject: [PATCH 08/12] Implemented better error handling on electronic
 stopping function - this is probable the model going forward, and it probably
 needs its own issue.

---
 src/lib.rs | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 31db22b5..51ba74c1 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -56,7 +56,7 @@ use pyo3::types::*;
 #[cfg(feature = "python")]
 use pythonize::*;
 #[cfg(feature = "python")]
-use pyo3::exceptions::PyValueError;
+use pyo3::exceptions::{PyValueError, PyRuntimeError};
 
 //Load internal modules
 pub mod material;
@@ -2143,7 +2143,6 @@ pub fn compound_reflection_coefficient<'py>(ion: &Bound<'py, PyDict>, targets: V
 
                 let mut residue = residue.lock().unwrap();
                 *residue = *residue + residue_part;
-
             }
         }
     });
@@ -2173,22 +2172,33 @@ fn moller_knuth_two_sum(a: f64, b: f64) -> (f64, f64) {
     let r = delta_a + delta_b;
     (s, r)
 }
+
 #[cfg(feature = "python")]
 #[pyfunction]
-#[pyo3(signature = (Za, Zb, Ma, Mb, E0, p, n_gl_points=100))]
-fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_gl_points: usize) -> (f64, f64, f64, f64) {
+#[pyo3(signature = (Za, Zb, Ma, Mb, E0, p, n_gl_points=100, interaction_potential="KR_C"))]
+fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_gl_points: usize, interaction_potential: &str) -> PyResult<(f64, f64, f64, f64)> {
     let E0 = E0*EV;
     let p = p*ANGSTROM;
 
-    let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 1000, 1E-12).unwrap();
+    let potential = match interaction_potential {
+        "KR_C" => InteractionPotential::KR_C,
+        "LENZ_JENSEN" => InteractionPotential::LENZ_JENSEN,
+        "MOLIERE" => InteractionPotential::MOLIERE,
+        "ZBL" => InteractionPotential::ZBL,
+        _ => return Err(PyValueError::new_err(format!("Unimplemented interaction potential {}; try 'KR_C'", interaction_potential)))
+    };
+
+    let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, potential, 1000, 1E-12).map_err(
+        |error| PyRuntimeError::new_err(format!("Rootfinder failed to find distance of closest approach; check input values."))
+    )?;
 
     //Compute center of mass deflection angle with each algorithm
-    let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C, n_gl_points);
-    let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C);
-    let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C);
-    let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C);
+    let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, potential, n_gl_points);
+    let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, potential);
+    let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, potential);
+    let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, potential);
 
-    (theta_gm, theta_gl, theta_mw, theta_magic)
+    Ok((theta_gm, theta_gl, theta_mw, theta_magic))
 }
 
 #[cfg(feature = "python")]

From df104cfdd656e9bd590f77b2e3b257fcbe115efc Mon Sep 17 00:00:00 2001
From: Jon Drobny <drobny2@illinois.edu>
Date: Wed, 5 Aug 2026 13:18:41 -0700
Subject: [PATCH 09/12] Add different interaction potentials to
 scattering_integrals()

---
 examples/test_scattering_integrals.py | 45 +++++++++++++++------------
 1 file changed, 25 insertions(+), 20 deletions(-)

diff --git a/examples/test_scattering_integrals.py b/examples/test_scattering_integrals.py
index 7e5d9c4e..2a43da15 100644
--- a/examples/test_scattering_integrals.py
+++ b/examples/test_scattering_integrals.py
@@ -10,7 +10,7 @@
 from formulas import *
 
 energies = np.logspace(0, 4, 4)
-impact_parameters = np.logspace(-3, 3, 100)
+impact_parameters = np.logspace(-3, 2, 100)
 
 ion = helium
 target = boron
@@ -24,23 +24,28 @@
 
 linestyles = ['-', '--', ':', '-.']
 
-for linestyle, energy in zip(linestyles, energies):
-    gm = np.zeros(100)
-    gl = np.zeros(100)
-    mw = np.zeros(100)
-    magic = np.zeros(100)
-    for index, p in enumerate(impact_parameters):
-        gm[index], gl[index], mw[index], magic[index] = scattering_integrals(Za, Zb, Ma, Mb, energy, p)
-    plt.semilogx(impact_parameters, gm, label=f'Gauss-Mehler, E={np.round(energy/1000)} keV', linestyle=linestyle)
-    plt.semilogx(impact_parameters, gl, label=f'Gauss-Legendre, E={np.round(energy/1000)} keV', linestyle=linestyle)
-    plt.semilogx(impact_parameters, mw, label=f'Mendenhall-Weller, E={np.round(energy/1000)} keV', linestyle=linestyle)
-    plt.semilogx(impact_parameters, magic, label=f'MAGIC, E={np.round(energy/1000)} keV', linestyle=linestyle)
-    plt.gca().set_prop_cycle(None)
-
-    np.testing.assert_allclose(gm, gl, atol=5e-3) # 0.5% seems reasonable? max is ~0.3%
-    np.testing.assert_allclose(gm, mw, atol=5e-3)
-    np.testing.assert_allclose(mw, gl, atol=5e-3)
-    plt.legend()
-    plt.xlabel('p [A]')
-    plt.ylabel('theta [rad]')
+for potential in ["KR_C", "MOLIERE", "ZBL"]:
+    plt.figure()
+    plt.title(f'Scattering Angles for {potential}')
+    for linestyle, energy in zip(linestyles, energies):
+        gm = np.zeros(100)
+        gl = np.zeros(100)
+        mw = np.zeros(100)
+        magic = np.zeros(100)
+        for index, p in enumerate(impact_parameters):
+            gm[index], gl[index], mw[index], magic[index] = scattering_integrals(Za, Zb, Ma, Mb, energy, p, interaction_potential=potential)
+
+        plt.semilogx(impact_parameters, gm, label=f'Gauss-Mehler, E={np.round(energy/1000, 3)} keV', linestyle=linestyle)
+        plt.semilogx(impact_parameters, gl, label=f'Gauss-Legendre, E={np.round(energy/1000, 3)} keV', linestyle=linestyle)
+        plt.semilogx(impact_parameters, mw, label=f'Mendenhall-Weller, E={np.round(energy/1000, 3)} keV', linestyle=linestyle)
+        plt.semilogx(impact_parameters, magic, label=f'MAGIC, E={np.round(energy/1000, 3)} keV', linestyle=linestyle)
+        plt.gca().set_prop_cycle(None)
+
+        np.testing.assert_allclose(gm, gl, atol=5e-3) # 0.5% seems reasonable? max is ~0.3%
+        np.testing.assert_allclose(gm, mw, atol=5e-3)
+        np.testing.assert_allclose(mw, gl, atol=5e-3)
+        plt.legend()
+        plt.xlabel('p [A]')
+        plt.ylabel('theta [rad]')
+
 if show_plots: plt.show()
\ No newline at end of file

From ebfe6f738a300fd9b2fc3c32d897143b045352d3 Mon Sep 17 00:00:00 2001
From: Jon Drobny <drobny2@illinois.edu>
Date: Wed, 5 Aug 2026 16:54:50 -0700
Subject: [PATCH 10/12] used macros to streamline both lib.rs geometry types
 and main.rs geometry types

---
 src/lib.rs  | 84 +++++++++++++++++++++++++++++++++++++++++------------
 src/main.rs | 71 ++++++++++++++------------------------------
 2 files changed, 87 insertions(+), 68 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 51ba74c1..74503769 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -2200,42 +2200,90 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g
 
     Ok((theta_gm, theta_gl, theta_mw, theta_magic))
 }
-
 #[cfg(feature = "python")]
-#[pyfunction]
-#[pyo3(signature=(input, geometry_mode="1D"))]
-fn rustbca_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<()> {
-
-    match geometry_mode {
-        "1D" => {
-            let input: <Mesh1D as geometry::Geometry>::InputFileFormat = depythonize(&input).unwrap();
+macro_rules! geometry_typed_loops {
+    ($geometry_type:ty, $input:expr, $python:expr) => {
+        {
+            let input: <$geometry_type as geometry::Geometry>::InputFileFormat = depythonize(&$input).unwrap();
             let (particle_input_array, material, options, output_units) = input::process_input_file(input);
             let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build().unwrap();
             pool.install( ||
-                physics::physics_loop::<Mesh1D>(particle_input_array, material, options, output_units)
+                physics::physics_loop::<$geometry_type>(particle_input_array, material, options, output_units)
             );
             Ok(())
-        },
-       _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode)))
+        }
     }
 }
 
 #[cfg(feature = "python")]
 #[pyfunction]
 #[pyo3(signature=(input, geometry_mode="1D"))]
-fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<Bound<'py, PyAny>> {
-
+fn rustbca_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<()> {
     match geometry_mode {
-        "1D" => {
-            let input: <Mesh1D as geometry::Geometry>::InputFileFormat = depythonize(&input).unwrap();
+        "0D" => geometry_typed_loops!(Mesh0D, input, python),
+        "1D" => geometry_typed_loops!(Mesh1D, input, python),
+        "2D" => geometry_typed_loops!(Mesh2D, input, python),
+        "HOMOGENEOUS2D" => geometry_typed_loops!(Mesh2D, input, python),
+        "SPHERE" => geometry_typed_loops!(Sphere, input, python),
+        #[cfg(feature="parry3d")]
+        "BALL" => geometry_typed_loops!(ParryBall, input, python),
+        #[cfg(feature="parry3d")]
+        "TRIMESH" => geometry_typed_loops!(ParryTriMesh, input, python),
+       _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode)))
+    }
+}
+
+/*
+Notes on macros - this is the first I have written, so I'm taking notes here as I go.
+macro_rules! makes a macro - here, the macro is called geometry_types_silent_loops
+macros pattern match an argument and replace it with anything you want
+I want it to take a tuple of a string (e.g., "1D") and a type (e.g., Mesh1D)
+and plop those into corresponding match arms.
+The first line tells the macro to expect an argument with that pattern.
+arguments are $<name>:<designator>. Designators:
+block
+expr is used for expressions
+ident is used for variable/function names
+item
+literal is used for literal constants
+pat (pattern)
+path
+stmt (statement)
+tt (token tree)
+ty (type)
+vis (visibility qualifier)
+*/
+#[cfg(feature = "python")]
+macro_rules! geometry_typed_silent_loops {
+    ($geometry_type:ty, $input:expr, $python:expr) => {
+        {
+            let input: <$geometry_type as geometry::Geometry>::InputFileFormat = depythonize(&$input).unwrap();
             let (particle_input_array, material, options, output_units) = input::process_input_file(input);
             let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build().unwrap();
-            let finished_particles = pool.install(||
-                physics::silent_physics_loop::<Mesh1D>(particle_input_array, material, options, output_units.clone())
+            let finished_particles = pool.install( ||
+                physics::silent_physics_loop::<$geometry_type>(particle_input_array, material, options, output_units.clone())
             );
             let finished_particles_container = physics::process_finished_particles_to_arrays(finished_particles, output_units);
-            Ok(pythonize(python, &finished_particles_container)?)
+            Ok(pythonize($python, &finished_particles_container)?)
         }
+    }
+}
+
+#[cfg(feature = "python")]
+#[pyfunction]
+#[pyo3(signature=(input, geometry_mode="1D"))]
+fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<Bound<'py, PyAny>> {
+
+    match geometry_mode {
+        "0D" => geometry_typed_silent_loops!(Mesh0D, input, python),
+        "1D" => geometry_typed_silent_loops!(Mesh1D, input, python),
+        "2D" => geometry_typed_silent_loops!(Mesh2D, input, python),
+        "HOMOGENEOUS2D" => geometry_typed_silent_loops!(Mesh2D, input, python),
+        "SPHERE" => geometry_typed_silent_loops!(Sphere, input, python),
+        #[cfg(feature="parry3d")]
+        "BALL" => geometry_typed_silent_loops!(ParryBall, input, python),
+        #[cfg(feature="parry3d")]
+        "TRIMESH" => geometry_typed_silent_loops!(ParryTriMesh, input, python),
         _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode)))
     }
 }
\ No newline at end of file
diff --git a/src/main.rs b/src/main.rs
index 8a3944a0..baf4a924 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -69,6 +69,20 @@ pub use crate::math::duff_orthonormal_basis;
 #[cfg(feature = "parry3d")]
 pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, ParryTriMeshInput, InputParryTriMesh};
 
+
+macro_rules! main_loop {
+    ($geometry_type:ident, $input_file:expr) => {
+        {
+            let (particle_input_array, material, options, output_units) = input::input::<$geometry_type>($input_file);
+            //Initialize threads with rayon
+            println!("Processing {} ions...", particle_input_array.len());
+            println!("Initializing with {} threads...", options.num_threads);
+            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
+            physics_loop::<$geometry_type>(particle_input_array, material, options, output_units);
+        }
+    }
+}
+
 fn main() {
 
     let args: Vec<String> = env::args().collect();
@@ -92,57 +106,14 @@ fn main() {
     };
 
      match geometry_type {
-        GeometryType::MESH0D => {
-            let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh0D>(input_file);
-            //Initialize threads with rayon
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<Mesh0D>(particle_input_array, material, options, output_units);
-        },
-        GeometryType::MESH1D => {
-            let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh1D>(input_file);
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<Mesh1D>(particle_input_array, material, options, output_units);
-        },
-        GeometryType::MESH2D => {
-            let (particle_input_array, material, options, output_units) = input::input::<geometry::Mesh2D>(input_file);
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<Mesh2D>(particle_input_array, material, options, output_units);
-        },
-        GeometryType::SPHERE => {
-            let (particle_input_array, material, options, output_units) = input::input::<Sphere>(input_file);
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<Sphere>(particle_input_array, material, options, output_units);
-        },
+        GeometryType::MESH0D => main_loop!(Mesh0D, input_file),
+        GeometryType::MESH1D => main_loop!(Mesh1D, input_file),
+        GeometryType::MESH2D => main_loop!(Mesh2D, input_file),
+        GeometryType::SPHERE => main_loop!(Sphere, input_file),
         #[cfg(feature = "parry3d")]
-        GeometryType::BALL => {
-            let (particle_input_array, material, options, output_units) = input::input::<ParryBall>(input_file);
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<ParryBall>(particle_input_array, material, options, output_units);
-        }
+        GeometryType::BALL => main_loop!(ParryBall, input_file),
         #[cfg(feature = "parry3d")]
-        GeometryType::TRIMESH => {
-            let (particle_input_array, material, options, output_units) = input::input::<ParryTriMesh>(input_file);
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<ParryTriMesh>(particle_input_array, material, options, output_units);
-        }
-        GeometryType::HOMOGENEOUS2D => {
-            let (particle_input_array, material, options, output_units) = input::input::<geometry::HomogeneousMesh2D>(input_file);
-            println!("Processing {} ions...", particle_input_array.len());
-            println!("Initializing with {} threads...", options.num_threads);
-            let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
-            physics_loop::<HomogeneousMesh2D>(particle_input_array, material, options, output_units);
-        }
+        GeometryType::TRIMESH => main_loop!(ParryTriMesh, input_file),
+        GeometryType::HOMOGENEOUS2D => main_loop!(HomogeneousMesh2D, input_file),
     }
 }

From 2b04847cb8dc6a469d53ac394386fec2422b1a2d Mon Sep 17 00:00:00 2001
From: Jon Drobny <drobny2@illinois.edu>
Date: Wed, 5 Aug 2026 17:14:50 -0700
Subject: [PATCH 11/12] Ran cargo fix

---
 src/main.rs    | 4 ++--
 src/physics.rs | 1 -
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/main.rs b/src/main.rs
index baf4a924..8c11ce0e 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -6,7 +6,6 @@ use std::{env, fmt};
 use std::mem::discriminant;
 
 //Progress bar crate - works with rayon
-use indicatif::{ProgressBar, ProgressStyle};
 
 //Error handling crate
 use anyhow::{Result, Context, anyhow};
@@ -105,7 +104,8 @@ fn main() {
         _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 (<input file name>) or 2 (<geometry type> <input file name>)"),
     };
 
-     match geometry_type {
+    // This invokes the above macro that expands into the physics loop invocation for each type
+    match geometry_type {
         GeometryType::MESH0D => main_loop!(Mesh0D, input_file),
         GeometryType::MESH1D => main_loop!(Mesh1D, input_file),
         GeometryType::MESH2D => main_loop!(Mesh2D, input_file),
diff --git a/src/physics.rs b/src/physics.rs
index 443444c2..94e39d0b 100644
--- a/src/physics.rs
+++ b/src/physics.rs
@@ -1,5 +1,4 @@
 use super::*;
-use rayon::iter::*;
 use indicatif::{ProgressBar, ProgressStyle};
 
 pub fn silent_physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) -> Vec<particle::Particle> {

From 3cb7fad70939837c168369d67d4cbaf328bb5ca3 Mon Sep 17 00:00:00 2001
From: Jon Drobny <drobny2@illinois.edu>
Date: Wed, 5 Aug 2026 17:36:25 -0700
Subject: [PATCH 12/12] Attempt to fix module imports across lib, main

---
 src/lib.rs     | 6 +++---
 src/main.rs    | 6 ++----
 src/physics.rs | 2 +-
 3 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 74503769..67ebc4a6 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -9,10 +9,10 @@ use std::alloc::{dealloc, Layout};
 use std::mem::align_of;
 
 //Parallelization - currently only used in python library functions
+//#[cfg(feature = "python")]
+//use rayon::ThreadPoolBuilder;
 #[cfg(feature = "python")]
-use rayon::prelude::*;
-#[cfg(feature = "python")]
-use rayon::*;
+use rayon::iter::{IndexedParallelIterator, ParallelExtend, IntoParallelIterator, ParallelIterator};
 
 //Error handling crate
 use anyhow::{Result, Context, anyhow};
diff --git a/src/main.rs b/src/main.rs
index 8c11ce0e..09208127 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -5,8 +5,6 @@
 use std::{env, fmt};
 use std::mem::discriminant;
 
-//Progress bar crate - works with rayon
-
 //Error handling crate
 use anyhow::{Result, Context, anyhow};
 
@@ -18,8 +16,8 @@ use serde::*;
 use hdf5::*;
 
 //Parallelization
-use rayon::prelude::*;
-use rayon::*;
+//use rayon::prelude::*;
+//use rayon::ThreadPoolBuilder;
 
 //I/O
 use std::fs::OpenOptions;
diff --git a/src/physics.rs b/src/physics.rs
index 94e39d0b..177b911c 100644
--- a/src/physics.rs
+++ b/src/physics.rs
@@ -1,5 +1,6 @@
 use super::*;
 use indicatif::{ProgressBar, ProgressStyle};
+use rayon::iter::{IndexedParallelIterator, ParallelExtend, IntoParallelIterator, ParallelIterator};
 
 pub fn silent_physics_loop<T: Geometry + Sync>(particle_input_array: Vec<particle::ParticleInput>, material: material::Material<T>, options: Options, output_units: OutputUnits) -> Vec<particle::Particle> {
 
@@ -54,7 +55,6 @@ impl FinishedParticlesContainer {
     }
 }
 
-
 pub fn process_finished_particles_to_arrays(finished_particles: Vec<particle::Particle>, output_units: OutputUnits) -> FinishedParticlesContainer {
     let mut finished_particles_container = FinishedParticlesContainer::new();
     for particle in finished_particles {