diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml
index 584476c..382c9b8 100644
--- a/.github/workflows/rustbca_compile_check.yml
+++ b/.github/workflows/rustbca_compile_check.yml
@@ -59,12 +59,12 @@ jobs:
         cargo run --release 0D examples/boron_nitride_0D.toml
         ./target/release/RustBCA 0D examples/titanium_dioxide_0D.toml
         ./target/release/RustBCA 1D examples/layered_geometry_1D.toml
-        cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
+        echo "layered geometry 1D:"; cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
         ./target/release/RustBCA examples/layered_geometry.toml
-        cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
+        echo "layered geometry 2D:"; cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
         ./target/release/RustBCA SPHERE examples/boron_nitride_sphere.toml
         cargo run --release --features parry3d TRIMESH examples/tungsten_twist_trimesh.toml
         ./target/release/RustBCA examples/boron_nitride_wire.toml
-        cat boron_nitride_summary.output
+        echo "boron nitride 2D:"; cat boron_nitride_summary.output
         ./target/release/RustBCA HOMOGENEOUS2D examples/boron_nitride_wire_homogeneous.toml
-        cat boron_nitride_h_summary.output
+        echo "boron nitride 2D Homogeneous:"; cat boron_nitride_h_summary.output
diff --git a/Cargo.toml b/Cargo.toml
index a335ee2..c945360 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -5,6 +5,8 @@ default-run = "RustBCA"
 authors = ["Jon Drobny <drobny2@illinois.edu>", "Jon Drobny <jdrobny@tae.com>"]
 edition = "2021"
 
+
+
 [[bin]]
 name = "RustBCA"
 path = "src/main.rs"
@@ -28,6 +30,7 @@ hdf5 = {version = "0.7.1", optional = true}
 rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
 ndarray = {version = "0.14.0", features = ["serde"], optional = true}
 parry3d-f64 = {optional = true, version="0.2.0"}
+parry2d-f64 = {optional = false, version="0.17.0"}
 
 [dependencies.pyo3]
 version = "0.19.0"
@@ -41,7 +44,7 @@ float-cmp = "0.8.0"
 lto = "fat"
 codegen-units = 1
 opt-level = 3
-debug = false
+debug = 1
 
 [features]
 hdf5_input = ["hdf5"]
diff --git a/README.md b/README.md
index 6e4223a..d5b5ee9 100644
--- a/README.md
+++ b/README.md
@@ -12,9 +12,9 @@ By discretizing the collision cascade into a sequence of binary collisions,
 between an energetic ion and a target material. This includes reflection,
 implantation, and transmission of the incident ion, as well as sputtering
 and displacement damage of the target. Generally, [BCA] codes can be
-valid for incident ion energies between approximately ~1 eV/nucleon 
+valid for incident ion energies between approximately ~1 eV/nucleon
 to <1 GeV/nucleon. Improvements to RustBCA have expanded the regime
-of validity for some quantities, such as reflection coefficients, below 
+of validity for some quantities, such as reflection coefficients, below
 1 eV/nucleon.
 
 Check out the `RustBCA` [Wiki] for detailed information, installation
@@ -160,7 +160,7 @@ Windows, MacOS, and Linux systems.
 * [rcpr], a CPR and polynomial rootfinder, required for using attractive-repulsive interaction potentials such as Lennard-Jones or Morse.
 * For manipulating input files and running associated scripts, the following are required:
   * [Python] 3.6+
-  * The [Python] libraries: `numpy`, `matplotlib`, `toml` (must build from source), `shapely`, and `scipy`.
+  * The [Python] libraries: `numpy`, `matplotlib`, `toml`, `shapely`, and `scipy`.
 
 ### Detailed instructions for Ubuntu 18.04 LTS
 
@@ -187,36 +187,28 @@ git clone https://github.com/uiri/toml.git
 cd toml
 python3 setup.py install
 ```
-7. (Optional) Install software for [rcpr]:
-```bash
-sudo apt-get install gcc gfortran build-essential cmake liblapack-dev libblas-dev liblapacke-dev
-```
-8. (Optional - should come with rustup) Install `cargo`:
+7. (Optional - should come with rustup) Install `cargo`:
 ```bash
 sudo apt-get install cargo
 ```
-9. Build `RustBCA`:
+8. Build `RustBCA`:
 ```bash
 git clone https://github.com/lcpp-org/rustBCA
 cd RustBCA
 cargo build --release
 ```
-10. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr` (with your choice of backend: `openblas`, `netlib`, or `intel-mkl`):
+9. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr`:
 ```bash
-cargo build --release --features cpr_rootfinder_netlib,hdf5_input
-cargo build --release --features cpr_rootfinder_openblas,hdf5_input
-cargo build --release --features cpr_rootfinder_intel_mkl,hdf5_input
+cargo build --release --features cpr_rootfinder,hdf5_input
  ```
-11. `input.toml` is the input file - see [Usage](https://github.com/lcpp-org/RustBCA/wiki/Usage,-Input-File,-and-Output-Files) for more information
-12. Run the required tests using:
+10. `input.toml` is the input file - see [Usage](https://github.com/lcpp-org/RustBCA/wiki/Usage,-Input-File,-and-Output-Files) for more information
+11. Run the required tests using:
 ```bash
 cargo test
 ```
-13. (Optional) Run the required and optional tests for the desired backend(s):
+12. (Optional) Run the required and optional tests for the desired backend(s):
 ```bash
-cargo test --features cpr_rootfinder_netlib
-cargo test --features cpr_rootfinder_openblas
-cargo test --features cpr_rootfinder_intel_mkl
+cargo test --features cpr_rootfinder
 ```
 
 ### Detailed instructions for Fedora 33
@@ -275,7 +267,7 @@ command line argument.
 ./RustBCA /path/to/input.toml
 ```
 
-Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `TRIMESH`, `SPHERE` - see the wiki for more details):
+Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `TRIMESH`, `SPHERE`, `HOMOGENEOUS2D` - see the wiki for more details):
 ```bash
 ./RustBCA 0D /path/to/input.toml
 ```
diff --git a/examples/boron_nitride_wire_homogeneous.toml b/examples/boron_nitride_wire_homogeneous.toml
index 57156b7..9afd461 100644
--- a/examples/boron_nitride_wire_homogeneous.toml
+++ b/examples/boron_nitride_wire_homogeneous.toml
@@ -1,4 +1,4 @@
-#Use with 0D geometry option only
+#Use with HOMOGENEOUS geometry option only
 [options]
 name = "boron_nitride_h_"
 weak_collision_order = 0
diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml
index 235cc5a..078c799 100644
--- a/examples/layered_geometry.toml
+++ b/examples/layered_geometry.toml
@@ -1,7 +1,7 @@
 #Use with 2D geometry option only
 [options]
 name = "2000.0eV_0.0001deg_He_TiO2_Al_Si"
-track_trajectories = false
+track_trajectories = true
 track_recoils = true
 track_recoil_trajectories = false
 write_buffer_size = 8000
@@ -54,7 +54,8 @@ points = [[0.0, -0.5], [0.01, -0.5], [0.04, -0.5], [0.5, -0.5], [0.5, 0.5], [0.0
 triangles = [[0, 1, 6], [0, 6, 7], [1, 2, 5], [1, 5, 6], [2, 3, 4], [2, 4, 5]]
 densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],]
 boundary = [0, 3, 4, 7]
-simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],]
+#simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],]
+simulation_boundary_points = [[0.6, 0.6], [-0.1, 0.6], [-0.1, -0.6], [0.6, -0.6]]
 electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
 energy_barrier_thickness = 2.2E-4
 
diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml
index f2194ab..736b39d 100644
--- a/examples/multiple_interaction_potentials.toml
+++ b/examples/multiple_interaction_potentials.toml
@@ -48,7 +48,7 @@ length_unit = "MICRON"
 triangles = [ [ 0.0, 0.01, 0.0, 0.5, 0.5, -0.5,], [ 0.0, 0.01, 0.01, -0.5, 0.5, -0.5,], [ 0.01, 0.01, 0.04, -0.5, 0.5, -0.5,], [ 0.01, 0.04, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.5, -0.5, 0.5, -0.5,],]
 densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],]
 material_boundary_points = [ [ 0.0, 0.5,], [ 0.5, 0.5,], [ 0.5, -0.5,], [ 0.0, -0.5,],]
-simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],]
+simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],]
 electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
 
 [material_parameters]
diff --git a/src/enums.rs b/src/enums.rs
index 908c040..d8582b6 100644
--- a/src/enums.rs
+++ b/src/enums.rs
@@ -3,13 +3,13 @@ use super::*;
 pub enum MaterialType {
     MESH0D(material::Material<geometry::Mesh0D>),
     MESH1D(material::Material<geometry::Mesh1D>),
-    MESH2D(material::Material<geometry::Mesh2D>),
+    MESH2D(material::Material<geometry::ParryMesh2D>),
     SPHERE(material::Material<sphere::Sphere>),
     #[cfg(feature = "parry3d")]
     BALL(material::Material<parry::ParryBall>),
     #[cfg(feature = "parry3d")]
     TRIMESH(material::Material<parry::ParryTriMesh>),
-    HOMOGENEOUS2D(material::Material<geometry::HomogeneousMesh2D>)
+    HOMOGENEOUS2D(material::Material<geometry::ParryHomogeneousMesh2D>)
 }
 
 #[derive(Deserialize, PartialEq, Clone, Copy)]
diff --git a/src/geometry.rs b/src/geometry.rs
index a915629..814f7e5 100644
--- a/src/geometry.rs
+++ b/src/geometry.rs
@@ -1,8 +1,11 @@
 use super::*;
 
-use geo::algorithm::contains::Contains;
-use geo::{Polygon, LineString, Point, point, Closest};
-use geo::algorithm::closest_point::ClosestPoint;
+use itertools::Itertools;
+
+use parry2d_f64::shape::{Polyline, TriMesh, FeatureId::*};
+use parry2d_f64::query::PointQuery;
+use parry2d_f64::math::Point as Point2d;
+use parry2d_f64::math::Isometry;
 
 ///Trait for a Geometry object - all forms of geometry must implement these traits to be used
 pub trait Geometry {
@@ -18,6 +21,7 @@ pub trait Geometry {
     fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool;
     fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool;
     fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64);
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64);
 }
 
 pub trait GeometryElement: Clone {
@@ -109,6 +113,10 @@ impl Geometry for Mesh0D {
     fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
         (0., y, z)
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        (-1., 0., 0.)
+    }
 }
 
 #[derive(Deserialize, Clone)]
@@ -265,6 +273,175 @@ impl Geometry for Mesh1D {
             (self.top, y, z)
         }
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        (-1., 0., 0.)
+    }
+}
+
+#[derive(Deserialize, Clone)]
+pub struct ParryHomogeneousMesh2DInput {
+    pub length_unit: String,
+    pub points: Vec<(f64, f64)>,
+    pub simulation_boundary_points: Vec<(f64, f64)>,
+    pub densities: Vec<f64>,
+    pub electronic_stopping_correction_factor: f64
+}
+
+#[derive(Clone)]
+pub struct ParryHomogeneousMesh2D {
+    pub boundary: Polyline,
+    pub simulation_boundary: Polyline,
+    pub energy_barrier_thickness: f64,
+    pub densities: Vec<f64>,
+    pub electronic_stopping_correction_factor: f64,
+    pub concentrations: Vec<f64>
+}
+
+
+impl Geometry for ParryHomogeneousMesh2D {
+
+    type InputFileFormat = InputHomogeneous2D;
+
+    fn new(input: &<<Self as Geometry>::InputFileFormat as GeometryInput>::GeometryInput) -> Self {
+        let length_unit: f64 = match input.length_unit.as_str() {
+            "MICRON" => MICRON,
+            "CM" => CM,
+            "MM" => MM,
+            "ANGSTROM" => ANGSTROM,
+            "NM" => NM,
+            "M" => 1.,
+            _ => input.length_unit.parse()
+                .expect(format!(
+                        "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
+                        &input.length_unit.as_str()
+                    ).as_str()),
+        };
+
+        let mut boundary_points_converted: Vec<Point2d<f64>> = input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect();
+        let number_boundary_points = boundary_points_converted.len() as u32;
+
+        for p in boundary_points_converted.iter().combinations(2) {
+            assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).")
+        }
+
+        let mut simulation_boundary_points_converted: Vec<Point2d<f64>> = input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect();
+        let number_simulation_boundary_points = simulation_boundary_points_converted.len() as u32;
+
+        let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize)
+        .map(|i| (simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted[i].x)*(simulation_boundary_points_converted[i].y + simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].y))
+        .sum::<f64>();
+
+        if test_simulation_boundary_ccw > 0.0 {
+            simulation_boundary_points_converted = simulation_boundary_points_converted.clone().into_iter().rev().collect::<Vec<Point2d<f64>>>();
+        }
+
+        for p in simulation_boundary_points_converted.iter().combinations(2) {
+            assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).")
+        }
+
+        let test_ccw = (0..number_boundary_points as usize)
+                .map(|i| (boundary_points_converted[(i + 1) % number_boundary_points as usize].x - boundary_points_converted[i].x)*(boundary_points_converted[i].y + boundary_points_converted[(i + 1) % number_boundary_points as usize].y))
+                .sum::<f64>();
+
+        if test_ccw > 0.0 {
+            boundary_points_converted = boundary_points_converted.clone().into_iter().rev().collect::<Vec<Point2d<f64>>>();
+        }
+
+        let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor;
+
+        let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
+
+        let total_density: f64 = densities.iter().sum();
+
+        let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
+
+        let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
+
+        let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::<Vec<[u32; 2]>>();
+        linked_boundary_points.push([number_boundary_points - 1, 0]);
+        let boundary2 = Polyline::new(boundary_points_converted, Some(linked_boundary_points));
+
+        let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::<Vec<[u32; 2]>>();
+        linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]);
+        let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points));
+
+        ParryHomogeneousMesh2D {
+            densities,
+            simulation_boundary: simulation_boundary2,
+            boundary: boundary2,
+            electronic_stopping_correction_factor,
+            energy_barrier_thickness,
+            concentrations
+        }
+    }
+
+    fn get_densities(&self,  x: f64, y: f64, z: f64) -> &Vec<f64> {
+        &self.densities
+    }
+
+    fn get_ck(&self,  x: f64, y: f64, z: f64) -> f64 {
+        self.electronic_stopping_correction_factor
+    }
+
+    fn get_total_density(&self,  x: f64, y: f64, z: f64) -> f64 {
+        self.densities.iter().sum::<f64>()
+    }
+
+    fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec<f64> {
+        &self.concentrations
+    }
+
+    fn inside(&self, x: f64, y: f64, z: f64) -> bool {
+        let p = Point2d::new(x, y);
+        if self.boundary.aabb(&Isometry::identity()).contains_local_point(&p) {
+            let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p);
+            point_projection.is_inside
+        } else {
+            false
+        }
+    }
+
+    fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool {
+        let p = Point2d::new(x, y);
+
+        let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p);
+        point_projection.is_inside
+
+    }
+    fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool {
+        if self.inside(x, y, z) {
+            true
+        } else {
+            let p = Point2d::new(x, y);
+            (self.boundary.distance_to_local_point(&p, true) as f64) < self.energy_barrier_thickness
+        }
+    }
+
+    fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let p = Point2d::new(x, y);
+        let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p);
+        let (x_, y_) = (point_projection.point.x, point_projection.point.y);
+        (x_ as f64, y_ as f64, z)
+    }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let p = Point2d::new(x, y);
+
+        let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p);
+        let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z);
+
+        let dx = x_intersect - x;
+        let dy = y_intersect - y;
+        let dz = z_intersect - z;
+        let mag = (dx*dx + dy*dy + dz*dz).sqrt();
+
+        if point_projection.is_inside {
+            (dx/mag, dy/mag, dz/mag)
+        } else {
+            (-dx/mag, -dy/mag, -dz/mag)
+        }
+    }
 }
 
 #[derive(Deserialize, Clone)]
@@ -276,6 +453,7 @@ pub struct HomogeneousMesh2DInput {
     pub electronic_stopping_correction_factor: f64
 }
 
+/*
 #[derive(Clone)]
 pub struct HomogeneousMesh2D {
     pub boundary: Polygon<f64>,
@@ -379,7 +557,12 @@ impl Geometry for HomogeneousMesh2D {
             panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
         }
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        panic!("Not implemented.")
+    }
 }
+*/
 
 /// Object that contains raw mesh input data.
 #[derive(Deserialize, Clone)]
@@ -394,6 +577,209 @@ pub struct Mesh2DInput {
     pub energy_barrier_thickness: f64,
 }
 
+#[derive(Clone)]
+pub struct ParryMesh2D {
+    trimesh: TriMesh,
+    densities: Vec<Vec<f64>>,
+    electronic_stopping_correction_factors: Vec<f64>,
+    concentrations: Vec<Vec<f64>>,
+    pub boundary: Polyline,
+    pub simulation_boundary: Polyline,
+    pub energy_barrier_thickness: f64
+}
+
+impl ParryMesh2D {
+    fn get_id(&self, x: f64, y: f64, z: f64) -> usize {
+        let p = Point2d::new(x, y);
+        let (point_projection, feature_id) = self.trimesh.project_local_point_and_get_feature(&p);
+        match feature_id {
+            Vertex(vertex_id) => {
+                for (triangle_id, triangle) in self.trimesh.triangles().enumerate() {
+                    if triangle.contains_local_point(&self.trimesh.vertices()[vertex_id as usize]) {
+                        return triangle_id as usize
+                    }
+                }
+                panic!("Geometry error: point ({}, {}) located on a vertex that is not associated with any triangle. Check mesh input.", x, y);
+            },
+            Face(triangle_id) => {
+                triangle_id as usize
+            }
+            Unknown => {
+                panic!("Geometry error: unknown feature detected near ({}, {}). Check mesh input.", x, y);
+            }
+        }
+    }
+}
+
+impl Geometry for ParryMesh2D {
+    type InputFileFormat = Input2D;
+
+    fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool {
+        if self.inside(x, y, z) {
+            true
+        } else {
+            let p = Point2d::new(x, y);
+            if self.trimesh.distance_to_local_point(&p, true) < self.energy_barrier_thickness {
+                true
+            } else {
+                false
+            }
+        }
+    }
+
+    fn inside_simulation_boundary (&self, x: f64, y: f64, z: f64) -> bool {
+        let p = Point2d::new(x, y);
+
+        let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p);
+
+        point_projection.is_inside
+    }
+
+    fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let p = Point2d::new(x, y);
+        let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p);
+        let (x_, y_) = (point_projection.point.x, point_projection.point.y);
+        (x_ as f64, y_ as f64, z)
+    }
+
+    fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec<f64> {
+        &self.densities[self.get_id(x, y, z)]
+    }
+
+    fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 {
+        self.electronic_stopping_correction_factors[self.get_id(x, y, z)]
+    }
+
+    /// Determine the total number density of the triangle that contains or is nearest to (x, y).
+    fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 {
+        self.get_densities(x, y, z).iter().sum::<f64>()
+    }
+
+    /// Find the concentrations of the triangle that contains or is nearest to (x, y).
+    fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec<f64> {
+        &self.concentrations[self.get_id(x, y, z)]
+    }
+
+    /// Determines whether the point (x, y) is inside the mesh.
+    fn inside(&self, x: f64, y: f64, z: f64) -> bool {
+        let p = Point2d::new(x, y);
+        if self.boundary.aabb(&Isometry::identity()).contains_local_point(&p) {
+            let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p);
+            point_projection.is_inside
+        } else {
+            false
+        }
+    }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let p = Point2d::new(x, y);
+
+        let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p);
+        let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z);
+
+        let dx = x_intersect - x;
+        let dy = y_intersect - y;
+        let dz = z_intersect - z;
+        let mag = (dx*dx + dy*dy + dz*dz).sqrt();
+
+        if point_projection.is_inside {
+            (dx/mag, dy/mag, dz/mag)
+        } else {
+            (-dx/mag, -dy/mag, -dz/mag)
+        }
+    }
+
+    fn new(geometry_input: &<<Self as Geometry>::InputFileFormat as GeometryInput>::GeometryInput) -> ParryMesh2D {
+        let length_unit: f64 = match geometry_input.length_unit.as_str() {
+            "MICRON" => MICRON,
+            "CM" => CM,
+            "MM" => MM,
+            "ANGSTROM" => ANGSTROM,
+            "NM" => NM,
+            "M" => 1.,
+            _ => geometry_input.length_unit.parse()
+                .expect(format!(
+                        "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
+                        &geometry_input.length_unit.as_str()
+                    ).as_str()),
+        };
+
+        let mut boundary_points_converted: Vec<Point2d<f64>> = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect();
+        let number_boundary_points = boundary_points_converted.len() as u32;
+
+        let mut simulation_boundary_points_converted: Vec<Point2d<f64>> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect();
+        let number_simulation_boundary_points = simulation_boundary_points_converted.len() as u32;
+
+        let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize)
+        .map(|i| (simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted[i].x)*(simulation_boundary_points_converted[i].y + simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].y))
+        .sum::<f64>();
+
+        dbg!(test_simulation_boundary_ccw > 0.0);
+
+        dbg!(&simulation_boundary_points_converted);
+
+        if test_simulation_boundary_ccw > 0.0 {
+            simulation_boundary_points_converted = simulation_boundary_points_converted.clone().into_iter().rev().collect::<Vec<Point2d<f64>>>();
+        }
+
+        dbg!(&simulation_boundary_points_converted);
+
+        for p in simulation_boundary_points_converted.iter().combinations(2) {
+            assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).")
+        }
+
+        let test_ccw = (0..number_boundary_points as usize)
+                .map(|i| (boundary_points_converted[(i + 1) % number_boundary_points as usize].x - boundary_points_converted[i].x)*(boundary_points_converted[i].y + boundary_points_converted[(i + 1) % number_boundary_points as usize].y))
+                .sum::<f64>();
+
+        if test_ccw > 0.0 {
+            boundary_points_converted = boundary_points_converted.clone().into_iter().rev().collect::<Vec<Point2d<f64>>>();
+        }
+
+        for p in boundary_points_converted.iter().combinations(2) {
+            assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).")
+        }
+
+        let triangles = geometry_input.triangles.iter().map(|(i, j, k)| [*i as u32, *j as u32, *k as u32]).collect::<Vec<[u32; 3]>>();
+        let points_converted = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect::<Vec<Point2d<f64>>>();
+
+        let trimesh = TriMesh::new(points_converted, triangles);
+
+        let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone();
+
+        let densities: Vec<Vec<f64>> = geometry_input.densities.iter().map(|density_list| density_list.iter().map(|density| density / (length_unit).powi(3)).collect::<Vec<f64>>()).collect::<Vec<Vec<f64>>>();
+
+
+        let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit;
+
+        let concentrations: Vec<Vec<f64>> = densities.iter().map(|density_list| density_list.iter().map(|density| density / density_list.iter().sum::<f64>() ).collect::<Vec<f64>>()).collect::<Vec<Vec<f64>>>();
+
+        for concentration in &concentrations {
+            assert!((concentration.iter().sum::<f64>() - 1.0).abs() < 1e-6);
+        }
+
+        let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::<Vec<[u32; 2]>>();
+        linked_boundary_points.push([number_boundary_points - 1, 0]);
+        let boundary2 = Polyline::new(boundary_points_converted, Some(linked_boundary_points));
+
+        let number_simulation_boundary_points = simulation_boundary_points_converted.clone().len() as u32;
+        let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::<Vec<[u32; 2]>>();
+        linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]);
+        let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points));
+
+        ParryMesh2D {
+            trimesh,
+            densities,
+            simulation_boundary: simulation_boundary2,
+            boundary: boundary2,
+            electronic_stopping_correction_factors,
+            energy_barrier_thickness,
+            concentrations
+        }
+    }
+}
+
+/*
 /// Triangular mesh for rustbca.
 #[derive(Clone)]
 pub struct Mesh2D {
@@ -491,26 +877,6 @@ impl Geometry for Mesh2D {
             cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck));
         }
 
-        /*
-        for ((coordinate_set, densities), ck) in triangles.iter().zip(densities).zip(electronic_stopping_correction_factors) {
-            let coordinate_set_converted = (
-                coordinate_set.0*length_unit,
-                coordinate_set.1*length_unit,
-                coordinate_set.2*length_unit,
-                coordinate_set.3*length_unit,
-                coordinate_set.4*length_unit,
-                coordinate_set.5*length_unit,
-            );
-
-            let total_density: f64 = densities.iter().sum();
-            let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
-
-            cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck));
-        }
-        */
-
-
-
         let mut boundary_points_converted = Vec::with_capacity(material_boundary_point_indices.len());
         for index in material_boundary_point_indices.iter() {
             boundary_points_converted.push((points[*index].0*length_unit, points[*index].1*length_unit));
@@ -549,8 +915,18 @@ impl Geometry for Mesh2D {
 
     fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
         if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+            if p.x() == x && p.y() == y {
+                dbg!(&self.boundary);
+                println!("single point");
+                panic!("({} {} {}) ({} {} {})", p.x(), p.y(), z, x, y, z);
+            }
             (p.x(), p.y(), z)
         } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+            if p.x() == x && p.y() == y {
+                dbg!(&self.boundary);
+                println!("intersection");
+                panic!("({} {} {}) ({} {} {})", p.x(), p.y(), z, x, y, z);
+            }
             (p.x(), p.y(), z)
         } else {
             panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
@@ -617,6 +993,10 @@ impl Geometry for Mesh2D {
     fn inside(&self, x: f64, y: f64, z: f64) -> bool {
         self.boundary.contains(&Point::new(x, y))
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        panic!("Not implemented.");
+    }
 }
 
 /// A mesh cell that contains a triangle and the local number densities and concentrations.
@@ -662,6 +1042,7 @@ impl GeometryElement for Cell2D {
         self.electronic_stopping_correction_factor
     }
 }
+*/
 
 #[derive(Clone)]
 pub struct Layer1D {
@@ -711,6 +1092,7 @@ impl GeometryElement for Layer1D {
     }
 }
 
+/*
 /// A triangle in 2D, with points (x1, y1), (x2, y2), (x3, y3), and the three line segments bewtween them.
 #[derive(Clone)]
 pub struct Triangle2D {
@@ -801,3 +1183,4 @@ impl Triangle2D {
         ((x - centroid.0)*(x - centroid.0) + (y - centroid.1)*(y - centroid.1)).sqrt()
     }
 }
+*/
\ No newline at end of file
diff --git a/src/input.rs b/src/input.rs
index 8e2b81d..af04e1f 100644
--- a/src/input.rs
+++ b/src/input.rs
@@ -172,11 +172,6 @@ fn one_u64() -> u64 {
     1
 }
 
-///This helper function is a workaround to issue #368 in serde
-fn three() -> usize {
-    3
-}
-
 fn zero_usize() -> usize{
     0
 }
diff --git a/src/lib.rs b/src/lib.rs
index 3a9f53e..595ef24 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -56,6 +56,7 @@ use pyo3::types::*;
 //Load internal modules
 pub mod material;
 pub mod particle;
+#[cfg(test)]
 pub mod tests;
 pub mod interactions;
 pub mod bca;
@@ -75,7 +76,7 @@ pub use crate::consts::*;
 pub use crate::structs::*;
 pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput};
 pub use crate::output::{OutputUnits};
-pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D};
+pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, ParryMesh2D};
 pub use crate::sphere::{Sphere, SphereInput, InputSphere};
 
 #[cfg(feature = "parry3d")]
diff --git a/src/main.rs b/src/main.rs
index 7da24c6..a66f4d7 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -58,7 +58,7 @@ pub use crate::consts::*;
 pub use crate::structs::*;
 pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput};
 pub use crate::output::{OutputUnits};
-pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, HomogeneousMesh2D};
+pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, ParryMesh2D, ParryHomogeneousMesh2D};
 pub use crate::sphere::{Sphere, SphereInput, InputSphere};
 pub use crate::physics::{physics_loop};
 
@@ -84,7 +84,7 @@ fn main() {
             "HOMOGENEOUS2D" => GeometryType::HOMOGENEOUS2D,
             _ => panic!("Unimplemented geometry {}.", args[1].clone())
         }),
-        _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 (<input file name>) or 2 (<geometry type> <input file name>)"),
+        _ => panic!("Too many command line arguments. RustBCA accepts 0 (defaulting to 'input.toml' and 2D mesh), 1 (<input file name> defaulting to 2D mesh), or 2 (<geometry type> <input file name>)"),
     };
 
      match geometry_type {
@@ -97,8 +97,8 @@ fn main() {
             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);
-            physics_loop::<Mesh2D>(particle_input_array, material, options, output_units);
+            let (particle_input_array, material, options, output_units) = input::input::<geometry::ParryMesh2D>(input_file);
+            physics_loop::<ParryMesh2D>(particle_input_array, material, options, output_units);
         },
         GeometryType::SPHERE => {
             let (particle_input_array, material, options, output_units) = input::input::<Sphere>(input_file);
@@ -115,8 +115,8 @@ fn main() {
             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);
-            physics_loop::<HomogeneousMesh2D>(particle_input_array, material, options, output_units);
+            let (particle_input_array, material, options, output_units) = input::input::<geometry::ParryHomogeneousMesh2D>(input_file);
+            physics_loop::<ParryHomogeneousMesh2D>(particle_input_array, material, options, output_units);
         }
     }
 }
diff --git a/src/material.rs b/src/material.rs
index e6931fa..1367f56 100644
--- a/src/material.rs
+++ b/src/material.rs
@@ -300,7 +300,7 @@ impl <T: Geometry> Material<T> {
 
             let stopping_power = match electronic_stopping_mode {
                 //Biersack-Varelas Interpolation
-                ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low*ck),
+                ElectronicStoppingMode::INTERPOLATED => ck/(1./S_high + 1./S_low),
                 //Oen-Robinson
                 ElectronicStoppingMode::LOW_ENERGY_LOCAL => S_low*ck,
                 //Lindhard-Scharff
@@ -339,14 +339,28 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
     let leaving = !inside_now & inside_old;
     let entering = inside_now & !inside_old;
 
+    assert!(!x.is_nan(), "Particle x is NaN before surface binding energy calculation.");
+    assert!(!y.is_nan(), "Particle y is NaN before surface binding energy calculation.");
+    assert!(!z.is_nan(), "Particle z is NaN before surface binding energy calculation.");
+
+    assert!(!cosx.is_nan(), "Particle dirx is NaN before surface binding energy calculation.");
+    assert!(!cosy.is_nan(), "Particle diry is NaN before surface binding energy calculation.");
+    assert!(!cosz.is_nan(), "Particle dirz is NaN before surface binding energy calculation.");
+
     if entering {
         if particle_1.backreflected {
             particle_1.backreflected = false;
         } else {
-            let (x2, y2, z2) = material.closest_point(x, y, z);
-            let dx = x2 - x;
-            let dy = y2 - y;
-            let dz = z2 - z;
+            let (x2, y2, z2) = material.closest_point(x_old, y_old, z_old);
+            assert!(!x2.is_nan(), "Closest x is NaN in normal calculation.");
+            assert!(!y2.is_nan(), "Closest y is NaN in normal calculation.");
+            assert!(!z2.is_nan(), "Closest z is NaN in normal calculation.");
+            let dx = x2 - x_old;
+            let dy = y2 - y_old;
+            let dz = z2 - z_old;
+            assert!(!dx.is_nan(), "dx is NaN in normal calculation.");
+            assert!(!dy.is_nan(), "dy is NaN in normal calculation.");
+            assert!(!dz.is_nan(), "dz is NaN in normal calculation.");
             let mag = (dx*dx + dy*dy + dz*dz).sqrt();
             let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag;
 
@@ -360,10 +374,17 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
     if leaving {
 
         let (x2, y2, z2) = material.closest_point(x, y, z);
+        assert!(!x2.is_nan(), "Closest x is NaN in normal calculation.");
+        assert!(!y2.is_nan(), "Closest y is NaN in normal calculation.");
+        assert!(!z2.is_nan(), "Closest z is NaN in normal calculation.");
         let dx = x2 - x;
         let dy = y2 - y;
         let dz = z2 - z;
+        assert!(!dx.is_nan(), "dx is NaN in normal calculation.");
+        assert!(!dy.is_nan(), "dy is NaN in normal calculation.");
+        assert!(!dz.is_nan(), "dz is NaN in normal calculation.");
         let mag = (dx*dx + dy*dy + dz*dz).sqrt();
+
         let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag;
 
         let leaving_energy = match material.surface_binding_model {
@@ -414,4 +435,6 @@ pub fn boundary_condition_planar<T: Geometry>(particle_1: &mut particle::Particl
     if !particle_1.stopped & !particle_1.left {
         surface_binding_energy(particle_1, material);
     }
+
+    assert!(!particle_1.pos.x.is_nan(), "Particle position is NaN: ({}, {}, {}) old: ({}, {}, {}), left? {} stopped? {}", x, y, z, particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z, particle_1.left, particle_1.stopped);
 }
diff --git a/src/parry.rs b/src/parry.rs
index 8686683..2a33721 100644
--- a/src/parry.rs
+++ b/src/parry.rs
@@ -124,6 +124,15 @@ impl Geometry for ParryBall {
         let (x_, y_, z_) = (point_projection.point.x, point_projection.point.y, point_projection.point.z);
         (x_ as f64, y_ as f64, z_ as f64)
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let r = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt();
+        let R = self.radius;
+        let ux = x/r;
+        let uy = y/r;
+        let uz = z/r;
+        (ux, uy, uz)
+    }
 }
 
 
@@ -263,6 +272,22 @@ impl Geometry for ParryTriMesh {
         let (x_, y_, z_) = (point_projection.point.x, point_projection.point.y, point_projection.point.z);
         (x_ as f64, y_ as f64, z_ as f64)
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let p = Point::new(x, y, z);
+        let point_projection = self.trimesh.project_local_point(&p, false);
+
+        let dx = point_projection.point.x - x;
+        let dy = point_projection.point.y - y;
+        let dz = point_projection.point.z - z;
+        let mag = (dx*dx + dy*dy + dz*dz).sqrt();
+
+        if point_projection.is_inside {
+            (dx/mag, dy/mag, dz/mag)
+        } else {
+            (-dx/mag, -dy/mag, -dz/mag)
+        }
+    }
 }
 
 fn inside_trimesh(trimesh: &TriMesh, x: f64, y: f64, z: f64) -> bool {
diff --git a/src/particle.rs b/src/particle.rs
index 9de3daa..ba3bdeb 100644
--- a/src/particle.rs
+++ b/src/particle.rs
@@ -65,7 +65,7 @@ pub struct ParticleInput {
 }
 
 /// Particle object. Particles in rustbca include incident ions and material atoms.
-#[derive(Clone)]
+#[derive(Clone, Debug)]
 pub struct Particle {
     pub m: f64,
     pub Z: f64,
@@ -271,6 +271,9 @@ impl Particle {
         //In order to keep average denisty constant, must add back previous asymptotic deflection
         let distance_traveled = mfp + self.asymptotic_deflection - asymptotic_deflection;
 
+        assert!(!distance_traveled.is_nan(), "Travel distance is NaN, mfp: {}, asymp. defl: {}, old asymp. defl: {}", mfp, self.asymptotic_deflection, asymptotic_deflection);
+        assert!(!self.dir_old.x.is_nan(), "Particle direction is NaN, dir: ({}, {}, {})", self.dir_old.x, self.dir_old.y, self.dir_old.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;
@@ -294,10 +297,20 @@ pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) {
     let u1x = (E/(E + Es)).sqrt()*particle.dir.x + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.x;
     let u1y = (E/(E + Es)).sqrt()*particle.dir.y + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.y;
     let u1z = (E/(E + Es)).sqrt()*particle.dir.z + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.z;
+
+    particle.E += Es;
+
+    if u1x.is_nan() {
+        println!("normal: ({}, {}, {})", normal.x, normal.y, normal.z);
+        dbg!(&particle);
+
+    }
+
     particle.dir.x = u1x;
     particle.dir.y = u1y;
     particle.dir.z = u1z;
-    particle.E += Es;
+
+    assert!(!u1x.is_nan(), "Particle direction after surface refraction is NaN. dir: ({}, {}, {}), old dir: ({}, {}, {}), normal: ({}, {}, {}), costheta: {}", u1x, u1y, u1z, particle.dir_old.x, particle.dir_old.y, particle.dir_old.z, normal.x, normal.y, normal.z, costheta);
 }
 
 /// Calcualte the refraction angle based on the surface binding energy of the material.
diff --git a/src/sphere.rs b/src/sphere.rs
index 92553b6..6f4055a 100644
--- a/src/sphere.rs
+++ b/src/sphere.rs
@@ -120,4 +120,13 @@ impl Geometry for Sphere {
         let uz = z/r;
         (ux*R, uy*R, uz*R)
     }
+
+    fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+        let r = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt();
+        let R = self.radius;
+        let ux = x/r;
+        let uy = y/r;
+        let uz = z/r;
+        (ux, uy, uz)
+    }
 }
diff --git a/src/structs.rs b/src/structs.rs
index 47f20a0..fa5993d 100644
--- a/src/structs.rs
+++ b/src/structs.rs
@@ -1,5 +1,5 @@
 /// 3D vector.
-#[derive(Clone, Copy)]
+#[derive(Clone, Copy, Debug)]
 pub struct Vector {
     pub x: f64,
     pub y: f64,
@@ -45,7 +45,7 @@ impl Vector {
 }
 
 /// TrajectoryElement is a trajectory-tracking object that includes x, y, z, and the current energy.
-#[derive(Clone)]
+#[derive(Clone, Debug)]
 pub struct TrajectoryElement {
     pub E: f64,
     pub x: f64,
@@ -65,7 +65,7 @@ impl TrajectoryElement {
 }
 
 /// Energy loss is an output tracker that tracks the separate nuclear and electronic energy losses.
-#[derive(Clone)]
+#[derive(Clone, Debug)]
 pub struct EnergyLoss {
     pub En: f64,
     pub Ee: f64,
diff --git a/src/tests.rs b/src/tests.rs
index fd263f0..5bd3200 100644
--- a/src/tests.rs
+++ b/src/tests.rs
@@ -3,6 +3,115 @@ use super::*;
 #[cfg(test)]
 use float_cmp::*;
 
+use parry2d_f64::shape::Polyline;
+use parry2d_f64::math::Isometry;
+use parry2d_f64::math::Point as Point2d;
+
+/*
+        let number_boundary_points = boundary_points_converted.clone().len() as u32;
+        let boundary = Polygon::new(LineString::from(boundary_points_converted), vec![]);
+        let mut linked_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::<Vec<[u32; 2]>>();
+        let boundary2 = Polyline::new(boundary_points_converted2, Some(linked_points));
+*/
+
+#[test]
+fn test_parry2d() {
+    let points: Vec<Point2d<f64>> = vec![Point2d::new(-1.0, -1.0), Point2d::new(-1.0, 1.0), Point2d::new(1.0, 1.0), Point2d::new(1.0, -1.0)];
+    //let mut indices: Vec<[u32; 2]> = vec![];
+    //let mut indices: Vec<[u32; 2]> = vec![[0, 1], [1, 2], [2, 3], [3, 0]];
+    let indices: Vec<[u32; 2]> = vec![[0, 3], [3, 2], [2, 1], [1, 0]];
+    //indices.reverse();
+
+    //(0.0000006261797114005236, 0.000006009591179670447)
+    let query_point = Point2d::new(0.0, 0.0);
+
+    let polyline = Polyline::new(points, Some(indices));
+
+    let (point_projection, (_, _)) = polyline.project_local_point_assuming_solid_interior_ccw(
+        query_point
+    );
+
+    assert!(point_projection.is_inside);
+
+    assert!(polyline.aabb(&Isometry::identity()).contains_local_point(&query_point));
+    //assert!(point_projection.is_inside);
+
+    let test_geometry = vec![(0.25, 0.2), (0.75, 0.2), (0.8, 0.8), (0.5, 0.5), (0.2, 0.8)];
+    let num_segments = 5;
+
+    let geometry_input_homogeneous_2D = geometry::HomogeneousMesh2DInput {
+        length_unit: "M".to_string(),
+        points: test_geometry.clone(),
+        densities: vec![0.03, 0.03],
+        simulation_boundary_points: vec![(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)],
+        electronic_stopping_correction_factor: 1.0,
+    };
+
+    let material_parameters = material::MaterialParameters{
+        energy_unit: "EV".to_string(),
+        mass_unit: "AMU".to_string(),
+        Eb: vec![0.0, 0.0],
+        Es: vec![2.0, 4.0],
+        Ec: vec![1.0, 1.0],
+        Ed: vec![0.0, 0.0],
+        Z: vec![29., 1.],
+        m: vec![63.54, 1.0008],
+        interaction_index: vec![0, 0],
+        surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET},
+        bulk_binding_model: BulkBindingModel::INDIVIDUAL,
+    };
+
+    let material_homogeneous_2D: material::Material<geometry::ParryHomogeneousMesh2D> = material::Material::<geometry::ParryHomogeneousMesh2D>::new(&material_parameters, &geometry_input_homogeneous_2D);
+
+    let query_points_outside = vec![[0.5, 0.1, 0.0], [0.9, 0.5, 0.0], [0.6, 0.7, 0.0], [0.4, 0.7, 0.0], [0.1, 0.5, 0.0]];
+    let query_points_inside = vec![[0.5, 0.21, 0.0], [0.75, 0.21, 0.0], [0.51, 0.5, 0.0], [0.49, 0.5, 0.0], [0.25, 0.21, 0.0]];
+
+    for query_point in query_points_outside.clone() {
+        assert!(!material_homogeneous_2D.inside(query_point[0], query_point[1], query_point[2]), "Point ({}, {}, {}) failed outsidedness check.", query_point[0], query_point[1], query_point[2]);
+    }
+
+    for query_point in query_points_inside.clone() {
+        assert!(material_homogeneous_2D.inside(query_point[0], query_point[1], query_point[2]), "Point ({}, {}, {}) failed outsidedness check.", query_point[0], query_point[1], query_point[2]);
+    }
+
+
+    let normals_outside = query_points_outside
+        .iter()
+        .map( |query_point| {
+            material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2])
+        }
+    ).collect::<Vec<(f64, f64, f64)>>();
+
+    let normals_inside = query_points_inside
+    .iter()
+    .map( |query_point| {
+        material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2])
+    }
+    ).collect::<Vec<(f64, f64, f64)>>();
+
+    let test_normals = (0..num_segments)
+        .zip(1..num_segments + 1)
+        .map(|(i, j)| {
+            let dx = test_geometry[i % num_segments].0 - test_geometry[j % num_segments].0;
+            let dy = test_geometry[i % num_segments].1 - test_geometry[j % num_segments].1;
+            let mag = (dx*dx + dy*dy).sqrt();
+            (-dy/mag, dx/mag, 0.0)
+        }).collect::<Vec<(f64, f64, f64)>>();
+
+    for (test_normal, normal) in test_normals.iter().zip(normals_outside) {
+        assert!(approx_eq!(f64, normal.0, test_normal.0, epsilon=1E-9), "Normal vector check failed. Calculated Normal: ({}, {}, {}) Test: ({}, {}, {})", normal.0, normal.1, normal.2, test_normal.0, test_normal.1, test_normal.2);
+        assert!(approx_eq!(f64, normal.1, test_normal.1, epsilon=1E-9));
+        assert!(approx_eq!(f64, normal.2, test_normal.2, epsilon=1E-9));
+    }
+
+    for (test_normal, normal) in test_normals.iter().zip(normals_inside) {
+        assert!(approx_eq!(f64, normal.0, test_normal.0, epsilon=1E-9), "Normal vector check failed. Calculated Normal: ({}, {}, {}) Test: ({}, {}, {})", normal.0, normal.1, normal.2, test_normal.0, test_normal.1, test_normal.2);
+        assert!(approx_eq!(f64, normal.1, test_normal.1, epsilon=1E-9));
+        assert!(approx_eq!(f64, normal.2, test_normal.2, epsilon=1E-9));
+    }
+
+}
+
 
 #[test]
 #[cfg(feature = "cpr_rootfinder")]
@@ -439,6 +548,61 @@ fn test_spherical_geometry() {
 
 }
 
+#[test]
+fn extended_test_2D_geometry() {
+    let mass = 1.;
+    let Z = 1.;
+    let E = 10.*EV;
+    let Ec = 1.*EV;
+    let Es = 5.76*EV;
+    let x = 1.5e-6;
+    let y = 10.0e-6;
+    let z = 0.;
+    let cosx = 0.0;
+    let cosy = -1.0;
+    let cosz = 0.;
+
+    let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0);
+
+    let material_parameters = material::MaterialParameters{
+        energy_unit: "EV".to_string(),
+        mass_unit: "AMU".to_string(),
+        Eb: vec![0.0, 0.0],
+        Es: vec![2.0, 4.0],
+        Ec: vec![1.0, 1.0],
+        Ed: vec![0.0, 0.0],
+        Z: vec![29., 1.],
+        m: vec![63.54, 1.0008],
+        interaction_index: vec![0, 0],
+        surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET},
+        bulk_binding_model: BulkBindingModel::INDIVIDUAL,
+    };
+
+    let thickness: f64 = 1000.;
+    let depth: f64 = 1000.;
+
+    let geometry_input_2D = geometry::Mesh2DInput {
+        length_unit: "MICRON".to_string(),
+        points: vec![(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (4.0, 1.0), (3.0, 1.0), (3.0, 10.0), (2.0, 10.0), (2.0, 1.0), (1.0, 1.0), (1.0, 10.0), (0.0, 10.0), (0.0, 1.0)],
+        triangles: vec![(13, 12, 11), (13, 11, 10), (0, 13, 10), (0, 10, 1), (1, 10, 9), (1, 9, 2), (9, 8, 7), (9, 7, 6), (2, 9, 6), (2, 6, 3), (3, 6, 5), (3, 5, 4)],
+        densities: vec![vec![3e10, 3e10]; 12],
+        boundary: vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13],
+        simulation_boundary_points: vec![(-0.1, -0.1), (4.1, -0.1), (4.1, 10.1), (-0.1, 10.1)],
+        energy_barrier_thickness: 4e-4,
+        electronic_stopping_correction_factors: vec![1.0; 12],
+    };
+
+    let material_2D: material::Material<geometry::ParryMesh2D> = material::Material::<geometry::ParryMesh2D>::new(&material_parameters, &geometry_input_2D);
+
+
+    assert!(material_2D.inside(1.5e-6, 0.999e-6, 0.0));
+    assert!(!material_2D.inside(1.5e-6, 1.0001e-6, 0.0));
+
+    material::surface_binding_energy(&mut particle_1, &material_2D);
+
+    material::boundary_condition_planar(&mut particle_1, &material_2D);
+}
+
 #[test]
 fn test_geometry() {
     let mass = 1.;
@@ -477,7 +641,7 @@ fn test_geometry() {
         triangles: vec![(0, 1, 2), (0, 2, 3)],
         densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]],
         boundary: vec![0, 1, 2, 3],
-        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)],
+        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)],
         energy_barrier_thickness: 10.,
         electronic_stopping_correction_factors: vec![1.0, 1.0],
     };
@@ -486,7 +650,7 @@ fn test_geometry() {
         length_unit: "ANGSTROM".to_string(),
         points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)],
         densities: vec![0.03, 0.03],
-        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)],
+        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)],
         electronic_stopping_correction_factor: 1.0,
     };
 
@@ -504,8 +668,8 @@ fn test_geometry() {
         electronic_stopping_correction_factor: 1.0
     };
 
-    let material_2D: material::Material<geometry::Mesh2D> = material::Material::<geometry::Mesh2D>::new(&material_parameters, &geometry_input_2D);
-    let material_homogeneous_2D: material::Material<geometry::HomogeneousMesh2D> = material::Material::<geometry::HomogeneousMesh2D>::new(&material_parameters, &geometry_input_homogeneous_2D);
+    let material_2D: material::Material<geometry::ParryMesh2D> = material::Material::<geometry::ParryMesh2D>::new(&material_parameters, &geometry_input_2D);
+    let material_homogeneous_2D: material::Material<geometry::ParryHomogeneousMesh2D> = material::Material::<geometry::ParryHomogeneousMesh2D>::new(&material_parameters, &geometry_input_homogeneous_2D);
     let material_1D: material::Material<geometry::Mesh1D> = material::Material::<geometry::Mesh1D>::new(&material_parameters, &geometry_input_1D);
     let mut material_0D: material::Material<geometry::Mesh0D> = material::Material::<geometry::Mesh0D>::new(&material_parameters, &geometry_input_0D);
     material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM;
@@ -529,7 +693,12 @@ fn test_geometry() {
     assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_0D.geometry.get_densities(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_0D.geometry.get_total_density(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_0D.geometry.get_concentrations(0., 0., 0.));
-    assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_0D.geometry.closest_point(-10., 0., 5.));
+    let (x2d, y2d, z2d) = material_2D.geometry.closest_point(0., 0., 0.);
+    let (x0d, y0d, z0d) = material_0D.geometry.closest_point(0., 0., 0.);
+    assert!(approx_eq!(f64, x2d, x0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d);
+    assert!(approx_eq!(f64, y2d, y0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d);
+    assert!(approx_eq!(f64, z2d, z0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d);
+
     assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_0D.geometry.get_densities(-10., 0., 5.));
     assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_0D.geometry.get_ck(-10., 0., 5.));
 
@@ -537,7 +706,14 @@ fn test_geometry() {
     assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_homogeneous_2D.geometry.get_densities(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_homogeneous_2D.geometry.get_total_density(0., 0., 0.));
     assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_homogeneous_2D.geometry.get_concentrations(0., 0., 0.));
-    assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_homogeneous_2D.geometry.closest_point(-10., 0., 5.));
+
+    let (x2d, y2d, z2d) = material_2D.geometry.closest_point(-10., 0., 5.);
+    let (x0d, y0d, z0d) = material_0D.geometry.closest_point(-10., 0., 5.);
+    assert!(approx_eq!(f64, x2d, x0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d);
+    assert!(approx_eq!(f64, y2d, y0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d);
+    assert!(approx_eq!(f64, z2d, z0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d);
+
+
     assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_homogeneous_2D.geometry.get_densities(-10., 0., 5.));
     assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_homogeneous_2D.geometry.get_ck(-10., 0., 5.));
 }
@@ -580,7 +756,7 @@ fn test_surface_binding_energy_barrier() {
         triangles: vec![(0, 1, 2), (0, 2, 3)],
         densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]],
         boundary: vec![0, 1, 2, 3],
-        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)],
+        simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)],
         energy_barrier_thickness: 10.,
         electronic_stopping_correction_factors: vec![1.0, 1.0],
     };
@@ -591,7 +767,7 @@ fn test_surface_binding_energy_barrier() {
         electronic_stopping_correction_factor: 1.0
     };
 
-    let material_2D: material::Material<geometry::Mesh2D> = material::Material::<geometry::Mesh2D>::new(&material_parameters, &geometry_input_2D);
+    let material_2D: material::Material<geometry::ParryMesh2D> = material::Material::<geometry::ParryMesh2D>::new(&material_parameters, &geometry_input_2D);
     let mut material_0D: material::Material<geometry::Mesh0D> = material::Material::<geometry::Mesh0D>::new(&material_parameters, &geometry_input_0D);
     material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM;
 
@@ -662,6 +838,8 @@ fn test_surface_binding_energy_barrier() {
     assert!(material_0D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.));
 }
 
+//Deprecated with geo
+/*
 #[test]
 fn test_triangle_contains() {
     let triangle_1 = geometry::Triangle2D::new((0., 2., 0., 2., 0., 0.));
@@ -685,6 +863,7 @@ fn test_triangle_distance_to() {
     assert!(approx_eq!(f64, triangle_1.distance_to(2., 0.), 0., epsilon=1E-12), "{}", triangle_1.distance_to(2., 0.));
     assert!(approx_eq!(f64, triangle_1.distance_to(0., 2.), 0., epsilon=1E-12), "{}", triangle_1.distance_to(0., 2.));
 }
+*/
 
 #[test]
 fn test_surface_refraction() {
@@ -805,20 +984,20 @@ fn test_momentum_conservation() {
             triangles: vec![(0, 1, 2), (0, 2, 3)],
             points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)],
             densities: vec![vec![0.06306], vec![0.06306]],
-            boundary: vec![0, 1, 2, 3],
-            simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)],
+            boundary: vec![3, 2, 1, 0],
+            simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)],
             electronic_stopping_correction_factors: vec![0.0, 0.0],
             energy_barrier_thickness: 0.,
         };
 
-        let material_1: material::Material<geometry::Mesh2D> = material::Material::<geometry::Mesh2D>::new(&material_parameters, &geometry_input);
+        let material_1: material::Material<geometry::ParryMesh2D> = material::Material::<geometry::ParryMesh2D>::new(&material_parameters, &geometry_input);
 
         for high_energy_free_flight_paths in vec![true, false] {
             for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] {
                 for scattering_integral in vec![ScatteringIntegral::MENDENHALL_WELLER, ScatteringIntegral::GAUSS_MEHLER{n_points: 10}, ScatteringIntegral::GAUSS_LEGENDRE] {
                     for root_finder in vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}] {
 
-                        println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral);
+                        //println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral);
 
                         let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0);
 
@@ -883,9 +1062,11 @@ fn test_momentum_conservation() {
 
                         let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options);
 
+                        /*
                         println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal,
                             binary_collision_geometries[0].impact_parameter/ANGSTROM,
                             binary_collision_geometries[0].mfp/ANGSTROM);
+                        */
 
                         let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1,
                             &binary_collision_geometries[0], &options);
@@ -898,12 +1079,13 @@ fn test_momentum_conservation() {
                         let binary_collision_result = bca::calculate_binary_collision(&particle_1,
                             &particle_2, &binary_collision_geometries[0], &options).unwrap();
 
+                        /*
                         println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV,
                             binary_collision_result.psi,
                             binary_collision_result.psi_recoil);
 
                         println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV);
-
+                        */
                         //Energy transfer to recoil
                         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);
 
@@ -924,11 +1106,13 @@ fn test_momentum_conservation() {
 
                         let final_momentum = mom1_1.add(&mom2_1);
 
+                        /*
                         println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV);
                         println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude());
                         println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude());
                         println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude());
                         println!();
+                        */
 
                         assert!(approx_eq!(f64, initial_momentum.x, final_momentum.x, epsilon = 1E-12));
                         assert!(approx_eq!(f64, initial_momentum.y, final_momentum.y, epsilon = 1E-12));