Algorithm Overview
Root solver: The inverter processes each output grid point independently (parallelized with OpenMP). For a target point \(\mathbf{x}_\text{target} = (x_0, x_1, \ldots)\), the inversion solves a nonlinear system of the form:
where \((x_0, x_1, \ldots)\) are the user-configured output grid axes (e.g. \(s, n_B, n_S, n_Q\) or \(e, n_B, y_S, y_Q\)).
The solver uses GSL’s derivative modified-Powell hybrid method (gsl_multiroot_fsolver_hybrids) to iteratively refine \((T, \mu_B, \mu_S, \mu_Q)\) until the residual norm \(|\mathbf{F}|\) falls below tolerance.
This method combines a Newton-type step with a trust-region (dogleg) globalization strategy. The Jacobian
is approximated internally via finite differences of \(\mathbf{F}\).
Seed generation: The seeds (or initial guesses) are generated for each EoS branch. The algorithm retrieves the \(k\) nearest table rows (in output space, measured by scaled Euclidean distance) as initial guesses. The last accepted solution for that branch is also tried as a warm-start seed.
Solution validation: A converged solution is accepted only if
(a) all coordinates lie within the EoS domain,
(b) optional physicality filters pass (\(0 \leq c_s^2 \leq 1\), bounded \(|\mu_i|/T\)), and
(c) the output-space envelope check passes (when enforce_output_envelope is enabled).
Phase selection: When multiple phase branches of the table each return a valid solution, the inverter selects the thermodynamically preferred branch: the mixed (unstable or meta-stable) phase is chosen when exactly one such solution exists; if the table has no phase data, the single stable solution is used directly.
Output: The accepted \((T, \mu_B, \mu_S, \mu_Q)\) and all derived quantities
are written to the output table. Failed points are filled with NaN or excluded,
depending on output.on_inversion_failure.