spdcalc/
utils.rs

1//! General utilities
2use crate::math::lerp;
3use crate::{Frequency, RIndex, Speed, Wavelength, Wavenumber, TWO_PI};
4use dim::ucum::{self, C_, ONE, RAD};
5use dim::{ucum::UCUM, Dimensioned};
6use na::{Unit, Vector3};
7use rayon::iter::IntoParallelIterator;
8
9/// Extension to facilitate getting the x, y, and z components of a dimensioned vector
10pub trait DimVector {
11  /// The type of the units
12  type Units;
13  /// Get the x component as a dimensioned value
14  fn x(&self) -> Self::Units;
15  /// Get the y component as a dimensioned value
16  fn y(&self) -> Self::Units;
17  /// Get the z component as a dimensioned value
18  fn z(&self) -> Self::Units;
19}
20
21impl<U> DimVector for UCUM<Vector3<f64>, U> {
22  type Units = UCUM<f64, U>;
23  fn x(&self) -> Self::Units {
24    UCUM::new(self.value_unsafe().x)
25  }
26  fn y(&self) -> Self::Units {
27    UCUM::new(self.value_unsafe().y)
28  }
29  fn z(&self) -> Self::Units {
30    UCUM::new(self.value_unsafe().z)
31  }
32}
33
34/// Get the dimensioned vector from a scalar value and a direction.
35pub fn dim_vector<U>(val: UCUM<f64, U>, dir: Unit<Vector3<f64>>) -> UCUM<Vector3<f64>, U> {
36  UCUM::<Vector3<f64>, U>::new(*val.value_unsafe() * *dir)
37}
38
39/// convert from celsius to kelvin
40pub fn from_celsius_to_kelvin(c: f64) -> ucum::Kelvin<f64> {
41  ucum::Kelvin::new(c + 273.15)
42}
43
44/// convert from kelvin to celsius
45pub fn from_kelvin_to_celsius(k: ucum::Kelvin<f64>) -> f64 {
46  *(k / ucum::K) - 273.15
47}
48
49/// Get the wavenumber of light at frequency in a medium with refractive index
50pub fn frequency_to_wavenumber(omega: Frequency, n: RIndex) -> Wavenumber {
51  n * omega / C_
52}
53
54/// Get the frequency of light with wavenumber in a medium with refractive index
55pub fn wavenumber_to_frequency(k: Wavenumber, n: RIndex) -> Frequency {
56  k * C_ / n
57}
58
59/// Get the frequency of light at wavelength in a medium with refractive index
60pub fn wavelength_to_frequency(lambda: Wavelength, n: RIndex) -> Frequency {
61  TWO_PI * RAD * C_ / (lambda * n)
62}
63
64/// Get the wavelength of light at frequency in a medium with refractive index
65pub fn frequency_to_wavelength(omega: Frequency, n: RIndex) -> Wavelength {
66  TWO_PI * RAD * C_ / (omega * n)
67}
68
69/// Get the phase velocity of light from frequency and wavenumber
70pub fn phase_velocity(omega: Frequency, k: Wavenumber) -> Speed {
71  omega / k
72}
73
74/// Get the frequency of light from its vacuum wavelength
75pub fn vacuum_wavelength_to_frequency(lambda: Wavelength) -> Frequency {
76  wavelength_to_frequency(lambda, ONE)
77}
78
79/// Get the vacuum wavelength of light from its frequency
80pub fn frequency_to_vacuum_wavelength(omega: Frequency) -> Wavelength {
81  frequency_to_wavelength(omega, ONE)
82}
83
84/// Utility for creating evenly spaced steps between two endpoints over floating point numbers
85///
86/// ## Example:
87/// ```
88/// use spdcalc::utils::Steps;
89///
90/// let arr : Vec<f64> = Steps(0., 1., 100).into_iter().collect();
91/// assert_eq!(arr.len(), 100);
92/// assert!((arr[0] - 0.).abs() < 1e-12);
93/// assert!((arr[99] - 1.).abs() < 1e-12);
94/// ```
95#[derive(Debug, Copy, Clone, Hash, Serialize, Deserialize)]
96pub struct Steps<T>(pub T, pub T, pub usize);
97
98impl<T> From<(T, T, usize)> for Steps<T> {
99  fn from(args: (T, T, usize)) -> Self {
100    Self(args.0, args.1, args.2)
101  }
102}
103
104impl<T> Steps<T>
105where
106  T: std::ops::Mul<f64, Output = T>
107    + std::ops::Add<T, Output = T>
108    + std::ops::Div<f64, Output = T>
109    + std::ops::Sub<T, Output = T>
110    + Copy,
111{
112  /// Starting value
113  #[inline(always)]
114  pub fn start(&self) -> T {
115    self.0
116  }
117  /// Ending value
118  #[inline(always)]
119  pub fn end(&self) -> T {
120    self.1
121  }
122  /// Number of steps
123  #[inline(always)]
124  pub fn steps(&self) -> usize {
125    self.2
126  }
127  /// Number of steps
128  #[inline(always)]
129  pub fn len(&self) -> usize {
130    self.steps()
131  }
132  /// Is this empty?
133  #[inline(always)]
134  pub fn is_empty(&self) -> bool {
135    self.steps() == 0
136  }
137  /// Get the range as a tuple
138  #[inline(always)]
139  pub fn range(&self) -> (T, T) {
140    (self.start(), self.end())
141  }
142  /// Get the number of divisions (steps - 1)
143  #[inline(always)]
144  pub fn divisions(&self) -> usize {
145    self.steps() - 1
146  }
147
148  /// Get the value at a given index
149  pub fn value(&self, index: usize) -> T {
150    let (start, end) = self.range();
151    // if we have only one step... then just set progress to be zero
152    if self.steps() > 1 {
153      let index = index as f64;
154      let d = self.divisions() as f64;
155      (start * (d - index) + end * index) / d
156    } else {
157      start
158    }
159  }
160
161  /// Get the width of the gap between each step.
162  ///
163  /// ## Example:
164  /// ```
165  /// use spdcalc::utils::Steps;
166  ///
167  /// let steps = Steps(0., 4., 3); // 3 steps: |0| --- |2| --- |4|
168  /// let dx = 2.;
169  /// assert!((steps.division_width() - dx).abs() < 1e-12);
170  /// ```
171  pub fn division_width(&self) -> T {
172    (self.end() - self.start()) / (self.divisions() as f64)
173  }
174}
175
176impl<T> IntoIterator for Steps<T>
177where
178  T: std::ops::Mul<f64, Output = T>
179    + std::ops::Add<T, Output = T>
180    + std::ops::Div<f64, Output = T>
181    + std::ops::Sub<T, Output = T>
182    + Copy,
183{
184  type Item = T;
185  type IntoIter = Iterator1D<T>;
186
187  fn into_iter(self) -> Self::IntoIter {
188    Iterator1D {
189      steps: self,
190      index: 0,
191      index_back: self.steps(),
192    }
193  }
194}
195
196impl<T> IntoParallelIterator for Steps<T>
197where
198  T: Send
199    + Sync
200    + std::ops::Mul<f64, Output = T>
201    + std::ops::Add<T, Output = T>
202    + std::ops::Div<f64, Output = T>
203    + std::ops::Sub<T, Output = T>
204    + Copy,
205{
206  type Item = T;
207  type Iter = ParIterator1D<T>;
208
209  fn into_par_iter(self) -> Self::Iter {
210    ParIterator1D { steps: self }
211  }
212}
213
214/// Iterator for 1D steps
215pub struct Iterator1D<T> {
216  steps: Steps<T>,
217  index: usize,
218  index_back: usize,
219}
220
221impl<T> Iterator for Iterator1D<T>
222where
223  T: std::ops::Mul<f64, Output = T>
224    + std::ops::Add<T, Output = T>
225    + std::ops::Div<f64, Output = T>
226    + std::ops::Sub<T, Output = T>
227    + Copy,
228{
229  type Item = T;
230
231  fn next(&mut self) -> Option<Self::Item> {
232    if self.index >= self.index_back {
233      return None;
234    }
235
236    let value = self.steps.value(self.index);
237    self.index += 1;
238
239    Some(value)
240  }
241}
242
243impl<T> ExactSizeIterator for Iterator1D<T>
244where
245  T: std::ops::Mul<f64, Output = T>
246    + std::ops::Add<T, Output = T>
247    + std::ops::Div<f64, Output = T>
248    + std::ops::Sub<T, Output = T>
249    + Copy,
250{
251  fn len(&self) -> usize {
252    self.steps.len()
253  }
254}
255
256impl<T> DoubleEndedIterator for Iterator1D<T>
257where
258  T: std::ops::Mul<f64, Output = T>
259    + std::ops::Add<T, Output = T>
260    + std::ops::Div<f64, Output = T>
261    + std::ops::Sub<T, Output = T>
262    + Copy,
263{
264  fn next_back(&mut self) -> Option<Self::Item> {
265    if self.index_back <= self.index {
266      return None;
267    }
268
269    self.index_back -= 1;
270    Some(self.steps.value(self.index_back))
271  }
272}
273
274/// Parallel iterator for 1D steps
275pub struct ParIterator1D<T> {
276  steps: Steps<T>,
277}
278
279impl<T> rayon::iter::plumbing::Producer for ParIterator1D<T>
280where
281  T: Send
282    + Sync
283    + std::ops::Mul<f64, Output = T>
284    + std::ops::Add<T, Output = T>
285    + std::ops::Div<f64, Output = T>
286    + std::ops::Sub<T, Output = T>
287    + Copy,
288{
289  type Item = T;
290  type IntoIter = Iterator1D<T>;
291
292  fn into_iter(self) -> Self::IntoIter {
293    self.steps.into_iter()
294  }
295
296  fn split_at(self, index: usize) -> (Self, Self) {
297    let (start, end) = self.steps.range();
298    let s1 = self.steps.value(index - 1);
299    let s2 = self.steps.value(index);
300    (
301      ParIterator1D {
302        steps: Steps(start, s1, index),
303      },
304      ParIterator1D {
305        steps: Steps(s2, end, self.steps.steps() - index),
306      },
307    )
308  }
309}
310
311impl<T> rayon::iter::ParallelIterator for ParIterator1D<T>
312where
313  T: Send
314    + Sync
315    + std::ops::Mul<f64, Output = T>
316    + std::ops::Add<T, Output = T>
317    + std::ops::Div<f64, Output = T>
318    + std::ops::Sub<T, Output = T>
319    + Copy,
320{
321  type Item = T;
322
323  fn drive_unindexed<C>(self, consumer: C) -> C::Result
324  where
325    C: rayon::iter::plumbing::UnindexedConsumer<Self::Item>,
326  {
327    rayon::iter::plumbing::bridge(self, consumer)
328  }
329
330  fn opt_len(&self) -> Option<usize> {
331    Some(self.steps.steps())
332  }
333}
334
335impl<T> rayon::iter::IndexedParallelIterator for ParIterator1D<T>
336where
337  T: Send
338    + Sync
339    + std::ops::Mul<f64, Output = T>
340    + std::ops::Add<T, Output = T>
341    + std::ops::Div<f64, Output = T>
342    + std::ops::Sub<T, Output = T>
343    + Copy,
344{
345  fn len(&self) -> usize {
346    self.steps.steps()
347  }
348
349  fn drive<C>(self, consumer: C) -> C::Result
350  where
351    C: rayon::iter::plumbing::Consumer<Self::Item>,
352  {
353    rayon::iter::plumbing::bridge(self, consumer)
354  }
355
356  fn with_producer<CB: rayon::iter::plumbing::ProducerCallback<Self::Item>>(
357    self,
358    callback: CB,
359  ) -> CB::Output {
360    callback.callback(self)
361  }
362}
363
364/// Utility for creating evenly spaced steps between two endpoints in 2 dimensions
365///
366/// ## Example:
367/// ```
368/// use spdcalc::utils::Steps2D;
369///
370/// let arr : Vec<(f64, f64)> = Steps2D((0., 1., 100), (0., 100., 100)).into_iter().collect();
371/// assert_eq!(arr.len(), 100 * 100);
372/// assert!((arr[0].0 - 0.).abs() < 1e-12);
373/// assert!((arr[99].0 - 1.).abs() < 1e-12);
374/// assert!((arr[0].1 - 0.).abs() < 1e-12);
375/// assert!((arr[99].1 - 0.).abs() < 1e-12);
376/// assert!((arr[9999].0 - 1.).abs() < 1e-12);
377/// assert!((arr[9999].1 - 100.).abs() < 1e-12);
378/// ```
379#[derive(Debug, Copy, Clone, Hash, Serialize, Deserialize)]
380pub struct Steps2D<T>(pub (T, T, usize), pub (T, T, usize));
381
382impl<T> Steps2D<T>
383where
384  T: std::ops::Div<f64, Output = T>
385    + std::ops::Sub<T, Output = T>
386    + std::ops::Mul<f64, Output = T>
387    + std::ops::Add<T, Output = T>
388    + Copy,
389{
390  /// Create a new instance
391  pub fn new(x: (T, T, usize), y: (T, T, usize)) -> Self {
392    Self(x, y)
393  }
394
395  /// Number of divisions for each axis
396  pub fn divisions(&self) -> (usize, usize) {
397    (
398      Steps::from(self.0).divisions(),
399      Steps::from(self.1).divisions(),
400    )
401  }
402
403  /// Range of each axis
404  pub fn ranges(&self) -> ((T, T), (T, T)) {
405    ((self.0 .0, self.0 .1), (self.1 .0, self.1 .1))
406  }
407
408  /// Total number of steps
409  pub fn len(&self) -> usize {
410    (self.0).2 * (self.1).2
411  }
412
413  /// Is it empty?
414  pub fn is_empty(&self) -> bool {
415    self.len() == 0
416  }
417
418  /// Same number of steps for each axis?
419  pub fn is_square(&self) -> bool {
420    self.0 .2 == self.1 .2
421  }
422
423  /// Get the value at the given index
424  pub fn value(&self, index: usize) -> (T, T) {
425    let cols = self.0 .2;
426    let rows = self.1 .2;
427    let (nx, ny) = get_2d_indices(index, cols);
428    let xt = if cols > 1 {
429      (nx as f64) / ((cols - 1) as f64)
430    } else {
431      0.
432    };
433    let yt = if rows > 1 {
434      (ny as f64) / ((rows - 1) as f64)
435    } else {
436      0.
437    };
438    let x = lerp(self.0 .0, self.0 .1, xt);
439    let y = lerp(self.1 .0, self.1 .1, yt);
440    (x, y)
441  }
442
443  /// Swap x and y
444  pub fn swapped(self) -> Steps2D<T> {
445    Steps2D(self.1, self.0)
446  }
447
448  /// Get the width of the gap between each step.
449  ///
450  /// ## Example:
451  /// ```
452  /// use spdcalc::utils::Steps2D;
453  ///
454  /// let steps = Steps2D((0., 4., 3), (0., 8., 3)); // 3 steps: |0| --- |2| --- |4| and |0| --- |4| --- |8|
455  /// let dx = 2.;
456  /// let dy = 4.;
457  /// assert!((steps.division_widths().0 - dx).abs() < 1e-12);
458  /// assert!((steps.division_widths().1 - dy).abs() < 1e-12);
459  /// ```
460  pub fn division_widths(&self) -> (T, T) {
461    (
462      Steps::from(self.0).division_width(),
463      Steps::from(self.1).division_width(),
464    )
465  }
466}
467
468impl<T> IntoIterator for Steps2D<T>
469where
470  T: std::ops::Div<f64, Output = T>
471    + std::ops::Sub<T, Output = T>
472    + std::ops::Mul<f64, Output = T>
473    + std::ops::Add<T, Output = T>
474    + Copy,
475{
476  type Item = (T, T);
477  type IntoIter = Iterator2D<T>;
478
479  fn into_iter(self) -> Self::IntoIter {
480    Iterator2D::new(self)
481  }
482}
483
484/// get the 2d indices (column, row) from the linear index
485pub fn get_2d_indices(index: usize, cols: usize) -> (usize, usize) {
486  ((index % cols), (index / cols))
487}
488
489/// get the 1d index corresponding to a (col, row) of a 2d lattice
490pub fn get_1d_index(col: usize, row: usize, cols: usize) -> usize {
491  assert!(col < cols);
492  row * cols + col
493}
494
495/// An iterator that will iterate through rows and columns, giving you the
496/// coordinates at every iteration. Like a 2d linspace.
497///
498/// ## Example
499///
500/// ```
501/// use spdcalc::utils::{Steps2D, Iterator2D};
502///
503/// let x_steps = (0., 100., 11); // 0-100 in 10 steps = [0, 10, 20, ..., 100]
504/// let y_steps = (0., 10., 6); // 0-10 in 5 steps = [0, 2, 4, 6, 8, 10]
505/// let steps = Steps2D::new(x_steps, y_steps);
506/// let grid : Vec<f64> = Iterator2D::new(steps).map(|(x, y)| {
507///    x * y
508/// }).collect();
509/// assert_eq!(grid[12], 20.);
510/// ```
511#[derive(Debug, Copy, Clone)]
512pub struct Iterator2D<T>
513where
514  T: std::ops::Div<f64, Output = T>
515    + std::ops::Sub<T, Output = T>
516    + std::ops::Mul<f64, Output = T>
517    + std::ops::Add<T, Output = T>
518    + Copy,
519{
520  steps: Steps2D<T>,
521  partition: (usize, usize),
522  index: usize,
523  index_back: usize,
524}
525
526impl<T> Iterator2D<T>
527where
528  T: std::ops::Div<f64, Output = T>
529    + std::ops::Sub<T, Output = T>
530    + std::ops::Mul<f64, Output = T>
531    + std::ops::Add<T, Output = T>
532    + Copy,
533{
534  /// Create a new 2d iterator
535  pub fn new(steps: Steps2D<T>) -> Self {
536    Self::new_partition(steps, 0, steps.len())
537  }
538
539  /// Create a new 2d iterator over a portion of the grid
540  pub fn new_partition(steps: Steps2D<T>, start: usize, end: usize) -> Self {
541    assert!(
542      start <= end,
543      "start is greater than end, start {}, end {}",
544      start,
545      end
546    );
547    Self {
548      steps,
549      partition: (start, end),
550      index: start,
551      index_back: end,
552    }
553  }
554
555  /// get the x and y coordinates corresponding to the specified index
556  pub fn get_xy(&self, index: usize) -> (T, T) {
557    self.steps.value(index)
558  }
559
560  /// Get the x step size
561  pub fn get_dx(&self) -> T {
562    self.steps.division_widths().0
563  }
564
565  /// Get the y step size
566  pub fn get_dy(&self) -> T {
567    self.steps.division_widths().1
568  }
569}
570
571impl<T> rayon::iter::IntoParallelIterator for Steps2D<T>
572where
573  T: Send
574    + Sync
575    + std::ops::Div<f64, Output = T>
576    + std::ops::Sub<T, Output = T>
577    + std::ops::Mul<f64, Output = T>
578    + std::ops::Add<T, Output = T>
579    + Copy,
580{
581  type Item = (T, T);
582  type Iter = ParIterator2D<T>;
583  fn into_par_iter(self) -> Self::Iter {
584    ParIterator2D {
585      it: self.into_iter(),
586    }
587  }
588}
589
590impl<T> Iterator for Iterator2D<T>
591where
592  T: std::ops::Div<f64, Output = T>
593    + std::ops::Sub<T, Output = T>
594    + std::ops::Mul<f64, Output = T>
595    + std::ops::Add<T, Output = T>
596    + Copy,
597{
598  type Item = (T, T); // x, y
599
600  fn next(&mut self) -> Option<Self::Item> {
601    if self.index >= self.index_back {
602      return None;
603    }
604
605    let item = self.get_xy(self.index);
606    self.index += 1;
607    Some(item)
608  }
609}
610
611impl<T> DoubleEndedIterator for Iterator2D<T>
612where
613  T: std::ops::Div<f64, Output = T>
614    + std::ops::Sub<T, Output = T>
615    + std::ops::Mul<f64, Output = T>
616    + std::ops::Add<T, Output = T>
617    + Copy,
618{
619  fn next_back(&mut self) -> Option<Self::Item> {
620    if self.index_back <= self.index {
621      return None;
622    }
623
624    self.index_back -= 1;
625    let item = self.get_xy(self.index_back);
626    Some(item)
627  }
628}
629
630impl<T> ExactSizeIterator for Iterator2D<T>
631where
632  T: std::ops::Div<f64, Output = T>
633    + std::ops::Sub<T, Output = T>
634    + std::ops::Mul<f64, Output = T>
635    + std::ops::Add<T, Output = T>
636    + Copy,
637{
638  fn len(&self) -> usize {
639    self.partition.1 - self.partition.0
640  }
641}
642
643/// Parallel iterator for 2D grid
644pub struct ParIterator2D<T>
645where
646  T: Send
647    + Sync
648    + std::ops::Mul<f64, Output = T>
649    + std::ops::Add<T, Output = T>
650    + std::ops::Div<f64, Output = T>
651    + std::ops::Sub<T, Output = T>
652    + Copy,
653{
654  it: Iterator2D<T>,
655}
656
657impl<T> rayon::iter::plumbing::Producer for ParIterator2D<T>
658where
659  T: Send
660    + Sync
661    + std::ops::Mul<f64, Output = T>
662    + std::ops::Add<T, Output = T>
663    + std::ops::Div<f64, Output = T>
664    + std::ops::Sub<T, Output = T>
665    + Copy,
666{
667  type Item = (T, T);
668  type IntoIter = Iterator2D<T>;
669
670  fn into_iter(self) -> Self::IntoIter {
671    self.it
672  }
673
674  fn split_at(self, index: usize) -> (Self, Self) {
675    let left = Iterator2D::new_partition(
676      self.it.steps,
677      self.it.partition.0,
678      self.it.partition.0 + index,
679    );
680    let right = Iterator2D::new_partition(
681      self.it.steps,
682      self.it.partition.0 + index,
683      self.it.partition.1,
684    );
685    (ParIterator2D { it: left }, ParIterator2D { it: right })
686  }
687}
688
689impl<T> rayon::iter::ParallelIterator for ParIterator2D<T>
690where
691  T: Send
692    + Sync
693    + std::ops::Mul<f64, Output = T>
694    + std::ops::Add<T, Output = T>
695    + std::ops::Div<f64, Output = T>
696    + std::ops::Sub<T, Output = T>
697    + Copy,
698{
699  type Item = (T, T);
700
701  fn drive_unindexed<C>(self, consumer: C) -> C::Result
702  where
703    C: rayon::iter::plumbing::UnindexedConsumer<Self::Item>,
704  {
705    rayon::iter::plumbing::bridge(self, consumer)
706  }
707
708  fn opt_len(&self) -> Option<usize> {
709    Some(self.it.len())
710  }
711}
712
713impl<T> rayon::iter::IndexedParallelIterator for ParIterator2D<T>
714where
715  T: Send
716    + Sync
717    + std::ops::Mul<f64, Output = T>
718    + std::ops::Add<T, Output = T>
719    + std::ops::Div<f64, Output = T>
720    + std::ops::Sub<T, Output = T>
721    + Copy,
722{
723  fn len(&self) -> usize {
724    self.it.len()
725  }
726
727  fn drive<C>(self, consumer: C) -> C::Result
728  where
729    C: rayon::iter::plumbing::Consumer<Self::Item>,
730  {
731    rayon::iter::plumbing::bridge(self, consumer)
732  }
733
734  fn with_producer<CB: rayon::iter::plumbing::ProducerCallback<Self::Item>>(
735    self,
736    callback: CB,
737  ) -> CB::Output {
738    callback.callback(self)
739  }
740}
741
742/// Transpose a 2D matrix represented as a 1D vector.
743pub fn transpose_vec<T: Clone>(vec: Vec<T>, num_cols: usize) -> Vec<T> {
744  // use swap to transpose in place
745  let mut vec = vec;
746  let len = vec.len();
747  let num_rows = len.div_ceil(num_cols);
748  for row in 0..num_rows {
749    for col in (row + 1)..num_cols {
750      let index1 = get_1d_index(row, col, num_cols);
751      let index2 = get_1d_index(col, row, num_cols);
752      vec.swap(index1, index2);
753    }
754  }
755  vec
756}
757
758#[cfg(test)]
759pub(crate) mod testing {
760  use crate::dim::f64prefixes::*;
761  use crate::dim::ucum::*;
762  use crate::{Apodization, PeriodicPoling, SPDC};
763
764  pub fn percent_diff(actual: f64, expected: f64) -> f64 {
765    if expected == 0. && actual == 0. {
766      0.
767    } else {
768      200. * ((expected - actual).abs() / (expected + actual))
769    }
770  }
771
772  macro_rules! assert_nearly_equal {
773    ($name:expr, $actual:expr, $expected:expr, $accept_percent_diff:expr) => {
774      let diff = crate::utils::testing::percent_diff($actual, $expected);
775      assert!(
776        diff.abs() < $accept_percent_diff,
777        "{} percent difference: {}%\nactual: {}\nexpected: {}",
778        $name,
779        diff,
780        $actual,
781        $expected
782      );
783    };
784    ($name:expr, $actual:expr, $expected:expr) => {
785      assert_nearly_equal!($name, $actual, $expected, 1e-6);
786    };
787  }
788
789  pub(crate) use assert_nearly_equal;
790
791  pub fn testing_props(pp: bool) -> SPDC {
792    let mut spdc = SPDC::from_json(serde_json::json!({
793      "crystal": {
794        "kind": "BBO_1",
795        "pm_type": "Type2_e_eo",
796        "theta_deg": -3.0,
797        "phi_deg": 1.0,
798        "length_um": 2_000.0,
799        "temperature_c": 20.0,
800        "counter_propagation": false
801      },
802      "signal": {
803        "phi_deg": 15.0,
804        "theta_deg": 0.5,
805        "wavelength_nm": 1550.0,
806        "waist_um": 100.0
807      },
808      "pump": {
809        "wavelength_nm": 775.0,
810        "waist_um": 100.0,
811        "bandwidth_nm": 5.35,
812        "average_power_mw": 1,
813      },
814      "deff_pm_per_volt": 1
815    }))
816    .unwrap();
817    if pp {
818      spdc.pp = PeriodicPoling::new(28.5 * MICRO * M, Apodization::Off)
819    }
820    spdc
821  }
822}
823
824#[cfg(test)]
825mod tests {
826  use super::*;
827  #[test]
828  fn single_step_test() {
829    let actual: Vec<f64> = Steps(3.3, 4., 1).into_iter().collect();
830    let expected = vec![3.3];
831    assert_eq!(actual, expected);
832  }
833
834  #[test]
835  fn transpose_test() {
836    let actual: Vec<f64> = transpose_vec(Steps(1., 4., 4).into_iter().collect(), 2);
837    let expected = vec![1., 3., 2., 4.];
838    assert_eq!(actual, expected);
839
840    let actual = transpose_vec(
841      vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
842      4,
843    );
844    let expected = vec![1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15, 4, 8, 12, 16];
845
846    assert_eq!(actual, expected);
847  }
848
849  #[test]
850  fn iterator_2d_test() {
851    let it = Iterator2D::new(Steps2D::new((0., 1., 2), (0., 1., 2)));
852
853    let actual: Vec<(f64, f64)> = it.collect();
854    let expected = vec![(0., 0.), (1., 0.), (0., 1.), (1., 1.)];
855
856    assert!(
857      actual.iter().zip(expected.iter()).all(|(a, b)| a == b),
858      "assertion failed. actual: {:?}, expected: {:?}",
859      actual,
860      expected
861    );
862  }
863
864  #[test]
865  fn iterator_2d_rev_test() {
866    let it = Iterator2D::new(Steps2D::new((0., 1., 2), (0., 1., 2)));
867
868    let actual: Vec<(f64, f64)> = it.rev().collect();
869    let expected = vec![(1., 1.), (0., 1.), (1., 0.), (0., 0.)];
870
871    assert!(
872      actual.iter().zip(expected.iter()).all(|(a, b)| a == b),
873      "assertion failed. actual: {:?}, expected: {:?}",
874      actual,
875      expected
876    );
877  }
878
879  #[test]
880  fn test_split_at() {
881    let s = Steps(0., 0.9, 10);
882    use rayon::iter::plumbing::Producer;
883    let (a, b) = ParIterator1D { steps: s }.split_at(5);
884    let a: Vec<f64> = a.into_iter().collect();
885    let b: Vec<f64> = b.into_iter().collect();
886    assert_eq!(a, vec![0., 0.1, 0.2, 0.30000000000000004, 0.4]);
887    assert_eq!(b, vec![0.5, 0.6, 0.7, 0.8, 0.9]);
888  }
889
890  #[test]
891  fn test_par_iter_2d() {
892    use rayon::prelude::*;
893    let s = Steps2D::new((2., 3., 2), (2., 3., 2));
894    let actual: f64 = s.into_par_iter().map(|(x, y)| x + y).sum();
895    assert_eq!(actual, 20.);
896  }
897}
spdcalc/utils.rs

spdcalc/
utils.rs
