r"""
Does the radial expansion change the discharge? 7 kV into CH4 through a 6.2 m cable
===================================================================================

A nanosecond pulse heats a thin channel at nearly constant volume, so the core
ends up over-pressured and pushes outward. The question this example answers is
whether that motion feeds back on the *electrical* behaviour, or whether it is
something that merely happens afterwards.

It is answered by running the **same** discharge through the same circuit
-- 7 kV trapezoid, 1 ohm source, 6.2 m of 75 ohm cable with reflections, 3.8 mm
of CH4 at one atmosphere, 300 ns -- under four models: the 1-D radial channel
with :math:`\rho` and :math:`u_r` solved (``hydrodynamics=True``), the same
channel with the density frozen, and a 0-D lumped reactor at constant volume
and on a polytropic expansion. Everything else is identical, so each difference
belongs to one modelling choice.

**Within the 1-D channel, the answer has two halves, and only the first is the
obvious one:**

* **The first current pulse does not care.** Peak current is 184.089 A with and
  without the sub-step, and the resistance minimum is 0.04937 vs 0.04942 ohm, a
  part in 2000. The channel breaks down in about 12 ns; the gas needs hundreds
  of nanoseconds to move a kernel radius. The expansion has not happened yet
  while the current is being drawn.
* **The cable echoes do care.** The 6.2 m line returns a reflection every
  :math:`2L/c = 65.3` ns, so the third arrives at 196 ns, by which time the
  channel has rarefied substantially on axis. The current it drives peaks a few
  nanoseconds later at 118.9 A into the expanded channel against 125.9 A into
  the frozen one, and over the 300 ns window the expanded channel takes
  3850 uJ against 4090 uJ, **5.9 % less**.

So the expansion is real -- 1461 m/s, 20.6 atm on axis at peak, 56.5 % of the
axis density gone by 300 ns -- and it is invisible to the pulse that caused it.
What it changes is what the *line* sees afterwards. Without a long cable there
would be nothing for it to change.

**Against the 0-D reactor, two more things separate:**

* **Geometry dominates the resistance.** The lumped cell holds one conductivity
  over a fixed 0.60 mm disc and bottoms out at 6.16 ohm; the 1-D channel
  integrates :math:`\int \sigma \, 2 \pi r \, \mathrm{d}r` over a core that
  ionizes far harder than the mean, and reaches 0.049 ohm. That is a factor of
  125, and it is the radial resolution, not the flow.
* **The current barely notices, because the line sets it.** Peak current is
  170.2 A in 0-D against 184.1 A in 1-D, 8 % apart despite those two decades of
  resistance: once the channel conducts, 75 ohms of cable is the larger
  impedance and the plasma stops being what limits the current.
* **A lumped expansion is not a small version of a radial one.** The polytropic
  cell grows by 15 % in volume over the window, a 13 % density drop, against
  56.5 % on the 1-D axis. A volume law driven by the cell's *mean* pressure
  against ambient cannot represent a rarefaction running inward from the core,
  which is the mechanism that empties the axis.

Conduction and species diffusion are off (``kappa = 0``), so what is plotted is
the pulse, the chemistry and the hydrodynamics alone. The 0-D radius is a
choice rather than a measurement, and the resistance comparison moves with it;
the derivation used is stated where it is set.

.. tags:: plasma, Cantera, two-temperature, transient, channel, hydrodynamics,
          circuit, NRP, 1D
"""  # noqa: D205

# %%
# Import the required libraries.
# ------------------------------

import pathlib
import tempfile
import time

import cantera as ct
import matplotlib.pyplot as plt
import numpy as np

import rizer.kinetics.extensible_rate  # noqa: F401 (registers the CH4 rates)
import rizer.misc.units as u
from rizer.electrical_model.circuit.transmission_line_circuit import (
    TransmissionLineResistiveCircuit,
)
from rizer.electrical_model.components.cable import IdealCable
from rizer.electrical_model.components.generator import TrapezoidalGenerator
from rizer.electrical_model.components.source_circuit import ResistiveSourceCircuit
from rizer.misc.ct_utils import partially_ionized_mole_fractions
from rizer.misc.plt_utils import set_mpl_style
from rizer.misc.utils import get_path_to_data
from rizer.models.nrp.isomass_2T_volume_reactor_cpp import Isomass2TVolumeReactor
from rizer.models.nrp.plasma_channel_split_2T1D import ChannelReactor1D
from rizer.transport.loaders import get_default_momentum_transfer_collision_frequencies

set_mpl_style()

# %%
# Discharge, circuit and grid.
# ----------------------------
#
# The applied field is not ``U_on / gap``: the circuit sets the plasma voltage
# against the channel's own resistance
# :math:`R_p = \mathrm{gap} / \int \sigma \,\mathrm{d}A`, so a more conductive
# channel loads the line and lowers its own field. The cable length is what
# makes this example about more than the first 12 ns: at ``c = 1.9e8`` m/s a
# 6.2 m line brings each reflection back after ``2 L / c``.

mechanism = str(get_path_to_data("mechanisms") / "Goutier2025" / "CH4_to_C2H2.yaml")

gap = 3.8e-3  # Inter-electrode gap [m]
kernel_radius = 500e-6  # Seed ionization width (FWHM) [m]
r_wall = 4.0e-3  # Domain outer radius [m]

p_0 = ct.one_atm  # Initial pressure [Pa]
Tg_0 = 1000.0  # Initial gas and electron temperature [K]
ne_0 = 1.0e19  # Peak seed electron density [m^-3]

U_on = 7000.0  # Generator voltage when on [V]
t_rise, t_on, t_fall = 3e-9, 6e-9, 3e-9  # Trapezoid timings [s]
cable_length, cable_speed = 6.2, 1.9e8  # [m], [m/s]

t_end = 300e-9  # End time [s]
dt_out = 1e-9  # Output cadence [s]
n_radial_points = 31  # Radial nodes [-]

# The 0D reactor lumps the channel into one cell and so needs a single radius,
# which a Gaussian seed does not uniquely supply. The top-hat radius carrying
# the same electron inventory at the peak density is used: for a Gaussian of
# standard deviation `sigma`, int(n_e * 2 pi r dr) = n_e0 * 2 pi sigma^2, which
# a disc of radius sigma*sqrt(2) at n_e0 matches. Any other defensible choice
# moves the 0D resistance, which is itself part of what the comparison shows.
seed_sigma = kernel_radius / np.sqrt(2.0 * np.log(2.0))  # [m]
radius_0d = seed_sigma * np.sqrt(2.0)  # [m]
#: Heat-capacity ratio of the seeded CH4 at the initial state [-], used as the
#: 0D polytropic index for its isentropic case. The native reactor takes a
#: number rather than the symbolic "gamma", so it is held fixed for the run
#: while the real ratio falls to about 1.07 by 5000 K.
polytropic_index_0d = 1.1272

round_trip = 2.0 * cable_length / cable_speed  # Cable echo period [s]
print(f"cable round trip 2L/c = {round_trip * 1e9:.1f} ns")


def make_circuit():
    """Build the NRP circuit (trapezoidal generator + ideal cable).

    Returns
    -------
    TransmissionLineResistiveCircuit
        The 7 kV source behind 1 ohm, feeding 6.2 m of 75 ohm line with
        reflections enabled.
    """
    return TransmissionLineResistiveCircuit(
        source=ResistiveSourceCircuit(
            TrapezoidalGenerator(
                U_on=U_on, U_off=0.0, t_rise=t_rise, t_on=t_on, t_fall=t_fall
            ),
            R_g=1.0,
        ),
        cable=IdealCable(L=cable_length, Z_c=75.0, c=cable_speed),
        include_reflections=True,
    )


# %%
# One seed, two runs.
# -------------------
#
# Both channels start from the identical Gaussian ionized core on the identical
# stretched grid, so nothing but ``hydrodynamics`` separates them. ``Y_amb`` is
# the ambient composition the outer ghost cell is held at; it is required once
# the sub-step is on, because the species fluxes only telescope into the mass
# flux at the wall when it sums to one.

plasma = ct.Solution(mechanism, "plasma", transport_model=None)
names = list(plasma.species_names)
mtcf = get_default_momentum_transfer_collision_frequencies(names)

r = r_wall * np.sinh(3.5 * np.linspace(0.0, 1.0, n_radial_points)) / np.sinh(3.5)
width = kernel_radius / np.sqrt(2.0 * np.log(2.0))
ne_profile = ne_0 * np.exp(-(r**2) / (2.0 * width**2)) + 1.0e12

gas = ct.Solution(mechanism, "plasma")
Y0 = np.zeros((n_radial_points, plasma.n_species))
for node in range(n_radial_points):
    x_e = ne_profile[node] * u.k_b * Tg_0 / p_0
    gas.TPX = Tg_0, p_0, partially_ionized_mole_fractions(x_e)
    Y0[node] = gas.Y
gas.TPX = Tg_0, p_0, "CH4:1"
Y_ambient = gas.Y.copy()

x_e = ne_0 * u.k_b * Tg_0 / p_0
plasma.TPX = Tg_0, p_0, partially_ionized_mole_fractions(x_e)
rho_0 = float(plasma.density)
t_out = np.arange(0.0, t_end + 1e-12, dt_out)


def solve(hydrodynamics):
    """Run the circuit-coupled channel with the radial sub-step on or off.

    Parameters
    ----------
    hydrodynamics : bool
        Whether to solve continuity, radial momentum, species advection and the
        two pool energies alongside the chemistry.

    Returns
    -------
    dict
        The recorded history, as :meth:`ChannelReactor1D.solve_split_circuit`
        returns it.
    """
    channel = ChannelReactor1D(
        mechanism,
        "plasma",
        r,
        rho_0,
        mtcf,
        names,
        reacting=True,
        kappa=0.0,
        T_amb=Tg_0,
        p_ext=p_0,
        hydrodynamics=hydrodynamics,
        Y_amb=Y_ambient,
    )
    start = time.time()
    history = channel.solve_split_circuit(
        t_out,
        make_circuit(),
        gap,
        np.full(n_radial_points, Tg_0),
        np.full(n_radial_points, Tg_0),
        Y0,
    )
    label = "expanding" if hydrodynamics else "frozen"
    print(f"  {label:<10} done in {time.time() - start:.1f} s")
    return history


print("running the same pulse twice ...")
expanding = solve(hydrodynamics=True)
frozen = solve(hydrodynamics=False)

t_ns = np.asarray(expanding["t"]) * 1e9
r_mm = r * 1e3
Tg = np.asarray(expanding["Tg"])
Te = np.asarray(expanding["Te"])
electron_density = np.asarray(expanding["ne"])
pressure = np.asarray(expanding["p"])
velocity = np.asarray(expanding["u_r"])
density = np.asarray(expanding["rho"])


def deposited(history):
    """Electrical energy delivered to the plasma over the window [J]."""
    return float(
        np.trapezoid(np.asarray(history["I"]) * np.asarray(history["Vp"]), t_out)
    )


print(
    f"  peak |I|  expanding {np.max(np.abs(expanding['I'])):.3f} A, "
    f"frozen {np.max(np.abs(frozen['I'])):.3f} A"
)
print(
    f"  min R_p   expanding {min(expanding['Rp']):.5f} ohm, "
    f"frozen {min(frozen['Rp']):.5f} ohm"
)
print(
    f"  energy    expanding {deposited(expanding) * 1e6:.1f} uJ, "
    f"frozen {deposited(frozen) * 1e6:.1f} uJ "
    f"({100 * (deposited(expanding) / deposited(frozen) - 1):+.1f} %)"
)
print(
    f"  peak Te {Te.max() * u.k_b / u.e:.2f} eV, peak n_e "
    f"{electron_density.max():.3e} m^-3, peak |u_r| {np.abs(velocity).max():.0f} m/s"
)
print(
    f"  axis density {density[0, 0]:.4f} -> {density[-1, 0]:.4f} kg/m^3 "
    f"({100 * (density[-1, 0] / density[0, 0] - 1):+.1f} %)"
)

# %%
# The same circuit seen by a 0D reactor.
# ---------------------------------------
#
# :class:`~rizer.models.nrp.isomass_2T_volume_reactor_cpp.Isomass2TVolumeReactor`
# lumps the channel into one cell of fixed mass and closes the identical
# circuit on :math:`R_p = \mathrm{gap} / (\sigma \pi R^2)`. It is run twice as
# well: at constant volume (``polytropic_index = inf``) and letting the cell
# expand against ambient on a polytropic law. That second run is the lumped
# analogue of the radial sub-step, which is what makes it worth plotting
# beside it.


def solve_0d(polytropic_index):
    """Run the 0D reactor on the same circuit, at a given polytropic index.

    Parameters
    ----------
    polytropic_index : float
        ``numpy.inf`` holds the volume fixed; a finite value lets the cell
        expand against ``p_ext`` on :math:`p V^n = \\mathrm{const}`.

    Returns
    -------
    dict of str to numpy.ndarray
        ``t`` [s], ``Tg``/``Te`` [K], ``ne`` [m^-3], ``Rp`` [ohm], ``I`` [A],
        ``Vp`` [V] and the cell volume ``V`` [m^3], on a 0.1 ns grid.
    """
    cell = ct.Solution(mechanism, "plasma")
    cell.TPX = Tg_0, p_0, partially_ionized_mole_fractions(x_e)
    volume = gap * np.pi * radius_0d**2
    reactor = Isomass2TVolumeReactor(
        cell,
        mechanism,
        "plasma",
        mtcf,
        mass=volume * float(cell.density),
        gap=gap,
        electric_circuit=make_circuit(),
        polytropic_index=polytropic_index,
        p_ext=p_0,
    )
    net = ct.ReactorNet([reactor])
    net.max_order = 5
    net.max_time_step = 1e-10
    net.atol, net.rtol, net.max_steps = 1e-13, 1e-7, 100000

    columns = dict(t=0.0, Tg=Tg_0, Te=Tg_0, ne=ne_0, Rp=0.0, I=0.0, Vp=0.0, V=volume)
    out = {key: [value] for key, value in columns.items()}
    start = time.time()
    for target in np.arange(1e-10, t_end + 1e-12, 1e-10):
        net.advance(target)
        state = np.asarray(reactor.get_state())
        electrons = reactor._r0d.electron_density(
            state[0],
            state[1],
            np.ascontiguousarray(state[3:]),
            reactor.plasma_mass / state[2],
        )
        resistance, voltage = reactor.plasma_resistance, reactor.plasma_voltage
        for key, value in dict(
            t=net.time,
            Tg=state[0],
            Te=state[1],
            ne=electrons,
            Rp=resistance,
            Vp=voltage,
            I=(voltage / resistance if resistance > 0 else 0.0),
            V=state[2],
        ).items():
            out[key].append(value)
    label = "constant volume" if np.isinf(polytropic_index) else "polytropic"
    print(f"  0D {label:<16} done in {time.time() - start:.1f} s")
    return {key: np.asarray(value) for key, value in out.items()}


zero_d_frozen = solve_0d(np.inf)
zero_d_expanding = solve_0d(polytropic_index_0d)

for name, run in (
    ("0D constant volume", zero_d_frozen),
    ("0D polytropic", zero_d_expanding),
):
    energy = float(np.trapezoid(run["I"] * run["Vp"], run["t"]))
    live = run["Rp"] > 0.0
    print(
        f"  {name:<20} peak |I| {np.max(np.abs(run['I'])):.2f} A, "
        f"min R_p {run['Rp'][live].min():.4f} ohm, "
        f"energy {energy * 1e6:.1f} uJ, V/V0 {run['V'][-1] / run['V'][0]:.2f}"
    )

# %%
# The current and the resistance: identical, then not.
# -----------------------------------------------------
#
# The first pulse lands on top of itself. Afterwards the two channels drift
# apart, and the cable echoes arriving at multiples of ``2L/c`` fall on
# measurably different plasmas. The dashed vertical lines are those echo
# arrivals, drawn from the cable length rather than read off the trace.

fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
echoes = np.arange(round_trip, t_end, round_trip) * 1e9

for ax, (label, key, scale) in zip(
    axes,
    (
        ("Current I [A]", "I", 1.0),
        (r"$R_p$ [$\Omega$]", "Rp", 1.0),
        ("Plasma voltage [kV]", "Vp", 1e-3),
    ),
    strict=True,
):
    ax.plot(t_ns, np.asarray(expanding[key]) * scale, label="1D expanding", lw=1.4)
    ax.plot(t_ns, np.asarray(frozen[key]) * scale, "--", label="1D frozen", lw=1.2)
    ax.plot(
        zero_d_expanding["t"] * 1e9,
        zero_d_expanding[key] * scale,
        ":",
        color="0.35",
        lw=1.2,
        label="0D polytropic",
    )
    for echo in echoes:
        ax.axvline(echo, color="0.8", lw=0.6, zorder=0)
    ax.set_xlabel("t [ns]")
    ax.set_ylabel(label)
axes[1].set_yscale("log")
axes[0].legend(fontsize="small")
fig.suptitle("Same circuit, same seed: the expansion only shows up after the pulse")

# %%
# Zoom on the breakdown: the expansion is not there yet.
# -------------------------------------------------------
#
# Over the first 20 ns the two traces are indistinguishable. This is the half
# of the answer that says the radial motion cannot be inferred from the first
# current pulse.

early = t_ns <= 20.0
fig, axes = plt.subplots(1, 2, figsize=(10, 3.6))
axes[0].plot(t_ns[early], np.asarray(expanding["I"])[early], lw=1.6, label="expanding")
axes[0].plot(t_ns[early], np.asarray(frozen["I"])[early], "--", lw=1.2, label="frozen")
axes[0].set_ylabel("Current I [A]")
axes[0].legend(fontsize="small")
relative = (np.asarray(expanding["I"])[early] - np.asarray(frozen["I"])[early]) / max(
    np.max(np.abs(np.asarray(frozen["I"])[early])), 1e-30
)
axes[1].plot(t_ns[early], 100 * relative, color="firebrick", lw=1.4)
axes[1].set_ylabel("Current difference [% of peak]")
for ax in axes:
    ax.set_xlabel("t [ns]")
fig.suptitle("First 20 ns: breakdown finishes before the gas moves")

# %%
# Radial profiles through the pulse and the expansion.
# -----------------------------------------------------
#
# The five fields the channel actually carries, at six times spanning the
# pulse, the relaxation and the rarefaction. The pressure panel is where the
# mechanism is clearest: a core at 20 atm at 10 ns, flat by 500 ns, with the
# velocity panel showing where the gas went.

snapshot_times = [0.0, 10e-9, 50e-9, 150e-9, 300e-9, t_end]
snapshots = sorted(
    {min(int(round(t / dt_out)), len(t_out) - 1) for t in snapshot_times}
)

fields = (
    (r"$T_\mathrm{H}$ [K]", Tg, None),
    (r"$T_\mathrm{e}$ [K]", Te, None),
    (r"$n_\mathrm{e}$ [m$^{-3}$]", electron_density, "log"),
    (r"$p$ [atm]", pressure / ct.one_atm, None),
    (r"$u_r$ [m s$^{-1}$]", velocity, None),
)
fig, axes = plt.subplots(1, 5, figsize=(21, 3.8))
for ax, (label, field, y_scale) in zip(axes, fields, strict=True):
    for frame in snapshots:
        ax.plot(r_mm, field[frame], label=f"{t_ns[frame]:.0f} ns", lw=1.3)
    if y_scale == "log":
        ax.set_yscale("log")
    ax.set_xlabel("r [mm]")
    ax.set_ylabel(label)
axes[0].legend(fontsize="small")
fig.suptitle("Radial profiles of the expanding channel")

# %%
# What the frozen run cannot show.
# ---------------------------------
#
# Overlaying the two at the same late time separates "the model has no flow"
# from "the flow does nothing". The frozen channel holds its density and keeps
# its heat on the axis; the expanding one has moved three quarters of the axis
# mass outward and cooled by doing the work.

late = len(t_out) - 1
frozen_Tg = np.asarray(frozen["Tg"])
frozen_ne = np.asarray(frozen["ne"])
frozen_rho = np.asarray(frozen["rho"])

fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
for ax, (label, a, b, y_scale) in zip(
    axes,
    (
        (r"$T_\mathrm{H}$ [K]", Tg[late], frozen_Tg[late], None),
        (r"$n_\mathrm{e}$ [m$^{-3}$]", electron_density[late], frozen_ne[late], "log"),
        (r"$\rho$ [kg m$^{-3}$]", density[late], frozen_rho[late], None),
    ),
    strict=True,
):
    ax.plot(r_mm, a, lw=1.6, label="expanding")
    ax.plot(r_mm, b, "--", lw=1.3, label="frozen")
    if y_scale == "log":
        ax.set_yscale("log")
    ax.set_xlabel("r [mm]")
    ax.set_ylabel(label)
axes[0].legend(fontsize="small")
fig.suptitle(f"At t = {t_ns[late]:.0f} ns: with and without the radial sub-step")

# %%
# What a lumped expansion cannot do.
# -----------------------------------
#
# The 0D polytropic cell expands against ambient pressure, which is the right
# idea in one dimension too few. Its density falls as ``V0/V``, driven by the
# *mean* pressure of the cell; the 1D core evacuates into a rarefaction running
# inward from a face that is not there in the lumped model. Plotting the two
# density ratios against each other is the sharpest statement of the
# difference, and the right-hand panels show what it costs: the 0D electron
# temperature and resistance follow a cell that never really rarefies.

fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))

axes[0].plot(t_ns, density[:, 0] / density[0, 0], lw=1.6, label="1D, on axis")
axes[0].plot(
    zero_d_expanding["t"] * 1e9,
    zero_d_expanding["V"][0] / zero_d_expanding["V"],
    ":",
    color="0.35",
    lw=1.5,
    label="0D polytropic cell",
)
axes[0].plot(
    zero_d_frozen["t"] * 1e9,
    zero_d_frozen["V"][0] / zero_d_frozen["V"],
    "-.",
    color="0.6",
    lw=1.2,
    label="0D constant volume",
)
axes[0].axhline(1.0, color="0.8", lw=0.6)
axes[0].set_ylabel(r"$\rho / \rho_0$ [-]")
axes[0].legend(fontsize="small")

axes[1].plot(t_ns, Te[:, 0], lw=1.6, label="1D, on axis")
axes[1].plot(
    zero_d_expanding["t"] * 1e9,
    zero_d_expanding["Te"],
    ":",
    color="0.35",
    lw=1.5,
    label="0D polytropic",
)
axes[1].set_ylabel(r"$T_\mathrm{e}$ [K]")
axes[1].legend(fontsize="small")

# `Rp_state`, not `Rp`: the 0D resistance is read after that model's own
# advance, so the comparable 1D quantity is the one belonging to the recorded
# state rather than the one the circuit was driven with over the step.
axes[2].plot(t_ns, np.asarray(expanding["Rp_state"]), lw=1.6, label="1D expanding")
axes[2].plot(t_ns, np.asarray(frozen["Rp_state"]), "--", lw=1.2, label="1D frozen")
axes[2].plot(
    zero_d_expanding["t"] * 1e9,
    zero_d_expanding["Rp"],
    ":",
    color="0.35",
    lw=1.5,
    label="0D polytropic",
)
axes[2].plot(
    zero_d_frozen["t"] * 1e9,
    zero_d_frozen["Rp"],
    "-.",
    color="0.6",
    lw=1.2,
    label="0D constant volume",
)
axes[2].set_yscale("log")
axes[2].set_ylabel(r"$R_p$ [$\Omega$]")
axes[2].legend(fontsize="small")

for ax in axes:
    ax.set_xlabel("t [ns]")
fig.suptitle("Radially resolved against lumped: the same pulse, two geometries")

plt.show()

# %%
# An animation of the expansion.
# -------------------------------
#
# Written to a file rather than embedded: the gallery does not execute the
# ``1dplasma`` examples, so this renders only when the script is run directly.

# sphinx_gallery_thumbnail_number = 3


def make_video(path=None):
    """Animate the radial profiles of the expanding channel over time.

    Parameters
    ----------
    path : str, optional
        Output path. A ``.gif`` of the same stem is written instead when no
        ffmpeg writer is available. The default writes into the system
        temporary directory rather than the working directory: a bare relative
        name drops a megabyte of binary beside the source whenever the example
        is run from the repository, and sphinx-gallery executes example
        scripts with ``__name__ == "__main__"``.

    Returns
    -------
    str or None
        The path actually written, or None when no writer succeeded.
    """
    if path is None:
        path = str(pathlib.Path(tempfile.gettempdir()) / "nrp_expansion.mp4")
    from matplotlib.animation import FFMpegWriter, FuncAnimation, PillowWriter

    panels = (
        (r"$T_\mathrm{H}$ [K]", Tg, None),
        (r"$T_\mathrm{e}$ [K]", Te, None),
        (r"$n_\mathrm{e}$ [m$^{-3}$]", electron_density, "log"),
        (r"$p$ [atm]", pressure / ct.one_atm, None),
        (r"$u_r$ [m s$^{-1}$]", velocity, None),
    )
    fig, axes = plt.subplots(1, 5, figsize=(21, 3.8))
    lines = []
    for ax, (label, field, y_scale) in zip(axes, panels, strict=True):
        (line,) = ax.plot(r_mm, field[0], lw=1.6, color="teal")
        lines.append(line)
        if y_scale == "log":
            ax.set_yscale("log")
            positive = field[field > 0.0]
            ax.set_ylim(positive.min() * 0.5, field.max() * 2.0)
        else:
            margin = 0.05 * (field.max() - field.min())
            ax.set_ylim(field.min() - margin, field.max() + margin)
        ax.set_xlim(0.0, r_mm.max())
        ax.set_xlabel("r [mm]")
        ax.set_ylabel(label)
    title = fig.suptitle("")

    def update(frame):
        """Redraw every panel at output frame `frame`."""
        for line, (_, field, _) in zip(lines, panels, strict=True):
            line.set_ydata(field[frame])
        title.set_text(f"7 kV CH4 channel, expanding:  t = {t_ns[frame]:5.0f} ns")
        return (*lines, title)

    anim = FuncAnimation(fig, update, frames=len(t_out), interval=40, blit=False)
    for writer, out in (
        (FFMpegWriter(fps=25), path),
        (PillowWriter(fps=25), path.replace(".mp4", ".gif")),
    ):
        try:
            anim.save(out, writer=writer)
            plt.close(fig)
            return out
        except Exception as exc:  # noqa: BLE001 - writer availability varies
            print(f"  ({type(writer).__name__} unavailable: {exc})")
    plt.close(fig)
    return None


if __name__ == "__main__":
    written = make_video()
    print(f"animation written to {written}" if written else "no animation writer")