The ρ_p/ρ_t mapping is 1:1, so the interpolation must always preserve
monotonicity, of which cubic splines generally make no guarantee.
Note: Linear interpolation does not provide even C¹ continuity, but
these data is not directly used in the numerical integration, so it
should be fine. Ideally this should be replaced with cubic splines
computed with the Fritsch–Carlson algorithm.
This adds a new `splines` module which implements a high-level interface
for creating and evaluating splines and rewrite almost all modules to
use it. Also, notably:
1. both `simplespline` and DIERCKX splines can now used with a uniform
interface
2. most complexity due to handling working space arrays is gone
3. memory management has been significantly simplified too
