automatica 1.0.0

//! # Polar plot
//!
//! Polar plot returns the iterator providing the complex numbers at the given
//! angular frequencies.
//!
//! Functions use angular frequencies as default inputs.

use std::ops::{Add, Div, Sub};

use crate::{
    plots::Plotter, units::RadiansPerSecond, Atan2, Complex, Const, Floor, Hypot, Log, MulAdd, One,
    Pow, Zero,
};

/// Struct representing a Polar plot.
#[derive(Clone, Debug)]
pub struct Polar<T, U: Plotter<T>> {
    /// Transfer function
    tf: U,
    /// Minimum angular frequency of the plot
    min_freq: RadiansPerSecond<T>,
    /// Maximum angular frequency of the plot
    max_freq: RadiansPerSecond<T>,
    /// Step between frequencies
    step: T,
}

impl<T, U> Polar<T, U>
where
    T: MulAdd<Output = T> + PartialOrd + Zero,
    U: Plotter<T>,
{
    /// Create a `Polar` plot struct
    ///
    /// # Arguments
    ///
    /// * `tf` - Transfer function to plot
    /// * `min_freq` - Minimum angular frequency of the plot
    /// * `max_freq` - Maximum angular frequency of the plot
    /// * `step` - Step between frequencies
    ///
    /// `step` shall be in logarithmic scale. Use 0.1 to have 10 point per decade
    ///
    /// # Panics
    ///
    /// Panics if the step is not strictly positive of the minimum frequency
    /// is not lower than the maximum frequency.
    pub fn new(
        tf: U,
        min_freq: RadiansPerSecond<T>,
        max_freq: RadiansPerSecond<T>,
        step: T,
    ) -> Self {
        assert!(step > T::zero());
        assert!(min_freq < max_freq);

        Self {
            tf,
            min_freq,
            max_freq,
            step,
        }
    }
}

impl<T, U> Polar<T, U>
where
    T: Const + PartialOrd + MulAdd<Output = T> + Zero,
    U: Plotter<T>,
{
    /// Create a `Polar` plot struct
    ///
    /// # Arguments
    ///
    /// * `tf` - Transfer function to plot
    /// * `min_freq` - Minimum angular frequency of the plot
    /// * `step` - Step between frequencies
    ///
    /// `step` shall be in logarithmic scale. Use 0.1 to have 10 point per decade
    ///
    /// # Panics
    ///
    /// Panics if the step is not strictly positive of the minimum frequency
    /// is not lower than pi.
    pub fn new_discrete(tf: U, min_freq: RadiansPerSecond<T>, step: T) -> Self {
        let pi = RadiansPerSecond(T::pi());
        assert!(step > T::zero());
        assert!(min_freq < pi);

        Self {
            tf,
            min_freq,
            max_freq: pi,
            step,
        }
    }
}

impl<T, U> IntoIterator for Polar<T, U>
where
    T: Add<Output = T>
        + Clone
        + Div<Output = T>
        + Floor
        + From<f32>
        + Log
        + MulAdd<Output = T>
        + Sub<Output = T>
        + One
        + PartialOrd
        + Pow<T>
        + Zero,
    U: Plotter<T>,
{
    type Item = Data<T>;
    type IntoIter = IntoIter<T, U>;

    fn into_iter(self) -> Self::IntoIter {
        let min = self.min_freq.0.log10();
        let max = self.max_freq.0.log10();
        let intervals = ((max - min.clone()) / self.step.clone()).floor();
        Self::IntoIter {
            tf: self.tf,
            intervals,
            step: self.step,
            base_freq_exp: min,
            index: T::zero(),
        }
    }
}

/// Struct for the Polar plot data point iteration.
#[derive(Clone, Debug)]
pub struct IntoIter<T, U: Plotter<T>> {
    /// Transfer function
    tf: U,
    /// Number of intervals of the plot
    intervals: T,
    /// Step between frequencies
    step: T,
    /// Start frequency exponent
    base_freq_exp: T,
    /// Current data index
    index: T,
}

/// Struct to hold the data returned by the Polar iterator.
#[derive(Clone, Copy, Debug)]
pub struct Data<T> {
    /// Frequency
    freq: T,
    /// Output
    output: Complex<T>,
}

impl<T> Data<T>
where
    T: Atan2 + Clone + Hypot,
{
    /// Get the frequency
    pub fn freq(&self) -> T {
        self.freq.clone()
    }

    /// Get the output
    pub fn output(&self) -> Complex<T> {
        self.output.clone()
    }

    /// Get the real part
    pub fn real(&self) -> T {
        self.output.re.clone()
    }

    /// Get the imaginary part
    pub fn imag(&self) -> T {
        self.output.im.clone()
    }

    /// Get the magnitude
    pub fn magnitude(&self) -> T {
        self.output.norm()
    }

    /// Get the phase
    pub fn phase(&self) -> T {
        self.output.arg()
    }
}

/// Implementation of the Iterator trait for `Polar` struct
impl<T, U> Iterator for IntoIter<T, U>
where
    T: Add<Output = T> + Clone + From<f32> + MulAdd<Output = T> + One + PartialOrd + Pow<T>,
    U: Plotter<T>,
{
    type Item = Data<T>;

    fn next(&mut self) -> Option<Self::Item> {
        if self.index > self.intervals {
            None
        } else {
            let freq_exponent = MulAdd::mul_add(
                self.step.clone(),
                self.index.clone(),
                self.base_freq_exp.clone(),
            );
            // Casting is safe for both f32 and f64, representation is exact.
            let omega = T::from(10.0_f32).powf(freq_exponent);
            self.index = self.index.clone() + T::one();
            Some(Data {
                freq: omega.clone(),
                output: self.tf.eval_point(omega),
            })
        }
    }
}

#[cfg(test)]
mod tests {
    use crate::transfer_function::{continuous::Tf, discrete::Tfz};

    use super::*;

    #[test]
    fn create_iterator() {
        let tf = Tf::new([2., 3.], [1., 1., 1.]);
        let iter = Polar::new(tf, RadiansPerSecond(10.), RadiansPerSecond(1000.), 0.1).into_iter();
        assert_relative_eq!(20., iter.intervals);
        assert_relative_eq!(1., iter.base_freq_exp);
        assert_relative_eq!(0., iter.index);
    }

    #[test]
    fn create_discrete() {
        let tf = Tfz::new([2., 3.], [1., 1., 1.]);
        let iter = Polar::new_discrete(tf, RadiansPerSecond(0.01), 0.1).into_iter();
        assert!(iter.last().unwrap().freq() <= std::f32::consts::PI);
    }

    #[test]
    fn data_struct() {
        let c = Complex::new(3., 4.);
        let p = Data {
            freq: 1.,
            output: c,
        };
        assert_eq!(c, p.output());
        assert_relative_eq!(1., p.freq());
        assert_relative_eq!(3., p.real());
        assert_relative_eq!(4., p.imag());
        assert_relative_eq!(5., p.magnitude());
        assert_relative_eq!(0.9273, p.phase(), max_relative = 0.00001);
    }

    #[test]
    fn iterator() {
        let tf = Tf::new([2., 3.], [1., 1., 1.]);
        let iter = Polar::new(tf, RadiansPerSecond(10.), RadiansPerSecond(1000.), 0.1).into_iter();
        // 20 steps -> 21 iteration
        assert_eq!(21, iter.count());
    }
}