Add check on maximum number density to catch errors.

Author: drobnyjt

src/consts.rs

@@ -44,4 +44,6 @@ pub const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2;
 /// Gauss-Legendre Points
 pub const GL_X: [f64; 5] = [0./2. + 1./2., -0.538469/2. + 1./2., 0.538469/2. + 1./2., -0.90618/2. + 1./2., 0.90618/2. + 1./2.];
 /// Gauss-Legendre Weights
-pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.];
+pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.];
+/// Enforced maximum number density
+pub const MAX_DENSITY: f64 = 1e32;

src/geometry.rs

@@ -65,8 +65,9 @@ impl Geometry for Mesh0D {
 
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
         assert!(!densities.is_empty(), "Input Error: density list empty.");
-
+        
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
 
         let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;
 
@@ -154,6 +155,7 @@ impl Geometry for Mesh1D {
         let densities: Vec<Vec<f64>> = geometry_input.densities
             .iter()
             .map( |row| row.iter().map(|element| element/(length_unit).powi(3)).collect() ).collect();
+    
 
         //Assert all layer density lists are equal length
         assert!(
@@ -170,6 +172,7 @@ impl Geometry for Mesh1D {
             layer_bottom += layer_thickness*length_unit;
 
             let total_density: f64 = densities.iter().sum();
+            assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
             let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             layers.push(Layer1D::new(layer_top, layer_bottom, densities, concentrations, ck));
@@ -308,6 +311,7 @@ impl Geometry for HomogeneousMesh2D {
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
 
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
 
         let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;
 
@@ -476,6 +480,7 @@ impl Geometry for Mesh2D {
                 y3*length_unit,
             );
             let total_density: f64 = densities.iter().sum();
+            assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
             let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck));

src/parry.rs

@@ -77,6 +77,7 @@ impl Geometry for ParryBall {
         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();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
         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 radius = input.radius*length_unit;
@@ -201,6 +202,7 @@ impl Geometry for ParryTriMesh {
         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();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
         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 points = input.vertices.iter().map(|p| Point::new(p[0]*length_unit , p[1]*length_unit , p[2]*length_unit)).collect();

src/sphere.rs

@@ -63,13 +63,14 @@ impl Geometry for Sphere {
             "NM" => NM,
             "M" => 1.,
             _ => input.length_unit.parse()
-                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
+                .unwrap_or_else(|_| panic!("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())),
         };
 
         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();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
         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 radius = input.radius*length_unit;