diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..d6f7256
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,194 @@
+# Generated by Cargo
+# will have compiled files and executables
+debug/
+target/
+
+# These are backup files generated by rustfmt
+**/*.rs.bk
+
+# MSVC Windows builds of rustc generate these, which store debugging information
+*.pdb
+
+# RustRover
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# UV
+#   Similar to Pipfile.lock, it is generally recommended to include uv.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#uv.lock
+
+# poetry
+#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
+#poetry.lock
+
+# pdm
+#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
+#pdm.lock
+#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
+#   in version control.
+#   https://pdm.fming.dev/latest/usage/project/#working-with-version-control
+.pdm.toml
+.pdm-python
+.pdm-build/
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# PyCharm
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/
+
+# Ruff stuff:
+.ruff_cache/
+
+# PyPI configuration file
+.pypirc
+.DS_Store
+Cargo.lock
diff --git a/examples/test_rustbca.py b/examples/test_rustbca.py
index dd3d2b5..bb65f64 100644
--- a/examples/test_rustbca.py
+++ b/examples/test_rustbca.py
@@ -85,6 +85,15 @@ def main():
     print(f'Particle reflection coefficient for {ion["symbol"]} on {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {R_N}. Thomas predicts {np.round(thomas_reflection(ion, target, energy), 3)}.')
     print(f'Energy reflection coefficient for {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {R_E}')
 
+    S_N, R_N, R_E = compound_sputtering_reflection_coefficient(
+        ion, [target, ion], [target['n'], 0.1*target['n']], energy, angle, num_samples)
+    print(
+        f'Particle sputtering yield for {ion["symbol"]} on {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {S_N}. Yamamura  predicts {np.round(yamamura(ion, target, energy), 3)}.')
+    print(
+        f'Particle reflection coefficient for {ion["symbol"]} on {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {R_N}. Thomas predicts {np.round(thomas_reflection(ion, target, energy), 3)}.')
+    print(
+        f'Energy reflection coefficient for {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {R_E}')
+
     vx0 = 1e5
     vy0 = 1e5
     vz0 = 0.0
diff --git a/src/lib.rs b/src/lib.rs
index 3a9f53e..823e785 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -94,6 +94,7 @@ pub fn libRustBCA(py: Python, m: &PyModule) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(sputtering_yield, m)?)?;
     m.add_function(wrap_pyfunction!(reflection_coefficient, m)?)?;
     m.add_function(wrap_pyfunction!(compound_reflection_coefficient, m)?)?;
+    m.add_function(wrap_pyfunction!(compound_sputtering_reflection_coefficient, m)?)?;
     m.add_function(wrap_pyfunction!(reflect_single_ion_py, m)?)?;
     #[cfg(feature = "parry3d")]
     m.add_function(wrap_pyfunction!(rotate_given_surface_normal_py, m)?)?;
@@ -106,6 +107,15 @@ pub fn libRustBCA(py: Python, m: &PyModule) -> PyResult<()> {
     Ok(())
 }
 
+pub fn find_index(Z: f64, Z2: &[f64]) -> Option<usize> {
+                for (index, &value) in Z2.iter().enumerate() {
+                    if value == Z {
+                return Some(index);
+            }
+        }
+        None
+    }
+
 #[derive(Debug)]
 #[repr(C)]
 pub struct InputSimpleBCA {
@@ -1948,7 +1958,7 @@ pub fn compound_reflection_coefficient(ion: &PyDict, targets: Vec<&PyDict>, targ
         m: m2,
         interaction_index: vec![0; number_target_species],
         surface_binding_model: SurfaceBindingModel::AVERAGE,
-        bulk_binding_model: BulkBindingModel::INDIVIDUAL,
+        bulk_binding_model: BulkBindingModel::AVERAGE,//INDIVIDUAL,
     };
 
     let geometry_input = geometry::Mesh0DInput {
@@ -1994,3 +2004,126 @@ pub fn compound_reflection_coefficient(ion: &PyDict, targets: Vec<&PyDict>, targ
 
     (num_reflected as f64 / num_samples as f64, energy_reflected / EV / energy / num_samples as f64)
 }
+
+
+#[cfg(feature = "python")]
+#[pyfunction]
+/// compound_sputtering_reflection_coefficient(ion, target_species, target_number_densities, energy, angle, num_samples)
+/// A routine the calculates the reflection coefficient of energetic ions incident upon materials using RustBCA.
+/// Args:
+///     ion: a dictionary with the keys Z (atomic number), m (atomic mass in AMU), Ec (cutoff energy in eV), Es (surface binding energy in eV)
+///     target_species: a list of dictionaries with the keys Z, m, Ec, Es, Eb (bulk binding energy in eV), n2 (number density in 1/m3)
+///     target_number_densities (list(f64)): number density of each target species in the compound
+///     energy: the incident energy of the ion in eV
+///     angle: incident angle of the ion in degrees from surface normal
+///     num_samples: number of ion trajectories to run; precision will go as 1/sqrt(N)
+/// Returns:
+///     S_N (f64): sputtering yield (number of particles sputtered / number of incident particles)
+///     R_N (f64): reflection coefficient (number of particles reflected / number of incident particles)
+///     R_E (f64): energy reflection coefficient (sum of reflected particle energies / total incident energy)
+pub fn compound_sputtering_reflection_coefficient(ion: &PyDict, targets: Vec<&PyDict>, target_number_densities: Vec<f64>, energy: f64, angle: f64, num_samples: usize) -> (Vec<f64>, f64, f64) {
+
+    const DELTA: f64 = 1e-6;
+
+    assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle);
+
+    let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists."));
+    let m1 = unpack(ion.get_item("m").expect("Cannot get ion mass from dictionary. Ensure ion['m'] exists."));
+    let Ec1 = unpack(ion.get_item("Ec").expect("Cannot get ion cutoff energy from dictionary. Ensure ion['Ec'] exists."));
+    let Es1 = unpack(ion.get_item("Es").expect("Cannot get ion surface binding energy from dictionary. Ensure ion['Es'] exists."));
+
+    let Z2: Vec<f64> = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Z").unwrap_or_else(|| panic!("Cannot get target Z from dictionary at index {}. Ensure target['Z'] exists.", index)))).collect();
+    let m2: Vec<f64> = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("m").unwrap_or_else(|| panic!("Cannot get target m from dictionary at index {}. Ensure target['m'] exists.", index)))).collect();
+    let Ec2: Vec<f64> = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Ec").unwrap_or_else(|| panic!("Cannot get target Ec from dictionary at index {}. Ensure target['Ec'] exists.", index)))).collect();
+    let Es2: Vec<f64> = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Es").unwrap_or_else(|| panic!("Cannot get target Es from dictionary at index {}. Ensure target['Es'] exists.", index)))).collect();
+    let Eb2: Vec<f64> = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Eb").unwrap_or_else(|| panic!("Cannot get target Eb from dictionary at index {}. Ensure target['Eb'] exists.", index)))).collect();
+    let Z2_ = Z2.clone(); 
+    let number_target_species = Z2.len();
+
+    let options = Options::default_options(true);
+
+    let y = 0.0;
+    let z = 0.0;
+
+    let ux = (angle/180.0*PI).cos() + DELTA;
+    let uy = (angle/180.0*PI).sin() - DELTA;
+    let uz = 0.0;
+
+    let mut direction = Vector::new(ux, uy, uz);
+    direction.normalize();
+
+    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; number_target_species],
+        Z: Z2,
+        m: m2,
+        interaction_index: vec![0; number_target_species],
+        surface_binding_model: SurfaceBindingModel::AVERAGE,
+        bulk_binding_model: BulkBindingModel::INDIVIDUAL,
+    };
+
+    let geometry_input = geometry::Mesh0DInput {
+        length_unit: "M".to_string(),
+        densities: target_number_densities,
+        electronic_stopping_correction_factor: 1.0
+    };
+
+    let m = material::Material::<Mesh0D>::new(&material_parameters, &geometry_input);
+
+    let x = -m.geometry.energy_barrier_thickness;
+
+    let num_reflected = Mutex::new(0);
+    let energy_reflected = Mutex::new(0.0);
+    let num_sputtered: Mutex<Vec<f64>> = Mutex::new(vec![0.0; targets.len()]);
+
+
+    (0..num_samples as u64).into_par_iter().for_each( |index| {
+
+        let p = particle::Particle::default_incident(
+            m1,
+            Z1,
+            energy,
+            Ec1,
+            Es1,
+            x,
+            ux,
+            uy,
+            uz
+        );
+
+        let output = bca::single_ion_bca(p, &m, &options);
+
+        for particle in output {
+            //println!("output");
+            if particle.E > 0.0 && particle.dir.x < 0.0 && particle.left {
+                if (particle.incident) {
+                    //println!("reflected");
+                    let mut num_reflected = num_reflected.lock().unwrap();
+                    *num_reflected += 1;
+                    let mut energy_reflected = energy_reflected.lock().unwrap();
+                    *energy_reflected += particle.E;
+                }
+                else{
+                    //println!("sputtered");
+                    let mut num_sputtered = num_sputtered.lock().unwrap();
+                    let ind = find_index(particle.Z, &Z2_).unwrap();
+
+                    num_sputtered[ind] += 1.0;
+                    //println!("num_sputtered: {:?}", *num_sputtered);
+                }
+            }
+        }
+    });
+
+
+
+    let num_reflected = *num_reflected.lock().unwrap();
+    let energy_reflected = *energy_reflected.lock().unwrap();
+    let num_sputtered = num_sputtered.lock().unwrap();
+    let num_sputtered_: Vec<f64> = num_sputtered.iter().map(|&x| x / num_samples as f64).collect();
+    (num_sputtered_, num_reflected as f64 / num_samples as f64, energy_reflected / EV / energy / num_samples as f64)
+}