Extend spectrum struct using ndarray.

Cargo.lock

@@ -30,6 +30,59 @@ dependencies = [
  "hashbrown",
 ]
 
+[[package]]
+name = "matrixmultiply"
+version = "0.3.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "090126dc04f95dc0d1c1c91f61bdd474b3930ca064c1edc8a849da2c6cbe1e77"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
+[[package]]
+name = "ndarray"
+version = "0.15.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32"
+dependencies = [
+ "matrixmultiply",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "rawpointer",
+]
+
+[[package]]
+name = "ndarray-interp"
+version = "0.1.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "16eb7477aebbc0c82bcc496d2d727e76875f800060a7eeac125d1f884a129750"
+dependencies = [
+ "ndarray",
+ "num-traits",
+ "thiserror",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.4.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "02e0d21255c828d6f128a1e41534206671e8c3ea0c62f32291e808dc82cff17d"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.45"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9"
+dependencies = [
+ "autocfg",
+ "num-traits",
+]
+
 [[package]]
 name = "num-traits"
 version = "0.2.15"
@@ -43,6 +96,8 @@ dependencies = [
 name = "opossum"
 version = "0.1.0"
 dependencies = [
+ "ndarray",
+ "ndarray-interp",
  "petgraph",
  "uom",
 ]
@@ -77,6 +132,12 @@ dependencies = [
  "proc-macro2",
 ]
 
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
 [[package]]
 name = "serde"
 version = "1.0.171"
@@ -105,6 +166,26 @@ dependencies = [
  "unicode-ident",
 ]
 
+[[package]]
+name = "thiserror"
+version = "1.0.43"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a35fc5b8971143ca348fa6df4f024d4d55264f3468c71ad1c2f365b0a4d58c42"
+dependencies = [
+ "thiserror-impl",
+]
+
+[[package]]
+name = "thiserror-impl"
+version = "1.0.43"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "463fe12d7993d3b327787537ce8dd4dfa058de32fc2b195ef3cde03dc4771e8f"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
 [[package]]
 name = "typenum"
 version = "1.16.0"

Cargo.toml

@@ -12,4 +12,8 @@ publish = false
 
 [dependencies]
 petgraph = {version="0.6.3", features = ["serde-1"]} # the graph library
-uom = "0"
\ No newline at end of file
+uom = "0"
+
+# for spectrum
+ndarray="0"
+ndarray-interp="0"
\ No newline at end of file

examples/spectrum_test.rs

@@ -9,7 +9,7 @@ fn main() -> Result<(), OpossumError> {
         Length::new::<nanometer>(400.0)..Length::new::<nanometer>(410.0),
         Length::new::<nanometer>(1.0),
     );
-    s.set_single_peak(Length::new::<nanometer>(405.0), Energy::new::<joule>(1.0));
+    s.set_single_peak(Length::new::<nanometer>(399.99999), Energy::new::<joule>(1.0))?;
     println!("{}", s);
     Ok(())
 }

src/error.rs

@@ -6,6 +6,7 @@ pub enum OpossumError {
     OpticGroup(String),
     OpticPort(String),
     Analysis(String),
+    Spectrum(String),
     Other(String),
 }
 
@@ -24,6 +25,9 @@ impl Display for OpossumError {
             OpossumError::Analysis(m) => {
                 write!(f, "Opossum Error::Analysis::{}", m)
             }
+            OpossumError::Spectrum(m) => {
+                write!(f, "Opossum Spectrum::{}", m)
+            }
             OpossumError::Other(m) => write!(f, "Opossum Error::Other::{}", m),
         }
     }
@@ -34,4 +38,4 @@ impl std::convert::From<String> for OpossumError {
     fn from(msg: String) -> Self {
         Self::Other(msg)
     }
-}
\ No newline at end of file
+}

src/spectrum.rs

+use ndarray::Array1;
 use std::fmt::Display;
 use std::ops::Range;
 use uom::fmt::DisplayStyle::Abbreviation;
 use uom::num_traits::Zero;
 use uom::si::energy::joule;
+use uom::si::length::meter;
 use uom::si::{
     f64::{Energy, Length},
     length::nanometer,
 };
+
+use crate::error::OpossumError;
+type Result<T> = std::result::Result<T, OpossumError>;
+
 pub struct Spectrum {
-    start: Length,
-    dlambda: Length,
-    data: Vec<Energy>,
+    data: Array1<Energy>,
+    lambdas: Array1<f64>,
 }
 
 impl Spectrum {
     pub fn new(range: Range<Length>, resolution: Length) -> Self {
+        let l = Array1::range(
+            range.start.get::<meter>(),
+            range.end.get::<meter>(),
+            resolution.get::<meter>(),
+        );
+        let length = l.len();
         Self {
-            start: range.start,
-            dlambda: resolution,
-            data: vec![Energy::zero(); ((range.end-range.start)/resolution).value as usize]
+            lambdas: l,
+            data: Array1::zeros(length),
         }
     }
-    pub fn set_single_peak(&mut self, wavelength: Length, energy: Energy) {
-        let index=((wavelength-self.start) / self.dlambda).value as usize;
-        self.data[index]=energy; 
+    pub fn set_single_peak(&mut self, wavelength: Length, energy: Energy) -> Result<()> {
+        let spectrum_range = self.lambdas.first().unwrap()..self.lambdas.last().unwrap();
+        if !spectrum_range.contains(&&wavelength.get::<meter>()) {
+            return Err(OpossumError::Other(
+                "wavelength is not in spectrum range".into(),
+            ));
+        }
+        if energy.is_sign_negative() {
+            return Err(OpossumError::Other("energy must be positive".into()));
+        }
+        let wavelength_in_meters = wavelength.get::<meter>();
+        let idx = self
+            .lambdas
+            .clone()
+            .into_iter()
+            .position(|w| w >= wavelength_in_meters);
+        if let Some(idx) = idx {
+            if idx == 0 {
+                self.data[0] = energy;
+            } else {
+                let lower_lambda = self.lambdas.get(idx - 1).unwrap();
+                let upper_lambda = self.lambdas.get(idx).unwrap();
+                let delta = upper_lambda - lower_lambda;
+                self.data[idx - 1] = energy * (1.0 - (wavelength_in_meters - lower_lambda) / delta);
+                self.data[idx] = energy * (wavelength_in_meters - lower_lambda) / delta;
+            }
+            Ok(())
+        } else {
+            Err(OpossumError::Other("insertion point not found".into()))
+        }
     }
     pub fn total_energy(&self) -> Energy {
-        let mut total_energy=Energy::zero();
+        let mut total_energy = Energy::zero();
         for data in self.data.iter() {
-            total_energy+=*data;
+            total_energy += *data;
         }
         total_energy
     }
@@ -39,14 +76,18 @@ impl Display for Spectrum {
         let fmt_length = Length::format_args(nanometer, Abbreviation);
         let fmt_energy = Energy::format_args(joule, Abbreviation);
         for value in self.data.iter().enumerate() {
-            let wavelength=self.start + value.0 as f64 * self.dlambda;
             write!(
                 f,
                 "{:7.2} -> {}\n",
-                fmt_length.with(wavelength),
+                fmt_length.with(Length::new::<meter>(self.lambdas[value.0])),
                 fmt_energy.with(*value.1)
-            ).unwrap();
+            )
+            .unwrap();
         }
-        write!(f, "\nTotal energy: {}",fmt_energy.with(self.total_energy()))
+        write!(
+            f,
+            "\nTotal energy: {}",
+            fmt_energy.with(self.total_energy())
+        )
     }
 }