"""Проверка ОДУ разряда против аналитического решения простого RLC-контура.

dL/dx заморожена в ноль (mu_eff=1.0), поэтому механика отвязана от контура
и [Q, I] эволюционируют как классический последовательный RLC-разряд
конденсатора. Сравниваем численное решение solve_ivp с учебным аналитическим
решением во всех трёх режимах затухания.
"""

import math

import numpy as np
import pytest
from scipy.integrate import solve_ivp

from gausse.physics.circuit import Q, StageCircuitParams, derivatives
from gausse.physics.inductance import CoilInductanceModel

L_H = 100e-6
C_F = 200e-6
V0 = 400.0
Q0 = C_F * V0


def _make_frozen_model() -> CoilInductanceModel:
    return CoilInductanceModel(
        l_air_h=L_H, coil_length_m=0.05, slug_length_m=0.03, mu_eff=1.0, smoothing_width_m=0.001
    )


def _integrate(r_ohm: float, t_span: tuple[float, float], t_eval: np.ndarray):
    model = _make_frozen_model()
    params = StageCircuitParams(
        capacitance_f=C_F,
        r_total_ohm=r_ohm,
        mass_kg=1.0,
        mu_eff=1.0,
        total_turns=100,
        coil_length_m=0.05,
        b_sat_tesla=1e9,  # эффективно отключает клэмп насыщения для этого теста
    )
    state0 = [Q0, 0.0, 0.0, 0.0]
    return solve_ivp(
        derivatives,
        t_span,
        state0,
        args=(model, params),
        t_eval=t_eval,
        method="RK45",
        rtol=1e-10,
        atol=1e-14,
    )


def test_underdamped_matches_analytical_solution():
    r_ohm = 0.05
    alpha = r_ohm / (2 * L_H)
    omega0 = 1 / math.sqrt(L_H * C_F)
    omega_d = math.sqrt(omega0**2 - alpha**2)

    t_eval = np.linspace(0, 5 / alpha, 200)
    sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)

    expected_q = Q0 * np.exp(-alpha * t_eval) * (
        np.cos(omega_d * t_eval) + (alpha / omega_d) * np.sin(omega_d * t_eval)
    )
    expected_i = Q0 * (omega0**2 / omega_d) * np.exp(-alpha * t_eval) * np.sin(omega_d * t_eval)

    assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4)
    assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=abs(expected_i).max() * 1e-3)


def test_overdamped_matches_analytical_solution():
    r_ohm = 5.0
    alpha = r_ohm / (2 * L_H)
    omega0 = 1 / math.sqrt(L_H * C_F)
    root = math.sqrt(alpha**2 - omega0**2)
    s1, s2 = -alpha + root, -alpha - root

    t_eval = np.linspace(0, 5 / abs(s1), 200)
    sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)

    a_coef = -Q0 * s2 / (s1 - s2)
    b_coef = Q0 * s1 / (s1 - s2)
    expected_q = a_coef * np.exp(s1 * t_eval) + b_coef * np.exp(s2 * t_eval)
    expected_i = -(a_coef * s1 * np.exp(s1 * t_eval) + b_coef * s2 * np.exp(s2 * t_eval))

    assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4)
    assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=max(abs(expected_i).max(), 1e-9) * 1e-3)


def test_critically_damped_matches_analytical_solution():
    alpha = 1 / math.sqrt(L_H * C_F)  # R = 2*sqrt(L/C) => alpha = omega0
    r_ohm = 2 * alpha * L_H

    t_eval = np.linspace(0, 5 / alpha, 200)
    sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)

    expected_q = Q0 * (1 + alpha * t_eval) * np.exp(-alpha * t_eval)
    expected_i = alpha**2 * Q0 * t_eval * np.exp(-alpha * t_eval)

    assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4)
    assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=max(abs(expected_i).max(), 1e-9) * 1e-3)


def test_energy_conserved_when_lossless():
    """С R=0 энергия конденсатора должна переходить в энергию катушки без потерь."""
    r_ohm = 0.0
    t_eval = np.linspace(0, 1e-3, 500)
    sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)

    cap_energy = sol.y[Q] ** 2 / (2 * C_F)
    inductor_energy = 0.5 * L_H * sol.y[1] ** 2
    total_energy = cap_energy + inductor_energy
    initial_energy = Q0**2 / (2 * C_F)

    assert np.allclose(total_energy, initial_energy, rtol=1e-4)