math_audio_xem_common/

source.rs

//! Sound source definitions for room acoustics simulations

use crate::types::Point3D;
use ndarray::Array2;
use serde::{Deserialize, Serialize};

/// Directivity pattern sampled on a grid
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct DirectivityGrid {
    /// Horizontal angles (azimuth) in degrees [0, 360) with step 10°
    pub horizontal_angles: Vec<f64>,
    /// Vertical angles (elevation) in degrees [0, 180] with step 10°
    pub vertical_angles: Vec<f64>,
    /// Magnitude at each (horizontal, vertical) angle pair
    /// Shape: [n_vertical, n_horizontal]
    pub magnitude: Array2<f64>,
}

impl DirectivityGrid {
    /// Create omnidirectional pattern (uniform radiation)
    pub fn omnidirectional() -> Self {
        let horizontal_angles: Vec<f64> = (0..36).map(|i| i as f64 * 10.0).collect();
        let vertical_angles: Vec<f64> = (0..19).map(|i| i as f64 * 10.0).collect();

        let magnitude = Array2::ones((vertical_angles.len(), horizontal_angles.len()));

        Self {
            horizontal_angles,
            vertical_angles,
            magnitude,
        }
    }

    /// Create a cardioid directivity pattern
    pub fn cardioid() -> Self {
        let horizontal_angles: Vec<f64> = (0..36).map(|i| i as f64 * 10.0).collect();
        let vertical_angles: Vec<f64> = (0..19).map(|i| i as f64 * 10.0).collect();

        let mut magnitude = Array2::zeros((vertical_angles.len(), horizontal_angles.len()));

        for (v_idx, &v_angle) in vertical_angles.iter().enumerate() {
            for (h_idx, &h_angle) in horizontal_angles.iter().enumerate() {
                // Cardioid pattern: 0.5 * (1 + cos(theta))
                let theta_rad = v_angle.to_radians();
                let phi_rad = h_angle.to_radians();
                // Forward direction is along +Y (phi=90, theta=90)
                let forward_dot = theta_rad.sin() * phi_rad.sin();
                magnitude[[v_idx, h_idx]] = 0.5 * (1.0 + forward_dot).max(0.0);
            }
        }

        Self {
            horizontal_angles,
            vertical_angles,
            magnitude,
        }
    }

    /// Interpolate directivity at arbitrary direction
    pub fn interpolate(&self, theta: f64, phi: f64) -> f64 {
        // Convert spherical to degrees
        let theta_deg = theta.to_degrees();
        let mut phi_deg = phi.to_degrees();

        // Normalize phi to [0, 360)
        while phi_deg < 0.0 {
            phi_deg += 360.0;
        }
        while phi_deg >= 360.0 {
            phi_deg -= 360.0;
        }

        // Find surrounding angles
        let h_idx = (phi_deg / 10.0).floor() as usize;
        let v_idx = (theta_deg / 10.0).floor() as usize;

        let h_idx = h_idx.min(self.horizontal_angles.len() - 1);
        let v_idx = v_idx.min(self.vertical_angles.len() - 1);

        let h_next = (h_idx + 1) % self.horizontal_angles.len();
        let v_next = (v_idx + 1).min(self.vertical_angles.len() - 1);

        // Bilinear interpolation
        let h_frac = (phi_deg / 10.0) - h_idx as f64;
        let v_frac = (theta_deg / 10.0) - v_idx as f64;

        let m00 = self.magnitude[[v_idx, h_idx]];
        let m01 = self.magnitude[[v_idx, h_next]];
        let m10 = self.magnitude[[v_next, h_idx]];
        let m11 = self.magnitude[[v_next, h_next]];

        let m0 = m00 * (1.0 - h_frac) + m01 * h_frac;
        let m1 = m10 * (1.0 - h_frac) + m11 * h_frac;

        m0 * (1.0 - v_frac) + m1 * v_frac
    }
}

/// Directivity model
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum Directivity {
    /// Sampled grid (static)
    Grid(DirectivityGrid),
    /// Classical frequency-dependent directional source
    Classical {
        /// Horizontal beamwidth (degrees) at high frequency
        h_angle: f64,
        /// Vertical beamwidth (degrees) at high frequency
        v_angle: f64,
    },
}

impl Directivity {
    /// Get amplitude scaling factor for a direction and frequency
    pub fn amplitude(&self, theta: f64, phi: f64, frequency: f64) -> f64 {
        match self {
            Directivity::Grid(grid) => grid.interpolate(theta, phi),
            Directivity::Classical { h_angle, v_angle } => {
                // Classical pattern: Gaussian-like beam narrowing with frequency
                // Transition from Omni (360) at 80Hz to Target at 800Hz (approx)

                let min_f = 80.0;
                let max_f = 800.0; // Transition end

                let t = if frequency <= min_f {
                    0.0
                } else if frequency >= max_f {
                    1.0
                } else {
                    (frequency.ln() - min_f.ln()) / (max_f.ln() - min_f.ln())
                };

                let h_width = 360.0 * (1.0 - t) + h_angle * t;
                let v_width = 360.0 * (1.0 - t) + v_angle * t;

                // Direction relative to forward axis (+Y: theta=90, phi=90)
                // We need angle from axis.
                // Axis vector: (0, 1, 0)
                // Point vector: (sin(theta)cos(phi), sin(theta)sin(phi), cos(theta))
                // Dot product = sin(theta)sin(phi)
                // Angle off axis 'alpha' = acos(dot)

                // But we have separate H and V widths.
                // H angle is deviation in XY plane (phi). Forward is 90.
                // V angle is deviation in Z plane (theta). Forward is 90.

                let mut phi_deg = phi.to_degrees();
                while phi_deg < 0.0 {
                    phi_deg += 360.0;
                }
                let delta_phi = (phi_deg - 90.0).abs().min((360.0 - (phi_deg - 90.0)).abs()); // Shortest distance to 90

                let theta_deg = theta.to_degrees();
                let delta_theta = (theta_deg - 90.0).abs();

                // Normalized deviations
                // x = angle / (width/2)
                let x_h = delta_phi / (h_width / 2.0);
                let x_v = delta_theta / (v_width / 2.0);

                // Combined distance in parameter space (elliptical cone)
                let r_sq = x_h * x_h + x_v * x_v;

                // Gaussian: m = 0.5^(r^2) -> -6dB at edge (r=1)
                0.5_f64.powf(r_sq)
            }
        }
    }
}

/// Crossover filter for frequency-limited sources
#[derive(Debug, Clone, Default, Serialize, Deserialize)]
pub enum CrossoverFilter {
    /// Full range (no filter)
    #[default]
    FullRange,
    /// Low-pass filter (Butterworth)
    Lowpass {
        /// Cutoff frequency (Hz)
        cutoff_freq: f64,
        /// Filter order (e.g., 2, 4)
        order: u32,
    },
    /// High-pass filter (Butterworth)
    Highpass {
        /// Cutoff frequency (Hz)
        cutoff_freq: f64,
        /// Filter order (e.g., 2, 4)
        order: u32,
    },
    /// Band-pass filter (combined Lowpass and Highpass)
    Bandpass {
        /// Low cutoff frequency (Hz)
        low_cutoff: f64,
        /// High cutoff frequency (Hz)
        high_cutoff: f64,
        /// Filter order
        order: u32,
    },
}

impl CrossoverFilter {
    /// Get amplitude multiplier at a given frequency
    pub fn amplitude_at_frequency(&self, frequency: f64) -> f64 {
        match self {
            CrossoverFilter::FullRange => 1.0,
            CrossoverFilter::Lowpass { cutoff_freq, order } => {
                let ratio = frequency / cutoff_freq;
                1.0 / (1.0 + ratio.powi(*order as i32 * 2)).sqrt()
            }
            CrossoverFilter::Highpass { cutoff_freq, order } => {
                let ratio = cutoff_freq / frequency;
                1.0 / (1.0 + ratio.powi(*order as i32 * 2)).sqrt()
            }
            CrossoverFilter::Bandpass {
                low_cutoff,
                high_cutoff,
                order,
            } => {
                let high_ratio = low_cutoff / frequency;
                let low_ratio = frequency / high_cutoff;
                let hp_response = 1.0 / (1.0 + high_ratio.powi(*order as i32 * 2)).sqrt();
                let lp_response = 1.0 / (1.0 + low_ratio.powi(*order as i32 * 2)).sqrt();
                hp_response * lp_response
            }
        }
    }
}

/// Sound source with position and directivity
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Source {
    /// Source position
    pub position: Point3D,
    /// Directivity pattern
    pub directivity: Directivity,
    /// Source amplitude (strength)
    pub amplitude: f64,
    /// Crossover filter
    pub crossover: CrossoverFilter,
    /// Source name (e.g., "Left Main", "Subwoofer")
    pub name: String,
}

impl Source {
    /// Create a new source
    pub fn new(position: Point3D, directivity: Directivity, amplitude: f64) -> Self {
        Self {
            position,
            directivity,
            amplitude,
            crossover: CrossoverFilter::FullRange,
            name: String::from("Source"),
        }
    }

    /// Create a simple omnidirectional source at position
    pub fn omnidirectional(position: Point3D, amplitude: f64) -> Self {
        Self::new(
            position,
            Directivity::Grid(DirectivityGrid::omnidirectional()),
            amplitude,
        )
    }

    /// Create a classical directional source
    pub fn classical(position: Point3D, h_angle: f64, v_angle: f64, amplitude: f64) -> Self {
        Self::new(
            position,
            Directivity::Classical { h_angle, v_angle },
            amplitude,
        )
    }

    /// Set crossover filter
    pub fn with_crossover(mut self, crossover: CrossoverFilter) -> Self {
        self.crossover = crossover;
        self
    }

    /// Set source name
    pub fn with_name(mut self, name: String) -> Self {
        self.name = name;
        self
    }

    /// Get directional amplitude towards a point at a specific frequency
    pub fn amplitude_towards(&self, point: &Point3D, frequency: f64) -> f64 {
        let dx = point.x - self.position.x;
        let dy = point.y - self.position.y;
        let dz = point.z - self.position.z;

        let r = (dx * dx + dy * dy + dz * dz).sqrt();
        if r < 1e-10 {
            return self.amplitude * self.crossover.amplitude_at_frequency(frequency);
        }

        let theta = (dz / r).acos();
        let phi = dy.atan2(dx);

        let directivity_factor = self.directivity.amplitude(theta, phi, frequency);
        let crossover_factor = self.crossover.amplitude_at_frequency(frequency);
        self.amplitude * directivity_factor * crossover_factor
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f64::consts::PI;

    #[test]
    fn test_omnidirectional_pattern() {
        let grid = DirectivityGrid::omnidirectional();
        // Should be 1.0 in all directions
        assert!((grid.interpolate(0.0, 0.0) - 1.0).abs() < 1e-6);
        assert!((grid.interpolate(PI / 2.0, PI) - 1.0).abs() < 1e-6);
        assert!((grid.interpolate(PI, 0.0) - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_crossover_lowpass() {
        let crossover = CrossoverFilter::Lowpass {
            cutoff_freq: 100.0,
            order: 2,
        };
        // Well below cutoff
        assert!((crossover.amplitude_at_frequency(10.0) - 1.0).abs() < 0.1);
        // At cutoff (-3dB point)
        let at_cutoff = crossover.amplitude_at_frequency(100.0);
        assert!(at_cutoff > 0.6 && at_cutoff < 0.8);
        // Well above cutoff
        assert!(crossover.amplitude_at_frequency(1000.0) < 0.1);
    }

    #[test]
    fn test_source_amplitude() {
        let source = Source::omnidirectional(Point3D::new(0.0, 0.0, 0.0), 1.0);
        let amp = source.amplitude_towards(&Point3D::new(1.0, 0.0, 0.0), 1000.0);
        assert!((amp - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_classical_source_beaming() {
        let source = Source::classical(Point3D::new(0.0, 0.0, 0.0), 60.0, 40.0, 1.0);

        // At low freq (50Hz), should be omni
        let amp_low = source.amplitude_towards(&Point3D::new(0.0, -1.0, 0.0), 50.0); // Rear (180 deg)
        // Forward is +Y. Rear is -Y.
        // Forward: theta=90, phi=90.
        // Rear: theta=90, phi=270 (-90).
        // Delta phi = 180.
        // If omni, width=360. x = 180/180 = 1. m = 0.5.
        // Wait, omni means width=360 implies +/- 180 deg?
        // Width is full width. +/- 180 covers full circle.
        // x = delta / 180.
        // Rear delta = 180. x=1. m=0.5 (-6dB).
        // So Omni via this formula is -6dB at rear.
        // True Omni should be 1.0 everywhere.
        // My Gaussian model with width=360 is not perfectly Omni.
        // But "Typically omni below 80Hz" usually allows some directivity or just "wide".
        // If I want true Omni, I should handle width >= 360 specially or increase width to infinity.
        // But 360 width (-6dB at +/-180) is acceptable for "Omni-like".
        assert!(amp_low > 0.4);

        // At high freq (2000Hz), should be narrow
        // H=60 -> +/- 30 deg is -6dB.
        // Side (90 deg off axis): x = 90/30 = 3. m = 0.5^9 = 0.0019.
        let amp_high_side = source.amplitude_towards(&Point3D::new(1.0, 0.0, 0.0), 2000.0); // Side (X axis, phi=0, delta=90)
        assert!(amp_high_side < 0.1);

        // On axis should be 1.0
        let amp_high_on = source.amplitude_towards(&Point3D::new(0.0, 1.0, 0.0), 2000.0); // Forward (+Y)
        assert!((amp_high_on - 1.0).abs() < 1e-6);
    }
}