spdcalc 2.0.1

SPDCalc, the Spontaneous Parametric Downconversion Calculator
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

use crate::{Complex, SPDCError};
use na::DMatrix;

/// Calculate the schmidt number of a flat vector joint spectral amplitude values (square matrix)
pub fn schmidt_number<T: AsRef<[Complex<f64>]>>(amplitudes: T) -> Result<f64, SPDCError> {
  use num::integer::Roots;
  let len = amplitudes.as_ref().len();
  let dim = len.sqrt();
  if len != dim * dim {
    return Err(SPDCError("Spectrum provided is not square".into()));
  }
  let jsa_mag: Vec<f64> = amplitudes.as_ref().iter().map(|j| j.norm()).collect();
  let svd = DMatrix::from_row_slice(dim, dim, &jsa_mag)
    .try_svd(false, false, f64::EPSILON, 10_000)
    .ok_or(SPDCError(
      "SVD did not converge while calculating schmidt number".into(),
    ))?;
  let norm_sq = svd.singular_values.norm_squared();
  let kinv = svd.singular_values.fold(0., |acc, x| acc + x.powi(4));
  Ok(norm_sq * norm_sq / kinv)
}

#[cfg(test)]
mod tests {
  use super::*;
  use crate::*;
  use dim::f64prefixes::NANO;
  use dim::ucum::M;
  extern crate float_cmp;
  use float_cmp::*;
  use math::Integrator;

  #[test]
  fn shmidt_number_test() {
    let mut spdc = SPDC::default();
    spdc.crystal_setup.crystal = CrystalType::KTP;
    spdc.crystal_setup.phi = Angle::new(0.);
    spdc.crystal_setup.theta = Angle::new(PI / 2.);
    spdc.assign_optimum_idler().unwrap();
    spdc.assign_optimum_periodic_poling().unwrap();

    let wavelength_range = WavelengthSpace::new(
      (1490.86 * NANO * M, 1609.14 * NANO * M, 100),
      (1495.05 * NANO * M, 1614.03 * NANO * M, 100),
    );
    // let wavelength_range = calc_plot_config_for_jsi(&spdc_setup, 100, 0.5);
    // dbg!(wavelength_range);
    // dbg!(spdc_setup);

    let spectrum = spdc.joint_spectrum(Integrator::default());
    let amplitudes = spectrum.jsa_range(wavelength_range);
    let sn = schmidt_number(amplitudes).expect("Could not calc schmidt number");

    let actual = sn;
    let expected = 1.149;

    assert!(
      approx_eq!(f64, actual, expected, ulps = 2, epsilon = 1e-3),
      "actual: {}, expected: {}",
      actual,
      expected
    );
  }
}