//! Pennes Bioheat Equation residual computation
//!
//! The Pennes equation models heat transfer in biological tissue:
//! ```text
//! ρc(∂T/∂t) = k∇²T + ωb·ρb·cb·(Ta - T) + Qm + Qs
//! ```
//!
//! Rearranged as residual (should equal zero):
//! ```text
//! R = ρc·∂T/∂t - k·∇²T - ωb·ρb·cb·(Ta - T) - Qm - Qs = 0
//! ```

use bioheat_shared::BioheatParams;
use serde::{Deserialize, Serialize};

/// Pennes bioheat equation residual computer
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct PennesResidual {
    /// Tissue and blood properties
    pub params: BioheatParams,
}

/// Derivatives of temperature field needed for residual computation
#[derive(Debug, Clone, Copy)]
pub struct TemperatureDerivatives {
    /// Temperature value T
    pub t: f32,
    /// Time derivative ∂T/∂t
    pub dt_dt: f32,
    /// First spatial derivatives
    pub dt_dx: f32,
    pub dt_dy: f32,
    pub dt_dz: f32,
    /// Second spatial derivatives (for Laplacian)
    pub d2t_dx2: f32,
    pub d2t_dy2: f32,
    pub d2t_dz2: f32,
}

impl TemperatureDerivatives {
    /// Compute the Laplacian ∇²T = ∂²T/∂x² + ∂²T/∂y² + ∂²T/∂z²
    #[must_use]
    pub fn laplacian(&self) -> f32 {
        self.d2t_dx2 + self.d2t_dy2 + self.d2t_dz2
    }

    /// Create derivatives for a constant temperature field
    #[must_use]
    pub fn constant(t: f32) -> Self {
        Self {
            t,
            dt_dt: 0.0,
            dt_dx: 0.0,
            dt_dy: 0.0,
            dt_dz: 0.0,
            d2t_dx2: 0.0,
            d2t_dy2: 0.0,
            d2t_dz2: 0.0,
        }
    }
}

impl PennesResidual {
    /// Create a new Pennes residual computer with given parameters
    #[must_use]
    pub fn new(params: BioheatParams) -> Self {
        Self { params }
    }

    /// Compute the Pennes equation residual at a single point
    ///
    /// # Arguments
    /// * `derivs` - Temperature and its derivatives at the point
    /// * `heat_source` - External heat source Qs (W/m³) from probe
    ///
    /// # Returns
    /// The residual value (should be zero for a correct solution)
    #[must_use]
    pub fn compute(&self, derivs: &TemperatureDerivatives, heat_source: f32) -> f32 {
        let tissue = &self.params.tissue;
        let blood = &self.params.blood;

        // Volumetric heat capacity ρc (J/m³/K)
        let rho_c = tissue.density * tissue.specific_heat;

        // Thermal conductivity k (W/m/K)
        let k = tissue.thermal_conductivity;

        // Perfusion coefficient ωb·ρb·cb (W/m³/K)
        let perfusion_coeff = self.params.perfusion_coefficient();

        // Metabolic heat Qm (W/m³)
        let q_m = tissue.metabolic_heat;

        // Arterial temperature Ta (°C)
        let t_a = blood.arterial_temperature;

        // Compute each term:
        // 1. Time derivative term: ρc·∂T/∂t
        let time_term = rho_c * derivs.dt_dt;

        // 2. Diffusion term: k·∇²T
        let diffusion_term = k * derivs.laplacian();

        // 3. Perfusion term: ωb·ρb·cb·(Ta - T)
        // Note: This is a cooling term when T > Ta
        let perfusion_term = perfusion_coeff * (t_a - derivs.t);

        // 4. Total source: Qm + Qs
        let source_term = q_m + heat_source;

        // Residual: ρc·∂T/∂t - k·∇²T - perfusion - Qm - Qs = 0
        time_term - diffusion_term - perfusion_term - source_term
    }

    /// Compute residual for steady-state (∂T/∂t = 0)
    #[must_use]
    pub fn compute_steady_state(&self, derivs: &TemperatureDerivatives, heat_source: f32) -> f32 {
        let mut steady_derivs = *derivs;
        steady_derivs.dt_dt = 0.0;
        self.compute(&steady_derivs, heat_source)
    }

    /// Compute residual for batch of points
    pub fn compute_batch(
        &self,
        derivs_batch: &[TemperatureDerivatives],
        heat_sources: &[f32],
    ) -> Vec<f32> {
        derivs_batch
            .iter()
            .zip(heat_sources)
            .map(|(d, &q)| self.compute(d, q))
            .collect()
    }

    /// Get the characteristic time scale (thermal diffusion time)
    /// τ = L² / α where α = k/(ρc) is thermal diffusivity
    #[must_use]
    pub fn thermal_time_scale(&self, length_scale: f32) -> f32 {
        let alpha = self.params.tissue.thermal_diffusivity();
        length_scale * length_scale / alpha
    }

    /// Get the perfusion time scale
    /// τ_p = (ρc) / (ωb·ρb·cb)
    #[must_use]
    pub fn perfusion_time_scale(&self) -> f32 {
        let rho_c = self.params.tissue.volumetric_heat_capacity();
        let perf_coeff = self.params.perfusion_coefficient();
        if perf_coeff > 1e-10 {
            rho_c / perf_coeff
        } else {
            f32::INFINITY // No perfusion
        }
    }
}

/// Analytical test case: 1D heat conduction with constant source
/// Useful for validating the residual computation
pub mod analytical {

    /// 1D steady-state solution with constant heat source
    /// For boundary conditions T(0) = T(L) = T_boundary
    /// and constant source Q, the solution is:
    /// T(x) = T_boundary + (Q/2k) * x * (L - x)
    pub struct SteadyStateBar {
        pub length: f32,
        pub t_boundary: f32,
        pub heat_source: f32,
        pub k: f32,
    }

    impl SteadyStateBar {
        /// Temperature at position x
        #[must_use]
        pub fn temperature(&self, x: f32) -> f32 {
            let q_over_2k = self.heat_source / (2.0 * self.k);
            self.t_boundary + q_over_2k * x * (self.length - x)
        }

        /// First derivative dT/dx
        #[must_use]
        pub fn dt_dx(&self, x: f32) -> f32 {
            let q_over_2k = self.heat_source / (2.0 * self.k);
            q_over_2k * (self.length - 2.0 * x)
        }

        /// Second derivative d²T/dx²
        #[must_use]
        pub fn d2t_dx2(&self) -> f32 {
            -self.heat_source / self.k
        }
    }
}

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

    #[test]
    fn test_pennes_residual_equilibrium() {
        // At equilibrium with no heat source and T = Ta, residual should be ~0
        let params = BioheatParams::liver_ablation();
        let residual = PennesResidual::new(params.clone());

        let derivs = TemperatureDerivatives {
            t: params.blood.arterial_temperature, // T = Ta
            dt_dt: 0.0,                           // Steady state
            dt_dx: 0.0,
            dt_dy: 0.0,
            dt_dz: 0.0,
            d2t_dx2: 0.0,
            d2t_dy2: 0.0,
            d2t_dz2: 0.0,
        };

        // Only metabolic heat contributes
        let r = residual.compute(&derivs, 0.0);
        assert!((r + params.tissue.metabolic_heat).abs() < 1.0);
    }

    #[test]
    fn test_laplacian() {
        let derivs = TemperatureDerivatives {
            t: 50.0,
            dt_dt: 0.0,
            dt_dx: 1.0,
            dt_dy: 2.0,
            dt_dz: 3.0,
            d2t_dx2: 1.0,
            d2t_dy2: 2.0,
            d2t_dz2: 3.0,
        };
        assert!((derivs.laplacian() - 6.0).abs() < 1e-6);
    }

    #[test]
    fn test_analytical_bar() {
        use analytical::SteadyStateBar;

        let bar = SteadyStateBar {
            length: 1.0,
            t_boundary: 37.0,
            heat_source: 1000.0,
            k: 0.5,
        };

        // At x=0 and x=L, T should equal boundary
        assert!((bar.temperature(0.0) - 37.0).abs() < 1e-6);
        assert!((bar.temperature(1.0) - 37.0).abs() < 1e-6);

        // Maximum at center
        let t_center = bar.temperature(0.5);
        assert!(t_center > 37.0);

        // Derivative at center should be 0
        assert!(bar.dt_dx(0.5).abs() < 1e-6);
    }

    #[test]
    fn test_time_scales() {
        let params = BioheatParams::liver_ablation();
        let residual = PennesResidual::new(params);

        // Thermal time scale for 1cm domain
        let tau_thermal = residual.thermal_time_scale(0.01);
        assert!(tau_thermal > 0.0);

        // Perfusion time scale
        let tau_perf = residual.perfusion_time_scale();
        assert!(tau_perf > 0.0);
        assert!(tau_perf < f32::INFINITY);
    }
}