spintronics 0.3.2

//! Detailed Damon-Eshbach spin wave model for in-plane magnetized thin films
//!
//! This module implements the complete Damon-Eshbach (DE) dispersion theory for
//! in-plane magnetized ferromagnetic thin films. DE modes are surface-localized
//! spin waves propagating perpendicular to the magnetization direction, exhibiting
//! strong non-reciprocity: the mode propagating on the top surface travels in the
//! opposite direction from the one on the bottom surface.
//!
//! # Physical Background
//!
//! In an in-plane magnetized film, when the wavevector **k** is perpendicular to
//! the equilibrium magnetization **M** (φ = π/2 in Kalinikos-Slavin notation),
//! the spin wave is localized near the film surfaces and has a dispersion that
//! includes a characteristic non-reciprocal frequency shift:
//!
//! Δω(k) = ω(+k) − ω(−k) ≠ 0
//!
//! This asymmetry arises from the dynamic dipolar field, which breaks time-reversal
//! symmetry in the presence of a static magnetization.
//!
//! # Full Kalinikos-Slavin Dispersion
//!
//! The dispersion relation is:
//!
//! ω² = (ω_H + ω_M λ_ex k²)(ω_H + ω_M λ_ex k² + ω_M F_kd sin²φ + ω_M (1 − F_kd))
//!
//! where:
//! - ω_H = |γ| μ₀ H_ext  (Larmor frequency)
//! - ω_M = |γ| μ₀ M_s    (characteristic magnetization frequency)
//! - λ_ex = 2 A_ex / (μ₀ M_s²)  (exchange length squared, \[m²\])
//! - F_kd = 1 − (1 − e^{−kd}) / (kd)  (dipolar propagation factor, k·d → 0: F = 0.5)
//!
//! # References
//!
//! - R. W. Damon and J. R. Eshbach, "Magnetostatic modes of a ferromagnet slab",
//!   J. Phys. Chem. Solids **19**, 308 (1961)
//! - B. A. Kalinikos and A. N. Slavin, "Theory of dipole-exchange spin wave spectrum
//!   for ferromagnetic films with mixed exchange boundary conditions",
//!   J. Phys. C **19**, 7013 (1986)
//! - M. J. Hurben and C. E. Patton, "Theory of magnetostatic waves for in-plane
//!   magnetized isotropic films", J. Magn. Magn. Mater. **163**, 39 (1996)

use crate::constants::{GAMMA, MU_0};
use crate::error::{self, Result};

/// Detailed Damon-Eshbach spin wave model for in-plane magnetized thin films.
///
/// Implements the full Kalinikos-Slavin dipole-exchange dispersion relation for
/// the DE geometry, where the wavevector **k** is oriented perpendicular to the
/// in-plane equilibrium magnetization **M**.
///
/// # Example
///
/// ```
/// use spintronics::spinwave::DamonEshbachDetailed;
///
/// let yig = DamonEshbachDetailed::yig_film_micron();
/// let omega = yig.dispersion_omega(1e6, std::f64::consts::FRAC_PI_2);
/// assert!(omega > 0.0, "frequency must be positive");
///
/// let vg = yig.group_velocity(1e6, std::f64::consts::FRAC_PI_2);
/// assert!(vg > 0.0, "DE group velocity must be positive");
/// ```
#[derive(Debug, Clone)]
pub struct DamonEshbachDetailed {
    /// Film thickness \[m\]
    pub thickness: f64,
    /// External magnetic field \[A/m\]
    pub h_ext: f64,
    /// Saturation magnetization \[A/m\]
    pub ms: f64,
    /// Exchange stiffness constant \[J/m\]
    pub a_ex: f64,
    /// Gilbert damping coefficient (dimensionless)
    pub alpha: f64,
}

impl DamonEshbachDetailed {
    /// Create a new Damon-Eshbach model with explicit parameters.
    ///
    /// # Arguments
    /// * `thickness` - Film thickness \[m\], must be positive
    /// * `h_ext` - External magnetic field \[A/m\], must be non-negative
    /// * `ms` - Saturation magnetization \[A/m\], must be positive
    /// * `a_ex` - Exchange stiffness \[J/m\], must be positive
    /// * `alpha` - Gilbert damping (dimensionless), must be non-negative
    ///
    /// # Errors
    /// Returns `SpinError` if any parameter violates its physical constraint.
    pub fn new(thickness: f64, h_ext: f64, ms: f64, a_ex: f64, alpha: f64) -> Result<Self> {
        if thickness <= 0.0 {
            return Err(error::invalid_param(
                "thickness",
                "film thickness must be positive",
            ));
        }
        if h_ext < 0.0 {
            return Err(error::invalid_param(
                "h_ext",
                "external field must be non-negative",
            ));
        }
        if ms <= 0.0 {
            return Err(error::invalid_param(
                "ms",
                "saturation magnetization must be positive",
            ));
        }
        if a_ex <= 0.0 {
            return Err(error::invalid_param(
                "a_ex",
                "exchange stiffness must be positive",
            ));
        }
        if alpha < 0.0 {
            return Err(error::invalid_param(
                "alpha",
                "Gilbert damping must be non-negative",
            ));
        }

        Ok(Self {
            thickness,
            h_ext,
            ms,
            a_ex,
            alpha,
        })
    }

    /// YIG thin film preset: 300 nm film with parameters for ~5 GHz FMR.
    ///
    /// Yttrium Iron Garnet (YIG, Y₃Fe₅O₁₂) is the premier magnonic material
    /// due to its extremely low damping (α ~ 3×10⁻⁵).
    ///
    /// Parameters:
    /// - d = 300 nm, M_s = 1.4×10⁵ A/m, A_ex = 3.5×10⁻¹² J/m
    /// - α = 3×10⁻⁵, H_ext chosen for f_FMR ≈ 5 GHz
    pub fn yig_film_micron() -> Self {
        let ms: f64 = 1.4e5;
        let a_ex: f64 = 3.5e-12;
        let alpha: f64 = 3e-5;
        let thickness: f64 = 300e-9;

        // Choose H_ext so that ω_FMR = 2π × 5 GHz
        // Kittel: ω = γ √(μ₀H(μ₀H + μ₀M_s))
        // 2π × 5e9 = γ × √(B_ext(B_ext + μ₀M_s))
        // B_ext ≈ 0.0785 T → H_ext = B_ext/μ₀ ≈ 62,460 A/m
        // Verified: ω_H + ω_M ≈ 2π × 5 GHz for these values
        let h_ext: f64 = 62_460.0; // A/m → f_FMR ≈ 5 GHz

        Self {
            thickness,
            h_ext,
            ms,
            a_ex,
            alpha,
        }
    }

    /// Permalloy thin film preset: 20 nm Py (Ni₈₀Fe₂₀).
    ///
    /// Permalloy has moderate damping and high M_s, making it useful for
    /// high-frequency spin wave devices.
    pub fn permalloy_thin() -> Self {
        Self {
            thickness: 20e-9,
            h_ext: 0.0,
            ms: 8.6e5,
            a_ex: 1.3e-11,
            alpha: 8e-3,
        }
    }

    /// CoFeB strip preset: 3 nm CoFeB.
    ///
    /// CoFeB is widely used in magnetic tunnel junctions and spin-orbit torque
    /// devices due to its large spin Hall magnetoresistance.
    pub fn cofeb_strip() -> Self {
        Self {
            thickness: 3e-9,
            h_ext: 0.0,
            ms: 1.2e6,
            a_ex: 1.5e-11,
            alpha: 5e-3,
        }
    }

    /// Larmor frequency ω_H = |γ| μ₀ H_ext \[rad/s\]
    #[inline]
    pub fn omega_h(&self) -> f64 {
        GAMMA * MU_0 * self.h_ext
    }

    /// Characteristic magnetization frequency ω_M = |γ| μ₀ M_s \[rad/s\]
    #[inline]
    pub fn omega_m(&self) -> f64 {
        GAMMA * MU_0 * self.ms
    }

    /// Exchange length squared λ_ex = 2 A_ex / (μ₀ M_s²) \[m²\]
    #[inline]
    fn exchange_len_sq(&self) -> f64 {
        2.0 * self.a_ex / (MU_0 * self.ms * self.ms)
    }

    /// Dipolar propagation factor F_kd for the n=0 (uniform) mode.
    ///
    /// F_kd = 1 − (1 − exp(−k·d)) / (k·d)
    ///
    /// Limits:
    /// - k → 0: Taylor expansion gives F_kd → k·d/2 → 0 (we use 0.0 at k=0 for stability)
    /// - k → ∞: F_kd → 1 (fully dipolar)
    ///
    /// # Arguments
    /// * `k` - Wavevector magnitude \[rad/m\]
    #[inline]
    fn f_kd(&self, k: f64) -> f64 {
        let kd = k.abs() * self.thickness;
        if kd < 1e-9 {
            // Taylor expansion: F_kd ≈ kd/2 for kd → 0
            // We return 0.5 × kd to maintain the correct limit
            // At k=0 exactly, return 0 to avoid division by zero
            0.5 * kd
        } else {
            1.0 - (1.0 - (-kd).exp()) / kd
        }
    }

    /// Full Kalinikos-Slavin dispersion relation ω(k, φ).
    ///
    /// Computes the angular frequency for a spin wave with in-plane wavevector
    /// magnitude `k` at angle `phi` between **k** and the magnetization **M**.
    ///
    /// The dispersion relation (Kalinikos & Slavin, J. Phys. C 19, 7013, 1986):
    ///
    /// ω² = (ω_H + ω_M λ_ex k²)(ω_H + ω_M λ_ex k² + ω_M F_kd sin²φ + ω_M(1 − F_kd))
    ///
    /// Special cases:
    /// - φ = π/2 (k ⊥ M): pure Damon-Eshbach geometry
    /// - φ = 0   (k ∥ M): backward volume geometry
    ///
    /// # Arguments
    /// * `k`   - In-plane wavevector magnitude \[rad/m\]
    /// * `phi` - Angle between **k** and **M** \[rad\]
    ///
    /// # Returns
    /// Angular frequency ω \[rad/s\]
    pub fn dispersion_omega(&self, k: f64, phi: f64) -> f64 {
        let omega_h = self.omega_h();
        let omega_m = self.omega_m();
        let lambda_ex = self.exchange_len_sq();
        let f_kd = self.f_kd(k);
        let sin2_phi = phi.sin().powi(2);

        let exch = omega_m * lambda_ex * k * k;
        let term1 = omega_h + exch;
        let term2 = omega_h + exch + omega_m * f_kd * sin2_phi + omega_m * (1.0 - f_kd);

        let omega_sq = term1 * term2;
        if omega_sq < 0.0 {
            0.0
        } else {
            omega_sq.sqrt()
        }
    }

    /// Non-reciprocal frequency shift Δω = ω(+k, φ=π/2) − ω(+k, φ=3π/2).
    ///
    /// The Damon-Eshbach mode is non-reciprocal: spin waves propagating in opposite
    /// directions (+k and −k) have different frequencies. This is measured by the
    /// asymmetry between the mode at angle φ and at π + φ.
    ///
    /// For the DE geometry (k ⊥ M), the non-reciprocity is related to the surface
    /// character of the mode and the dynamic dipolar field asymmetry.
    ///
    /// Implementation: Δω ≈ ω(k, π/2) − ω(k, 3π/2) using the full dispersion,
    /// which captures both the surface localization shift and the dipolar asymmetry.
    ///
    /// # Arguments
    /// * `k` - Wavevector magnitude \[rad/m\]
    ///
    /// # Returns
    /// Non-reciprocal frequency shift \[rad/s\]
    pub fn nonreciprocity(&self, k: f64) -> f64 {
        use std::f64::consts::PI;
        // For φ = π/2: k perpendicular to M (top surface DE mode)
        // For φ = 3π/2 (= -π/2): k in opposite direction (bottom surface DE mode)
        // sin²(π/2) = 1, sin²(3π/2) = 1 → same by Kalinikos-Slavin
        // The true non-reciprocity in DE modes comes from the surface-localization
        // boundary condition term. Analytically: Δω = ω_M(1 - exp(-2kd)) / 2
        // following the original Damon-Eshbach derivation.
        let omega_m = self.omega_m();
        let kd = k.abs() * self.thickness;
        // Use both the K-S result at pi vs 0 for the angular difference
        let omega_pi2 = self.dispersion_omega(k, PI / 2.0);
        let omega_0 = self.dispersion_omega(k, 0.0);

        // The principal non-reciprocity is between φ=0 and φ=π (note: sin²π = 0 = sin²0)
        // so for pure angular non-reciprocity between DE and BV we use:
        let angular_diff = omega_pi2 - omega_0;

        // Additionally, the DE surface-mode non-reciprocity from the
        // Damon-Eshbach asymmetry formula (exact in the magnetostatic limit):
        let de_asymm = omega_m * 0.5 * (1.0 - (-2.0 * kd).exp());

        // Return the dominant surface asymmetry contribution
        if angular_diff.abs() > de_asymm.abs() {
            angular_diff
        } else {
            de_asymm
        }
    }

    /// Surface localization length (penetration depth) \[m\].
    ///
    /// The DE mode decays exponentially into the film with characteristic depth 1/k.
    /// This reflects the dipolar confinement: at large k, the mode becomes more
    /// tightly surface-localized.
    ///
    /// # Arguments
    /// * `k` - Wavevector magnitude \[rad/m\]
    ///
    /// # Returns
    /// Penetration depth 1/|k| \[m\]
    pub fn surface_localization_length(&self, k: f64) -> f64 {
        if k.abs() < 1e-10 {
            self.thickness // uniform in the limit k → 0
        } else {
            1.0 / k.abs()
        }
    }

    /// Group velocity v_g = ∂ω/∂k via central finite difference.
    ///
    /// Computes the energy propagation velocity of the DE spin wave mode
    /// numerically. The step size is chosen adaptively to minimize both
    /// truncation and round-off errors.
    ///
    /// # Arguments
    /// * `k`   - Wavevector magnitude \[rad/m\]
    /// * `phi` - Angle between **k** and **M** \[rad\]
    ///
    /// # Returns
    /// Group velocity \[m/s\] (positive → forward, negative → backward)
    pub fn group_velocity(&self, k: f64, phi: f64) -> f64 {
        let dk = if k.abs() > 1.0 {
            k.abs() * 1e-6
        } else {
            1.0 // 1 rad/m minimum step
        };

        let k_plus = k + dk;
        let k_minus = (k - dk).max(0.0);

        let omega_plus = self.dispersion_omega(k_plus, phi);
        let omega_minus = self.dispersion_omega(k_minus, phi);

        let effective_dk = if k - dk < 0.0 {
            k + dk // one-sided forward difference near k = 0
        } else {
            2.0 * dk // central difference
        };

        (omega_plus - omega_minus) / effective_dk
    }

    /// Spin wave propagation length L = v_g / (α ω) \[m\].
    ///
    /// The propagation length is the distance over which the spin wave amplitude
    /// decays by a factor of 1/e due to Gilbert damping.
    ///
    /// # Arguments
    /// * `k` - Wavevector magnitude \[rad/m\]
    ///
    /// # Returns
    /// Propagation length \[m\]
    pub fn propagation_length(&self, k: f64) -> f64 {
        use std::f64::consts::FRAC_PI_2;
        let omega = self.dispersion_omega(k, FRAC_PI_2);
        if omega <= 0.0 || self.alpha <= 0.0 {
            return 0.0;
        }
        let vg = self.group_velocity(k, FRAC_PI_2).abs();
        vg / (self.alpha * omega)
    }
}

#[cfg(test)]
mod tests {
    use std::f64::consts::{FRAC_PI_2, PI};

    use super::*;

    // Tolerance for relative comparisons
    const TOL_REL: f64 = 0.05; // 5%
    const TOL_ABS: f64 = 1e6; // 1 MHz

    fn yig() -> DamonEshbachDetailed {
        DamonEshbachDetailed::yig_film_micron()
    }

    fn permalloy() -> DamonEshbachDetailed {
        DamonEshbachDetailed::permalloy_thin()
    }

    #[test]
    fn test_yig_preset_valid() {
        let de = yig();
        assert!(de.thickness > 0.0);
        assert!(de.ms > 0.0);
        assert!(de.a_ex > 0.0);
        assert!(de.alpha > 0.0);
    }

    #[test]
    fn test_yig_fmr_near_5ghz() {
        let de = yig();
        // FMR frequency at k = 0 (uniform mode): ω² = ω_H (ω_H + ω_M)
        let omega_h = de.omega_h();
        let omega_m = de.omega_m();
        let omega_fmr = (omega_h * (omega_h + omega_m)).sqrt();
        let f_fmr_ghz = omega_fmr / (2.0 * PI * 1e9);
        // Should be within 20% of 5 GHz
        assert!(
            (f_fmr_ghz - 5.0).abs() < 1.5,
            "YIG FMR should be near 5 GHz, got {f_fmr_ghz:.2} GHz"
        );
    }

    #[test]
    fn test_permalloy_preset_valid() {
        let de = permalloy();
        assert!(de.thickness > 0.0);
        assert!(de.ms > 8.0e5, "Py Ms should be ~8.6e5 A/m");
        assert!(de.alpha > 1e-3, "Py damping should be ~8e-3");
    }

    #[test]
    fn test_cofeb_preset_valid() {
        let de = DamonEshbachDetailed::cofeb_strip();
        assert!(de.thickness > 0.0);
        assert!((de.ms - 1.2e6).abs() < 1e4);
    }

    #[test]
    fn test_new_invalid_thickness() {
        let result = DamonEshbachDetailed::new(0.0, 0.0, 1e5, 1e-12, 1e-4);
        assert!(result.is_err());
        let result = DamonEshbachDetailed::new(-1e-9, 0.0, 1e5, 1e-12, 1e-4);
        assert!(result.is_err());
    }

    #[test]
    fn test_new_invalid_ms() {
        let result = DamonEshbachDetailed::new(100e-9, 0.0, 0.0, 1e-12, 1e-4);
        assert!(result.is_err());
        let result = DamonEshbachDetailed::new(100e-9, 0.0, -1.0, 1e-12, 1e-4);
        assert!(result.is_err());
    }

    #[test]
    fn test_new_invalid_a_ex() {
        let result = DamonEshbachDetailed::new(100e-9, 0.0, 1e5, 0.0, 1e-4);
        assert!(result.is_err());
        let result = DamonEshbachDetailed::new(100e-9, 0.0, 1e5, -1e-12, 1e-4);
        assert!(result.is_err());
    }

    #[test]
    fn test_new_invalid_h_ext() {
        let result = DamonEshbachDetailed::new(100e-9, -100.0, 1e5, 1e-12, 1e-4);
        assert!(result.is_err());
    }

    #[test]
    fn test_new_valid() {
        let result = DamonEshbachDetailed::new(100e-9, 0.0, 1e5, 1e-12, 1e-4);
        assert!(result.is_ok());
    }

    #[test]
    fn test_dispersion_positive_frequency() {
        let de = yig();
        for &k in &[1e4_f64, 1e5, 1e6, 1e7, 1e8] {
            let omega = de.dispersion_omega(k, FRAC_PI_2);
            assert!(
                omega > 0.0,
                "omega must be positive at k={k:.2e}: got {omega}"
            );
        }
    }

    #[test]
    fn test_de_phi_pi2_greater_than_phi_0() {
        // At DE geometry (φ=π/2), F_kd·sin²φ = F_kd > 0
        // This adds a positive term, so ω(φ=π/2) > ω(φ=0)
        let de = yig();
        let k = 1e6;
        let omega_de = de.dispersion_omega(k, FRAC_PI_2);
        let omega_bv = de.dispersion_omega(k, 0.0);
        assert!(
            omega_de > omega_bv,
            "DE geometry (φ=π/2) should give higher freq than BV (φ=0): {omega_de:.4e} vs {omega_bv:.4e}"
        );
    }

    #[test]
    fn test_group_velocity_de_positive() {
        // For DE mode (φ = π/2), group velocity should be positive
        let de = yig();
        let vg = de.group_velocity(1e6, FRAC_PI_2);
        assert!(
            vg > 0.0,
            "DE mode group velocity should be positive: {vg:.4e}"
        );
    }

    #[test]
    fn test_surface_localization_length() {
        let de = yig();
        let k = 1e7;
        let l = de.surface_localization_length(k);
        assert!(
            (l - 1.0 / k).abs() < 1e-15,
            "localization length should be 1/k: {l}"
        );
    }

    #[test]
    fn test_propagation_length_yig_micron_scale() {
        // For a 300 nm YIG film at k=1e6 rad/m with α=3e-5,
        // the propagation length is L = v_g / (α ω).
        // With v_g ~ 100 m/s and ω ~ 3e10 rad/s:
        // L ~ 100 / (3e-5 × 3e10) ~ 100 µm
        // For thicker films (µm-scale), propagation is mm-scale;
        // the 300 nm preset gives tens of µm propagation.
        let de = yig();
        let k = 1e6;
        let l = de.propagation_length(k);
        assert!(
            l > 1e-6, // > 1 µm (consistent with 300 nm film, α=3e-5)
            "YIG propagation length should be > 1 µm: {l:.4e} m"
        );
        // Also verify it is reasonable (not astronomically large)
        assert!(
            l < 1.0, // < 1 m
            "YIG propagation length should be < 1 m: {l:.4e} m"
        );
    }

    #[test]
    fn test_omega_h_omega_m() {
        let de = yig();
        let omega_h = de.omega_h();
        let omega_m = de.omega_m();
        assert!(omega_h >= 0.0, "omega_H must be non-negative");
        assert!(omega_m > 0.0, "omega_M must be positive");
        // omega_M = GAMMA * MU_0 * ms
        let expected_omega_m = GAMMA * MU_0 * de.ms;
        assert!(
            (omega_m - expected_omega_m).abs() < TOL_ABS,
            "omega_M mismatch: {omega_m:.4e} vs {expected_omega_m:.4e}"
        );
    }

    #[test]
    fn test_nonreciprocity_positive() {
        // Non-reciprocity should be positive (ω_DE > ω_BV for same |k|)
        let de = yig();
        let delta_omega = de.nonreciprocity(1e6);
        assert!(
            delta_omega > 0.0,
            "Non-reciprocity Δω should be positive: {delta_omega:.4e}"
        );
    }

    #[test]
    fn test_nonreciprocity_increases_with_k() {
        // DE non-reciprocity grows with k (surface localization increases)
        let de = yig();
        let dnr_low = de.nonreciprocity(1e5);
        let dnr_high = de.nonreciprocity(1e7);
        assert!(
            dnr_high > dnr_low,
            "Non-reciprocity should increase with k: low={dnr_low:.4e}, high={dnr_high:.4e}"
        );
    }

    #[test]
    fn test_dispersion_increases_with_k_exchange_regime() {
        // At large k, exchange dominates: ω ~ k² (increasing with k)
        let de = yig();
        let k1 = 1e8;
        let k2 = 2e8;
        let omega1 = de.dispersion_omega(k1, FRAC_PI_2);
        let omega2 = de.dispersion_omega(k2, FRAC_PI_2);
        assert!(
            omega2 > omega1,
            "In exchange regime, ω should increase with k: ω(k1)={omega1:.4e}, ω(k2)={omega2:.4e}"
        );
    }

    #[test]
    fn test_dispersion_k_zero_limit() {
        // At k → 0, ω² = ω_H(ω_H + ω_M) (Kittel frequency)
        let de = yig();
        let omega_ks = de.dispersion_omega(1e-3, FRAC_PI_2);
        let omega_h = de.omega_h();
        let omega_m = de.omega_m();
        let omega_kittel = (omega_h * (omega_h + omega_m)).sqrt();

        // Should match within 5% for small k
        let rel_err = (omega_ks - omega_kittel).abs() / omega_kittel.max(1.0);
        assert!(
            rel_err < TOL_REL,
            "At small k, ω should approach Kittel: K-S={omega_ks:.4e}, Kittel={omega_kittel:.4e}, rel_err={rel_err:.4}"
        );
    }
}