rate_equations.py — ODE right-hand side for RadCluster_2_1.
REFERENCE ONLY — not used by the solver; the production RHS is the C++
implementation. May be stale. Scheduled for replacement by the Stage-3
graph-walker. The `ode_system` / `_rhs_case1` / `_rhs_case2` methods below
are dead code: simulation.py integrates exclusively via the C++ solver and
never calls these Python RHS methods. They have drifted from the C++
(Python RHS is mono-defect only) and are kept solely as a physics reference
until the Stage-3 graph-walker replaces them.
Implements the full cluster dynamics master equations (Eqs. 152, 155, 157)
and dispatches to the appropriate He-reduction mode based on physics_option.
Physics reference
-----------------
Ghoniem, N.M. (2026), Sections 5-8 (Rate_Equations.pdf):
Eqs. 152, 155, 157 (master equations)
Eq. 174 (He Case 1 — fusion/fast_eq)
Eq. 175 (He Case 2 — fission/decoupled)
He kinetics (he_kinetics)
-------------------------
'dynamic' — free He c_h integrated as a full ODE variable (Eq. 157)
'quasi_steady_state' — c_h eliminated; computed algebraically from dc_h/dt = 0
at each RHS call (valid because E_m_h = 0.06 eV is small)
State vector y for full_CD modes
---------------------------------
Case 2 — fission (decoupled), dynamic:
y[0..I-1] : SIA clusters c_n, n=1..I
y[I..I+V-1] : void/bubble marginal c_m = Σ_ℓ c_{m,ℓ}, m=1..V
y[I+V] : total He in bubbles Q_tot = Σ_{m,ℓ} ℓ·c_{m,ℓ}
y[I+V+1] : free He c_h
N_eq = I + V + 2
Case 2 — fission (decoupled), quasi_steady_state:
y[0..I-1], y[I..I+V-1], y[I+V] as above
c_h from QSS: c_h = (G_He + β·Q_tot) / (Σ K_HeV·c_v + k2_He)
N_eq = I + V + 1
Case 1 — fusion (mean-field), dynamic:
y[0..I-1] : SIA clusters c_n
y[I..I+V-1] : marginal c_m^tot
y[I+V..I+2V-1] : He content per class Q_m = Σ_ℓ ℓ·c_{m,ℓ}
y[I+2V] : free He c_h
N_eq = I + 2V + 1
Case 1 — fusion (mean-field), quasi_steady_state:
As above without c_h in state; c_h from QSS.
N_eq = I + 2V
Concentration floor (C_floor)
------------------------------
Any state variable below C_floor is clamped to C_floor before evaluating
rate terms. C_floor is set by inp.reactions['C_floor'] (default 1e-15).
