//! Ionic models for cardiac action potential simulation.
//!
//! Implements various ionic models from simple to detailed:
//! - Mitchell-Schaeffer (2 variables)
//! - FitzHugh-Nagumo (2 variables)
//! - Aliev-Panfilov (2 variables)
//! - ten Tusscher-Panfilov (17 variables)

use cardiosim_shared::IonicModel;

/// Ionic state for a cell.
#[derive(Debug, Clone)]
pub struct IonicState {
    model: IonicModel,
    /// Recovery variable (w or h)
    pub recovery: f32,
    /// Additional state variables for complex models
    pub aux_vars: Vec<f32>,
}

impl IonicState {
    /// Create a new ionic state.
    #[must_use]
    pub fn new(model: IonicModel) -> Self {
        let aux_vars = match model {
            IonicModel::TenTusscherPanfilov | IonicModel::OHaraRudy => {
                // Many gating variables
                vec![0.0; 15]
            }
            _ => vec![],
        };

        Self {
            model,
            recovery: 0.0,
            aux_vars,
        }
    }

    /// Compute ionic currents.
    #[must_use]
    pub fn compute_currents(&self, voltage: f32) -> (f32, f32) {
        match self.model {
            IonicModel::MitchellSchaeffer => self.mitchell_schaeffer(voltage),
            IonicModel::FitzHughNagumo => self.fitzhugh_nagumo(voltage),
            IonicModel::AlievPanfilov => self.aliev_panfilov(voltage),
            IonicModel::TenTusscherPanfilov => self.ten_tusscher(voltage),
            IonicModel::OHaraRudy => self.ohara_rudy(voltage),
        }
    }

    /// Mitchell-Schaeffer model.
    /// Simple 2-variable model with good action potential shape.
    fn mitchell_schaeffer(&self, v: f32) -> (f32, f32) {
        // Parameters
        let tau_in = 0.3;
        let tau_out = 6.0;
        let tau_open = 120.0;
        let tau_close = 150.0;
        let v_gate = 0.13;

        // Normalize voltage to [0, 1]
        let u = (v + 85.0) / 120.0;
        let u = u.clamp(0.0, 1.0);

        let h = self.recovery;

        // Currents
        let j_in = h * u * u * (1.0 - u) / tau_in;
        let j_out = -u / tau_out;
        let j_ion = j_in + j_out;

        // Gate dynamics
        let dh = if u < v_gate {
            (1.0 - h) / tau_open
        } else {
            -h / tau_close
        };

        // Scale back to mV
        (j_ion * 120.0, dh)
    }

    /// FitzHugh-Nagumo model.
    /// Classic excitable media model.
    fn fitzhugh_nagumo(&self, v: f32) -> (f32, f32) {
        // Parameters
        let a = 0.7;
        let b = 0.8;
        let tau = 12.5;
        let epsilon = 0.08;

        // Normalize
        let u = (v + 85.0) / 120.0;
        let u = u.clamp(-0.5, 1.5);
        let w = self.recovery;

        // Cubic nullcline
        let du = u - u * u * u / 3.0 - w;
        let dw = epsilon * (u + a - b * w);

        (du * 120.0 / tau, dw)
    }

    /// Aliev-Panfilov model.
    /// Modified FHN with more realistic restitution.
    fn aliev_panfilov(&self, v: f32) -> (f32, f32) {
        // Parameters
        let k = 8.0;
        let a = 0.15;
        let epsilon0 = 0.002;
        let mu1 = 0.2;
        let mu2 = 0.3;

        // Normalize
        let u = (v + 85.0) / 120.0;
        let u = u.clamp(0.0, 1.0);
        let w = self.recovery;

        // Dynamics
        let du = -k * u * (u - a) * (u - 1.0) - u * w;

        let epsilon = epsilon0 + mu1 * w / (u + mu2);
        let dw = epsilon * (-w - k * u * (u - a - 1.0));

        (du * 120.0, dw)
    }

    /// ten Tusscher-Panfilov model (simplified).
    /// More detailed model with major ionic currents.
    fn ten_tusscher(&self, v: f32) -> (f32, f32) {
        // Simplified version with main currents

        // Reversal potentials
        let e_na = 70.0;
        let e_k = -88.0;
        let e_ca = 120.0;

        // Maximum conductances (mS/cm²)
        let g_na = 14.838;
        let g_k1 = 5.405;
        let g_to = 0.294;
        let g_cal = 0.000175;

        // Gating (simplified)
        let m_inf = 1.0 / (1.0 + ((-v - 40.0) / 9.0).exp());
        let h_inf = 1.0 / (1.0 + ((v + 70.0) / 7.0).exp());
        let j_inf = h_inf;

        let d_inf = 1.0 / (1.0 + ((-v - 5.0) / 6.0).exp());
        let f_inf = 1.0 / (1.0 + ((v + 20.0) / 7.0).exp());

        let xr1_inf = 1.0 / (1.0 + ((-v - 26.0) / 7.0).exp());

        // Use recovery as simplified gating
        let h = 0.5 + 0.5 * self.recovery;

        // Currents
        let i_na = g_na * m_inf.powi(3) * h * j_inf * (v - e_na);
        let i_cal = g_cal * d_inf * f_inf * (v - e_ca);
        let i_k1 = g_k1 * (v - e_k) / (1.0 + (0.1 * (v - e_k)).exp());
        let i_to = g_to * xr1_inf * self.recovery * (v - e_k);

        let i_ion = i_na + i_cal + i_k1 + i_to;

        // Recovery dynamics (simplified)
        let tau_h = 5.0 + 100.0 / (1.0 + ((v + 40.0) / 10.0).exp());
        let dh = (h_inf - self.recovery) / tau_h;

        (-i_ion, dh)
    }

    /// O'Hara-Rudy model (simplified).
    fn ohara_rudy(&self, v: f32) -> (f32, f32) {
        // Very simplified version
        // Full model has 41 state variables

        let e_na = 70.0;
        let e_k = -88.0;

        let g_na_fast = 75.0;
        let g_na_late = 0.0075;

        // Fast sodium
        let m_inf = 1.0 / (1.0 + (-(v + 39.57) / 9.871).exp());
        let h_inf = 1.0 / (1.0 + ((v + 82.9) / 6.086).exp());

        let h = 0.5 + 0.5 * self.recovery;

        let i_na_fast = g_na_fast * m_inf.powi(3) * h.powi(2) * (v - e_na);
        let i_na_late = g_na_late * m_inf.powi(3) * (1.0 - h) * (v - e_na);

        // Potassium
        let xk1_inf = 1.0 / (1.0 + ((v + 2.55 + e_k) / 34.2).exp());
        let i_k1 = 0.1 * xk1_inf * (v - e_k);

        let i_ion = i_na_fast + i_na_late + i_k1;

        // Recovery
        let tau_h = 2.0 + 50.0 / (1.0 + ((v + 50.0) / 10.0).exp());
        let dh = (h_inf - self.recovery) / tau_h;

        (-i_ion, dh)
    }
}

/// Action potential metrics.
#[derive(Debug, Clone)]
pub struct APMetrics {
    /// Peak voltage (mV)
    pub v_max: f32,
    /// Resting voltage (mV)
    pub v_rest: f32,
    /// Maximum upstroke velocity (mV/ms)
    pub dv_dt_max: f32,
    /// APD at 50% repolarization (ms)
    pub apd50: f32,
    /// APD at 90% repolarization (ms)
    pub apd90: f32,
}

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

    #[test]
    fn test_ionic_state_creation() {
        let state = IonicState::new(IonicModel::MitchellSchaeffer);
        assert_eq!(state.recovery, 0.0);
    }

    #[test]
    fn test_mitchell_schaeffer() {
        let state = IonicState::new(IonicModel::MitchellSchaeffer);
        let (dv, dw) = state.compute_currents(-85.0);
        // At rest, should have small currents
        assert!(dv.abs() < 10.0);
        assert!(dw.abs() < 1.0);
    }

    #[test]
    fn test_mitchell_schaeffer_upstroke() {
        let mut state = IonicState::new(IonicModel::MitchellSchaeffer);
        state.recovery = 0.8; // Gate open
        let (dv, _) = state.compute_currents(-40.0);
        // Should have positive dv (upstroke)
        assert!(dv > 0.0);
    }

    #[test]
    fn test_fitzhugh_nagumo() {
        let state = IonicState::new(IonicModel::FitzHughNagumo);
        let (dv, dw) = state.compute_currents(-85.0);
        assert!(dv.is_finite());
        assert!(dw.is_finite());
    }

    #[test]
    fn test_aliev_panfilov() {
        let state = IonicState::new(IonicModel::AlievPanfilov);
        let (dv, dw) = state.compute_currents(-85.0);
        assert!(dv.is_finite());
        assert!(dw.is_finite());
    }

    #[test]
    fn test_ten_tusscher() {
        let state = IonicState::new(IonicModel::TenTusscherPanfilov);
        let (dv, dw) = state.compute_currents(-85.0);
        assert!(dv.is_finite());
        assert!(dw.is_finite());
    }

    #[test]
    fn test_ohara_rudy() {
        let state = IonicState::new(IonicModel::OHaraRudy);
        let (dv, dw) = state.compute_currents(-85.0);
        assert!(dv.is_finite());
        assert!(dw.is_finite());
    }

    #[test]
    fn test_action_potential_cycle() {
        // Simulate one action potential
        let mut state = IonicState::new(IonicModel::MitchellSchaeffer);
        let mut v = -85.0;
        let dt = 0.1;

        // Apply stimulus
        v = -40.0;
        state.recovery = 0.8;

        let mut v_max = v;
        for _ in 0..1000 {
            let (dv, dw) = state.compute_currents(v);
            v += dt * dv;
            state.recovery += dt * dw;

            if v > v_max {
                v_max = v;
            }
        }

        // Should have depolarized
        assert!(v_max > 0.0);
    }
}