// gate.sy -- H-H Gating Dynamics
// Hodgkin-Huxley gates, FitzHugh-Nagumo glial modulation.
// No sigmoid anywhere. Alpha/beta rate constants only.
//
// Aislinn McCubbins, April 2026

// === H-H GATE ===
// dm/dt = (m_inf - m) / tau
// m_inf = alpha / (alpha + beta)  -- NOT sigmoid
// tau = 1 / (alpha + beta)

temporal struct HHGate {
    state m: Float = 0.05,
    state h: Float = 0.9,
    const use_inactivation: Bool = true,
}

impl Evolve for HHGate {
    fn evolve(&mut self, dt: Duration) {
        // Subclasses override with specific alpha/beta
    }
}

impl HHGate {
    fn gate_value(&self) -> Float {
        if self.use_inactivation {
            self.m * self.h
        } else {
            self.m
        }
    }

    fn reset(&mut self) {
        self.m = 0.05
        self.h = 0.9
    }
}

// === FITZHHUGH-NAGUMO GLIAL GATE ===
// Emergent thresholds from bifurcation dynamics.
// No sigmoid. The fold IS the threshold.

temporal struct GlialGate {
    state v_mu: Float = -1.2,
    state w_mu: Float = -0.6,
    state calcium: Float = 0.0,
    state energy: Float = 1.0,

    const epsilon: Float = 0.05,
    const a: Float = 0.7,
    const b: Float = 0.8,
    const tau_calcium: Float = 100.0,
    const depletion_rate: Float = 0.02,
    const recovery_rate: Float = 0.005,
}

impl Evolve for GlialGate {
    fn evolve(&mut self, dt: Duration) {
        let dt_s = dt.seconds()

        // Calcium decay
        self.calcium = self.calcium * 0.99

        // FHN dynamics -- cubic nullcline, emergent threshold
        let dv = (self.v_mu - self.v_mu * self.v_mu * self.v_mu / 3.0 - self.w_mu + self.calcium * 0.5) * dt_s
        let dw = self.epsilon * (self.v_mu + self.a - self.b * self.w_mu) * dt_s

        self.v_mu = self.v_mu + dv
        self.w_mu = self.w_mu + dw

        // Energy recovery
        self.energy = self.energy + self.recovery_rate * (1.0 - self.energy) * dt_s
    }
}

impl GlialGate {
    fn record_activity(&mut self, intensity: Float) {
        self.calcium = self.calcium + intensity
        self.energy = self.energy - self.depletion_rate * intensity
        if self.energy < 0.1 {
            self.energy = 0.1
        }
    }

    fn route_availability(&self) -> Float {
        // V > 0 means past upper fold -> route suppressed
        // Linear mapping, no sigmoid
        let gate = 0.5 - 0.25 * self.v_mu
        if gate < 0.0 { return 0.0 }
        if gate > 1.0 { return 1.0 }
        gate
    }

    fn modulation(&self) -> Float {
        // How much to slow gate time constants
        // 1.0 = normal, up to 10.0 = 10x slower
        let mu_max = 10.0
        let max_slowdown = 10.0
        let mu_clamped = self.v_mu
        if mu_clamped < 0.0 { mu_clamped = 0.0 }
        if mu_clamped > mu_max { mu_clamped = mu_max }
        1.0 + (max_slowdown - 1.0) * (mu_clamped / mu_max)
    }
}
