jnp.interp expects monotonic increasing input. abundances, weak_rates etc interpolate on photon temperature, so this breaks for reheating scenarios. Possible solution is to use interpax.interp1d in everything except thermo (where we're interpolating pretabulated results that we know are monotonic); performance needs to be tested