oxiphoton 0.1.1

/// THz generation mechanisms:
///   - Optical rectification (OR) in nonlinear crystals via difference-frequency mixing
///   - Photoconductive antenna (PCA) emitters
///   - Free-electron laser (FEL) THz sources
#[cfg(test)]
use crate::error::OxiPhotonError;

// ─── Physical constants ────────────────────────────────────────────────────
const C0: f64 = 2.997_924_58e8; // m/s  (speed of light in vacuum)
const EPS0: f64 = 8.854_187_817e-12; // F/m
#[allow(dead_code)]
const H_PLANCK: f64 = 6.626_070_15e-34; // J·s
const E_CHARGE: f64 = 1.602_176_634e-19; // C

// ─── Optical Rectification ────────────────────────────────────────────────

/// Optical rectification THz source — difference-frequency mixing of a
/// femtosecond pump pulse inside a χ⁽²⁾ nonlinear crystal.
///
/// The THz field is generated by the second-order polarisation driven by
/// the intensity envelope of the ultrashort pump.
#[derive(Debug, Clone)]
pub struct OpticalRectification {
    /// Human-readable crystal identifier (e.g. "ZnTe", "LiNbO3").
    pub crystal_name: String,
    /// Effective second-order nonlinear coefficient d_eff (pm V⁻¹).
    pub d_eff_pm_per_v: f64,
    /// Crystal interaction length (mm).
    pub crystal_length_mm: f64,
    /// Pump centre wavelength (nm).
    pub pump_wavelength_nm: f64,
    /// Pump pulse duration FWHM (fs).
    pub pump_pulse_duration_fs: f64,
    /// Pump pulse energy (μJ).
    pub pump_energy_uj: f64,
    /// Real refractive index at the optical pump wavelength.
    pub refractive_index_optical: f64,
    /// Real refractive index at THz frequencies.
    pub refractive_index_thz: f64,
}

impl OpticalRectification {
    /// General constructor.
    #[allow(clippy::too_many_arguments)]
    pub fn new(
        crystal_name: impl Into<String>,
        d_eff_pm_per_v: f64,
        crystal_length_mm: f64,
        pump_wavelength_nm: f64,
        pump_pulse_duration_fs: f64,
        pump_energy_uj: f64,
        refractive_index_optical: f64,
        refractive_index_thz: f64,
    ) -> Self {
        Self {
            crystal_name: crystal_name.into(),
            d_eff_pm_per_v,
            crystal_length_mm,
            pump_wavelength_nm,
            pump_pulse_duration_fs,
            pump_energy_uj,
            refractive_index_optical,
            refractive_index_thz,
        }
    }

    /// ZnTe — the classic optical rectification material for 800 nm pumping.
    ///
    /// Parameters from literature:
    ///   d_eff = 68 pm/V, n_opt ≈ 2.85, n_THz ≈ 3.17 at 1 THz.
    pub fn znte() -> Self {
        Self::new(
            "ZnTe", 68.0,  // pm/V
            2.0,   // mm — reasonable for lab use
            800.0, // nm — Ti:Sapphire
            100.0, // fs
            1.0,   // μJ
            2.85, 3.17,
        )
    }

    /// LiNbO₃ with tilted-wavefront pumping scheme (yields higher efficiency).
    ///
    /// d_eff ≈ 168 pm/V, but phase matching requires pulse-front tilt.
    pub fn linbo3_thz() -> Self {
        Self::new(
            "LiNbO3", 168.0, 5.0,    // mm — longer crystal for LiNbO3
            1030.0, // nm — Yb laser pump
            300.0,  // fs
            100.0,  // μJ
            2.26, 4.96,
        )
    }

    /// Phase mismatch Δk = (ω_THz / c) · |n_opt − n_THz|  (rad m⁻¹).
    ///
    /// # Arguments
    /// * `freq_thz` — THz frequency (THz).
    pub fn phase_mismatch(&self, freq_thz: f64) -> f64 {
        let omega_thz = 2.0 * std::f64::consts::PI * freq_thz * 1e12; // rad/s
        omega_thz / C0 * (self.refractive_index_optical - self.refractive_index_thz).abs()
    }

    /// Coherence length L_c = π / |Δk|  (mm).
    ///
    /// Diverges (returned as `f64::MAX`) when the phase mismatch is zero
    /// (perfect phase matching).
    pub fn coherence_length_mm(&self, freq_thz: f64) -> f64 {
        let delta_k = self.phase_mismatch(freq_thz);
        if delta_k.abs() < 1e-30 {
            return f64::MAX;
        }
        let lc_m = std::f64::consts::PI / delta_k;
        lc_m * 1e3 // convert m → mm
    }

    /// Sinc phase-matching factor sinc²(Δk·L/2).
    fn sinc_factor(&self, freq_thz: f64) -> f64 {
        let delta_k = self.phase_mismatch(freq_thz);
        let l_m = self.crystal_length_mm * 1e-3;
        let arg = delta_k * l_m / 2.0;
        if arg.abs() < 1e-15 {
            1.0
        } else {
            (arg.sin() / arg).powi(2)
        }
    }

    /// THz electric field amplitude (V m⁻¹) — simplified plane-wave model:
    ///
    /// E_THz ∝ d_eff · L · E_pump² / (n_THz · c)
    ///
    /// The pump field amplitude is estimated from pulse energy, duration and
    /// beam area (assumed 1 mm² Gaussian spot).
    pub fn thz_field_amplitude_v_per_m(&self, freq_thz: f64) -> f64 {
        let d_eff_si = self.d_eff_pm_per_v * 1e-12; // pm/V → m/V
        let l_m = self.crystal_length_mm * 1e-3;

        // Pump intensity estimate: I = E / (τ · A)
        let tau_s = self.pump_pulse_duration_fs * 1e-15;
        let area_m2 = 1e-6; // 1 mm² beam spot
        let energy_j = self.pump_energy_uj * 1e-6;
        let intensity = energy_j / (tau_s * area_m2); // W/m²

        // Pump field: I = ½ c ε₀ n_opt E²
        let e_pump = (2.0 * intensity / (C0 * EPS0 * self.refractive_index_optical)).sqrt();

        // E_THz (simplified): ω_THz · d_eff · L · E_pump² / (n_THz · c)
        let omega_thz = 2.0 * std::f64::consts::PI * freq_thz * 1e12;
        let e_thz = omega_thz * d_eff_si * l_m * e_pump.powi(2) / (self.refractive_index_thz * C0);

        e_thz * self.sinc_factor(freq_thz).sqrt()
    }

    /// Estimated THz pulse energy (nJ).
    ///
    /// Rough expression: U_THz ≈ η · U_pump.
    pub fn thz_pulse_energy_nj(&self) -> f64 {
        let eta = self.conversion_efficiency();
        self.pump_energy_uj * eta * 1e3 // μJ → nJ via η
    }

    /// Optical-to-THz conversion efficiency (dimensionless, typically 10⁻³–10⁻²).
    ///
    /// Analytical estimate based on nonlinear interaction length and phase matching.
    pub fn conversion_efficiency(&self) -> f64 {
        let freq_thz = 1.0; // evaluate at 1 THz
        let d_eff_si = self.d_eff_pm_per_v * 1e-12;
        let l_m = self.crystal_length_mm * 1e-3;
        let omega_thz = 2.0 * std::f64::consts::PI * freq_thz * 1e12;

        // η ∝ (d_eff · ω · L)² / (n² · c² · ε₀)
        let tau_s = self.pump_pulse_duration_fs * 1e-15;
        let energy_j = self.pump_energy_uj * 1e-6;
        let area_m2 = 1e-6;
        let intensity = energy_j / (tau_s * area_m2);

        let numerator = (d_eff_si * omega_thz * l_m).powi(2) * intensity;
        let denominator =
            self.refractive_index_thz.powi(2) * self.refractive_index_optical * C0.powi(3) * EPS0;

        let eta_raw = numerator / denominator;
        // Phase-matching sinc² factor at 1 THz
        eta_raw * self.sinc_factor(freq_thz)
    }

    /// Optimal crystal length for maximum THz output at a given THz frequency.
    ///
    /// This equals the coherence length: L_opt = L_c = π / |Δk|  (mm).
    pub fn optimal_length_mm(&self, freq_thz: f64) -> f64 {
        self.coherence_length_mm(freq_thz)
    }

    /// THz bandwidth (THz) set by the pump pulse duration.
    ///
    /// For a transform-limited Gaussian pulse: Δν ≈ 0.44 / τ_FWHM.
    pub fn thz_bandwidth_thz(&self) -> f64 {
        let tau_ps = self.pump_pulse_duration_fs * 1e-3; // fs → ps
        0.44 / tau_ps
    }
}

// ─── Photoconductive Antenna ──────────────────────────────────────────────

/// Substrate material used in a photoconductive antenna.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum PcaSubstrate {
    /// Low-temperature-grown GaAs — classic emitter/detector (τ_c ≈ 0.3 ps).
    LowTemperatureGaAs,
    /// ErAs:GaAs — ultrafast variant.
    ErAsGaAs,
    /// InGaAs — suitable for 1550 nm excitation.
    InGaAs,
    /// Silicon — used for high-power applications.
    Si,
}

impl PcaSubstrate {
    /// Default carrier lifetime (ps) for each substrate.
    pub fn default_carrier_lifetime_ps(self) -> f64 {
        match self {
            Self::LowTemperatureGaAs => 0.3,
            Self::ErAsGaAs => 0.15,
            Self::InGaAs => 1.0,
            Self::Si => 10.0,
        }
    }

    /// Default electron mobility (cm² V⁻¹ s⁻¹).
    pub fn default_mobility_cm2_per_vs(self) -> f64 {
        match self {
            Self::LowTemperatureGaAs => 200.0,
            Self::ErAsGaAs => 100.0,
            Self::InGaAs => 3000.0,
            Self::Si => 1400.0,
        }
    }

    /// Refractive index of the substrate at THz frequencies.
    pub fn refractive_index_thz(self) -> f64 {
        match self {
            Self::LowTemperatureGaAs | Self::ErAsGaAs => 3.56,
            Self::InGaAs => 3.60,
            Self::Si => 3.42,
        }
    }
}

/// Photoconductive antenna (PCA) — THz emitter or detector driven by an
/// ultrashort optical pulse.
///
/// On the emitter side the optical pulse generates free carriers across the
/// biased gap; the fast transient photocurrent radiates THz.
#[derive(Debug, Clone)]
pub struct PhotoconductiveAntenna {
    /// Substrate material.
    pub substrate: PcaSubstrate,
    /// Antenna gap (μm).
    pub gap_um: f64,
    /// Full dipole length (μm).
    pub antenna_length_um: f64,
    /// Applied DC bias voltage (V).
    pub bias_voltage_v: f64,
    /// Photocarrier lifetime τ_c (ps).
    pub carrier_lifetime_ps: f64,
    /// Electron mobility (cm² V⁻¹ s⁻¹).
    pub mobility_cm2_per_vs: f64,
}

impl PhotoconductiveAntenna {
    /// General constructor.
    #[allow(clippy::too_many_arguments)]
    pub fn new(
        substrate: PcaSubstrate,
        gap_um: f64,
        antenna_length_um: f64,
        bias_voltage_v: f64,
        carrier_lifetime_ps: f64,
        mobility_cm2_per_vs: f64,
    ) -> Self {
        Self {
            substrate,
            gap_um,
            antenna_length_um,
            bias_voltage_v,
            carrier_lifetime_ps,
            mobility_cm2_per_vs,
        }
    }

    /// Standard LT-GaAs photoconductive antenna — widely used lab source.
    pub fn lt_gaas_standard() -> Self {
        let sub = PcaSubstrate::LowTemperatureGaAs;
        Self::new(
            sub,
            5.0,  // μm gap
            30.0, // μm dipole
            30.0, // V bias
            sub.default_carrier_lifetime_ps(),
            sub.default_mobility_cm2_per_vs(),
        )
    }

    /// Radiation impedance of the dipole antenna (Ω).
    ///
    /// For a Hertz dipole much shorter than λ the impedance is dominated by
    /// the geometric factor; we use the textbook value of 73 Ω for a
    /// half-wave resonant dipole as an upper reference.
    pub fn radiation_impedance_ohm(&self) -> f64 {
        73.0 // Ω — half-wave dipole
    }

    /// Photocurrent density J = e · μ · n_e · E_bias  (A m⁻²).
    ///
    /// # Arguments
    /// * `carrier_density` — instantaneous photo-carrier density (m⁻³).
    pub fn photocurrent_density_a_per_m2(&self, carrier_density: f64) -> f64 {
        let mu_si = self.mobility_cm2_per_vs * 1e-4; // cm²/Vs → m²/Vs
        let e_bias = self.bias_voltage_v / (self.gap_um * 1e-6); // V/m
        E_CHARGE * mu_si * carrier_density * e_bias
    }

    /// Simulate the THz E-field waveform emitted by the PCA.
    ///
    /// The photocurrent density is proportional to the convolution of the
    /// optical intensity envelope with an exponential carrier decay of
    /// time-constant τ_c.  The THz field is proportional to dJ/dt, here
    /// computed as a finite-difference derivative.
    ///
    /// # Arguments
    /// * `optical_pulse` — optical intensity envelope sampled at `dt_ps` intervals.
    /// * `dt_ps`         — time step (ps).
    ///
    /// Returns a `Vec<f64>` of THz E-field samples (relative units) with the
    /// same length as `optical_pulse`.
    pub fn thz_field_waveform(&self, optical_pulse: &[f64], dt_ps: f64) -> Vec<f64> {
        let n = optical_pulse.len();
        if n == 0 {
            return Vec::new();
        }

        // Bias field (V/m)
        let e_bias = self.bias_voltage_v / (self.gap_um * 1e-6);
        let mu_si = self.mobility_cm2_per_vs * 1e-4; // m²/Vs
        let tau_ps = self.carrier_lifetime_ps;
        let dt = dt_ps; // keep units in ps throughout

        // Convolve optical intensity with exponential carrier decay to get J(t)
        let mut j_current = vec![0.0_f64; n];
        for (i, j_val) in j_current.iter_mut().enumerate().take(n) {
            let mut acc = 0.0_f64;
            for (k, &pulse_k) in optical_pulse.iter().enumerate().take(i + 1) {
                let t_delay = (i - k) as f64 * dt;
                let decay = (-t_delay / tau_ps).exp();
                acc += pulse_k * decay * dt;
            }
            // Scale: J ∝ e · μ · E_bias · ∫ G(t') e^(-(t-t')/τ) dt'
            *j_val = E_CHARGE * mu_si * e_bias * acc;
        }

        // E_THz ∝ dJ/dt  (central finite differences, forward at boundaries)
        let mut e_thz = vec![0.0_f64; n];
        for i in 0..n {
            let dj = if i == 0 {
                j_current[1] - j_current[0]
            } else if i == n - 1 {
                j_current[n - 1] - j_current[n - 2]
            } else {
                (j_current[i + 1] - j_current[i - 1]) / 2.0
            };
            e_thz[i] = dj / (dt * 1e-12); // ps → s for correct SI units
        }
        e_thz
    }

    /// Estimate of radiated THz power (mW).
    ///
    /// Proportional to the square of the photocurrent, which scales linearly
    /// with optical pump power.
    pub fn radiated_power_mw(&self, optical_power_mw: f64) -> f64 {
        // Approximate responsivity: ~μA/mW → mA scale photocurrent
        let responsivity_a_per_w = 0.01; // A/W — typical PCA
        let i_photo = responsivity_a_per_w * optical_power_mw * 1e-3; // A
        let z_rad = self.radiation_impedance_ohm();
        let power_w = 0.5 * i_photo.powi(2) * z_rad;
        power_w * 1e3 // W → mW
    }

    /// Antenna resonance frequency (THz).
    ///
    /// f_res = c / (2 · L · n_sub)  where L is the dipole full length.
    pub fn resonance_frequency_thz(&self, n_substrate: f64) -> f64 {
        let l_m = self.antenna_length_um * 1e-6;
        C0 / (2.0 * l_m * n_substrate) * 1e-12 // Hz → THz
    }

    /// 3 dB bandwidth limited by the carrier lifetime (THz).
    ///
    /// BW ≈ 1 / (π · τ_c)
    pub fn bandwidth_thz(&self) -> f64 {
        let tau_s = self.carrier_lifetime_ps * 1e-12;
        1.0 / (std::f64::consts::PI * tau_s) * 1e-12 // Hz → THz
    }
}

// ─── Free-Electron Laser THz Source ──────────────────────────────────────

/// Free-electron laser (FEL) THz source — narrowband, high-power, coherent.
#[derive(Debug, Clone)]
pub struct FelThzSource {
    /// Centre frequency (THz).
    pub frequency_thz: f64,
    /// Average output power (W).
    pub power_w: f64,
    /// Linewidth FWHM (MHz).
    pub bandwidth_mhz: f64,
    /// Macro-pulse duration (ps).
    pub pulse_duration_ps: f64,
}

impl FelThzSource {
    /// General constructor.
    pub fn new(
        frequency_thz: f64,
        power_w: f64,
        bandwidth_mhz: f64,
        pulse_duration_ps: f64,
    ) -> Self {
        Self {
            frequency_thz,
            power_w,
            bandwidth_mhz,
            pulse_duration_ps,
        }
    }

    /// Free-space wavelength λ = c / f  (mm).
    pub fn wavelength_mm(&self) -> f64 {
        C0 / (self.frequency_thz * 1e12) * 1e3 // m → mm
    }

    /// Coherence length L_c = c / Δν  (mm).
    ///
    /// Δν is the spectral linewidth in Hz.
    pub fn coherence_length_mm(&self) -> f64 {
        let delta_nu_hz = self.bandwidth_mhz * 1e6;
        C0 / delta_nu_hz * 1e3 // m → mm
    }

    /// Returns `true` when the coherence time exceeds the pulse duration —
    /// i.e. the emission is temporally coherent.
    pub fn is_coherent(&self) -> bool {
        let coherence_time_ps = self.coherence_length_mm() / (C0 * 1e3) * 1e12;
        coherence_time_ps > self.pulse_duration_ps
    }
}

// ─── Tests ────────────────────────────────────────────────────────────────

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

    #[test]
    fn test_or_coherence_length_at_phase_match() {
        // Construct a crystal with identical optical and THz indices → perfect
        // phase matching → Δk = 0 → coherence length should diverge.
        let crystal =
            OpticalRectification::new("PerfectMatch", 68.0, 2.0, 800.0, 100.0, 1.0, 3.17, 3.17);
        let lc = crystal.coherence_length_mm(1.0);
        assert!(
            lc > 1e10,
            "coherence length should diverge when Δk = 0, got {lc}"
        );
    }

    #[test]
    fn test_or_bandwidth() {
        // 100 fs pulse → BW ≈ 0.44 / 0.1 ps = 4.4 THz
        let crystal = OpticalRectification::znte();
        let bw = crystal.thz_bandwidth_thz();
        let expected = 0.44 / (100e-3_f64); // ps units
        assert!(
            (bw - expected).abs() < 1e-6,
            "bandwidth = {bw}, expected {expected}"
        );
    }

    #[test]
    fn test_or_conversion_efficiency_small() {
        let crystal = OpticalRectification::znte();
        let eta = crystal.conversion_efficiency();
        assert!(eta < 0.01, "efficiency should be <1%, got {eta}");
        assert!(eta > 0.0, "efficiency must be positive, got {eta}");
    }

    #[test]
    fn test_pca_bandwidth() {
        // Standard LT-GaAs: τ_c = 0.3 ps → BW = 1/(π × 0.3e-12) Hz
        let pca = PhotoconductiveAntenna::lt_gaas_standard();
        let bw = pca.bandwidth_thz();
        let expected = 1.0 / (std::f64::consts::PI * 0.3e-12) * 1e-12;
        assert!(
            (bw - expected).abs() / expected < 1e-9,
            "PCA bandwidth = {bw} THz, expected ≈{expected} THz"
        );
    }

    #[test]
    fn test_pca_resonance_frequency_positive() {
        let pca = PhotoconductiveAntenna::lt_gaas_standard();
        let n_sub = PcaSubstrate::LowTemperatureGaAs.refractive_index_thz();
        let f_res = pca.resonance_frequency_thz(n_sub);
        assert!(
            f_res > 0.0,
            "resonance frequency must be positive, got {f_res}"
        );
    }

    #[test]
    fn test_pca_waveform_length() {
        let pca = PhotoconductiveAntenna::lt_gaas_standard();
        let pulse: Vec<f64> = (0..256)
            .map(|i| {
                let t = (i as f64 - 128.0) * 0.05;
                (-t * t).exp()
            })
            .collect();
        let waveform = pca.thz_field_waveform(&pulse, 0.05);
        assert_eq!(
            waveform.len(),
            pulse.len(),
            "waveform length must match input"
        );
    }

    #[test]
    fn test_fel_wavelength() {
        // 1 THz → λ = c/f = 0.3 mm
        let fel = FelThzSource::new(1.0, 100.0, 10.0, 10.0);
        let lambda_mm = fel.wavelength_mm();
        let expected_mm = C0 / 1e12 * 1e3;
        assert!(
            (lambda_mm - expected_mm).abs() < 1e-6,
            "λ = {lambda_mm} mm, expected {expected_mm} mm"
        );
    }

    // Suppress unused-import warning that the OxiPhotonError import would trigger
    // if no fallible methods were tested here.
    fn _assert_error_type(_: OxiPhotonError) {}
}