fdtd-coremem/crates/coremem/src/stim/time_varying.rs

//! time-varying portions of a Stimulus

use crate::real::{self, Real};
use crate::stim::FieldMags;

pub trait TimeVarying<R> {
    fn at(&self, t_sec: R) -> FieldMags<R>;
}

pub trait TimeVaryingExt<R>: Sized {
    fn shifted(self, new_start: R) -> Shifted<R, Self> {
        Shifted::new(self, new_start)
    }
    fn gated(self, from: R, to: R) -> Gated<R, Self> {
        Gated::new(Pulse::new(from, to), self)
    }
    fn scaled<T: TimeVarying<R>>(self, scale: T) -> Scaled<Self, T> {
        Scaled::new(self, scale)
    }
    fn summed<T: TimeVarying<R>>(self, with: T) -> Summed<Self, T> {
        Summed::new(self, with)
    }
}

impl<R, T> TimeVaryingExt<R> for T {}

impl<R: Real> TimeVarying<R> for FieldMags<R> {
    fn at(&self, _t_sec: R) -> FieldMags<R> {
        *self
    }
}

// assumed to represent the E field
impl<R: Real> TimeVarying<real::Finite<R>> for real::Finite<R> {
    fn at(&self, _t_sec: real::Finite<R>) -> FieldMags<real::Finite<R>> {
        FieldMags::new_e(*self)
    }
}
impl TimeVarying<f32> for f32 {
    fn at(&self, _t_sec: f32) -> FieldMags<f32> {
        FieldMags::new_e(*self)
    }
}
impl TimeVarying<f64> for f64 {
    fn at(&self, _t_sec: f64) -> FieldMags<f64> {
        FieldMags::new_e(*self)
    }
}

// Vec<T> at any `t_sec` behaves as the sum of all its components at that time.
impl<R: Real, T: TimeVarying<R>> TimeVarying<R> for Vec<T> {
    fn at(&self, t_sec: R) -> FieldMags<R> {
        self.iter().fold(FieldMags::default(), |acc, i| acc + i.at(t_sec))
    }
}

pub struct UnitEH;
impl<R: Real> TimeVarying<R> for UnitEH {
    fn at(&self, _t_sec: R) -> FieldMags<R> {
        FieldMags::new_eh(R::one(), R::one())
    }
}

/// E field which changes magnitude sinusoidally as a function of t
#[derive(Clone)]
pub struct Sinusoid<R> {
    omega: R,
}

impl<R: Real> Sinusoid<R> {
    pub fn new(freq: R) -> Self {
        Self {
            omega: freq * R::two_pi(),
        }
    }
    pub fn from_wavelength(lambda: R) -> Self {
        Self::new(lambda.inv())
    }
    pub fn freq(&self) -> R {
        self.omega / R::two_pi()
    }
    pub fn wavelength(&self) -> R {
        self.freq().inv()
    }
    pub fn one_cycle(self) -> Gated<R, Self> {
        let wl = self.wavelength();
        self.gated(R::zero(), wl)
    }
    pub fn half_cycle(self) -> Gated<R, Self> {
        let wl = self.wavelength();
        self.gated(R::zero(), R::half() * wl)
    }
}

impl<R: Real> TimeVarying<R> for Sinusoid<R> {
    fn at(&self, t_sec: R) -> FieldMags<R> {
        let v = (t_sec * self.omega).sin();
        FieldMags::new_eh(v, v)
    }
}

/// E field that decays exponentially over t.
/// zero for all t < 0
#[derive(Clone)]
pub struct Exp<R> {
    tau: R,
}

impl<R: Real + TimeVarying<R>> Exp<R> {
    pub fn new(half_life: R) -> Self {
        let tau = R::ln2()/half_life;
        Self { tau }
    }

    pub fn new_at(amp: R, start: R, half_life: R) -> Shifted<R, Gated<R, Scaled<Self, R>>> {
        Self::new(half_life).scaled(amp).gated(R::zero(), half_life*100f32.cast::<R>()).shifted(start)
    }
}

impl<R: Real> TimeVarying<R> for Exp<R> {
    fn at(&self, t_sec: R) -> FieldMags<R> {
        let a = if t_sec < R::zero() {
            // queries for very negative `t_sec` could cause `a` to explode
            // and IEEE 754 makes exp(LARGE) be infinity.
            // later, these queries are gated with a multiply-by-zero,
            // but 0 times INF is NaN.
            // so make this zero-valued before the moment of interest.
            // (an alternative would be to set it to 1.0).
            R::zero()
        } else {
            (t_sec * -self.tau).exp()
        };
        FieldMags::new_eh(a, a)
    }
}

/// pulses E=1.0 and H=1.0 over the provided duration.
/// this is used as a building block to gate some VectorField over a specific time.
#[derive(Clone)]
pub struct Pulse<R> {
    start: R,
    end: R,
}

impl<R> Pulse<R> {
    pub fn new(start: R, end: R) -> Self {
        Self { start, end }
    }
}
impl<R: Real> Pulse<R> {
    fn contains(&self, t: R) -> bool {
        t >= self.start && t < self.end
    }
}

impl<R: Real> TimeVarying<R> for Pulse<R> {
    fn at(&self, t: R) -> FieldMags<R> {
        if self.contains(t) {
            FieldMags::new_eh(R::one(), R::one())
        } else {
            FieldMags::new_eh(R::zero(), R::zero())
        }
    }
}


pub type Gated<R, T> = Scaled<Pulse<R>, T>;

#[derive(Clone)]
pub struct Shifted<R, T> {
    start_at: R,
    inner: T,
}

impl<R, T> Shifted<R, T> {
    pub fn new(inner: T, start_at: R) -> Self {
        Self { inner, start_at }
    }
}

impl<R: Real, T: TimeVarying<R>> TimeVarying<R> for Shifted<R, T> {
    fn at(&self, t_sec: R) -> FieldMags<R> {
        self.inner.at(t_sec - self.start_at)
    }
}


#[derive(Clone)]
pub struct Scaled<A, B>(A, B);

impl<A, B> Scaled<A, B> {
    pub fn new(a: A, b: B) -> Self {
        Self(a, b)
    }
}

impl<R: Real, A: TimeVarying<R>, B: TimeVarying<R>> TimeVarying<R> for Scaled<A, B> {
    fn at(&self, t_sec: R) -> FieldMags<R> {
        self.0.at(t_sec).elem_mul(self.1.at(t_sec))
    }
}

#[derive(Clone)]
pub struct Summed<A, B>(A, B);

impl<A, B> Summed<A, B> {
    pub fn new(a: A, b: B) -> Self {
        Self(a, b)
    }
}

impl<R: Real, A: TimeVarying<R>, B: TimeVarying<R>> TimeVarying<R> for Summed<A, B> {
    fn at(&self, t_sec: R) -> FieldMags<R> {
        self.0.at(t_sec) + self.1.at(t_sec)
    }
}


#[cfg(test)]
mod test {
    use super::*;

    macro_rules! assert_approx_eq {
        ($x:expr, $e:expr, $h:expr) => {
            let x = $x;
            let e = $e;
            let h = $h;
            let diff_e = (x.e - e).abs();
            assert!(diff_e <= 0.001, "{:?} != {:?}", x, e);
            let diff_h = (x.h - h).abs();
            assert!(diff_h <= 0.001, "{:?} != {:?}", x, h);
        }
    }
    
    #[test]
    fn sinusoid() {
        let s = Sinusoid::new(1000.0);
        assert_eq!(s.at(0.0), FieldMags::default());
        assert_approx_eq!(s.at(0.00025), 1.0, 1.0);
        assert_approx_eq!(s.at(0.00050), 0.0, 0.0);
        assert_approx_eq!(s.at(0.00075), -1.0, -1.0);
    }
    
    #[test]
    fn sinusoid_from_wavelength() {
        let s = Sinusoid::from_wavelength(0.001);
        assert_eq!(s.at(0.0), FieldMags::default());
        assert_approx_eq!(s.at(0.00025), 1.0, 1.0);
        assert_approx_eq!(s.at(0.00050), 0.0, 0.0);
        assert_approx_eq!(s.at(0.00075), -1.0, -1.0);
    }
}