Better interpolation accuracy in cosmology module

At the moment we have a simple linear interpolation for some quantities in the cosmology module. Specifically the function interp_table() does it in a very straightforward way and is used for various quantities in this file.

Is there scope to get additional accuracy here for free by moving to a higher-order interpolation scheme or by using splines without hindering the performance too much?

The expressions for the integrals are in the theory/cosmology sub-directory.