EOS: clarify the meaning of ieos 13 and 14 in sources & maybe ensure that R_ref = R_0 = 1

[tdavidcl]

While matching Shamrock and Phantom for ieos 13 and 14 I noticed some issues in docstrings and maybe with the disc setup (maybe only in an edge case where R_ref != R_in).

Currently the code refer to the Farris 2014 EOS which only corresponds to q=1/2

$$c_s = \frac{H_0}{r_0}\left(\frac{G M_1}{r_1} + \frac{G M_2}{r_2}\right)$$

However the extension of that EOS to q != 1/2 was only introduced in Ragussa et al 2016 if I'm right with

$$c_s = \frac{H_0}{r_0} \left( \frac{G M_1}{r_1} + \frac{G M_2}{r_2} \right)^{q}$$

But as is the units are broken if q is not 1/2 so you need to compensate with $r_0 \Omega_0$

$$c_s = \frac{H_0}{r_0} \frac{1}{(r_0 \Omega_0)^{2q - 1}} \left( \frac{G M_1}{r_1} + \frac{G M_2}{r_2} \right)^{q}$$ $$= c_{s0} \frac{1}{(r_0 \Omega_0)^{q}} \left( \frac{G M_1}{r_1} + \frac{G M_2}{r_2} \right)^{q}$$ $$= c_{s0} \frac{1}{r_0^{q}} \left[ \frac{1}{\sum_i M_i} \sum_i \frac{M_i}{r_i} \right]^{q}$$

So in phantom polyk should be $(c_{s0}/{r_0^{q}})^2$ rather than

phantom/src/setup/set_disc.f90

Line 269 in f6d5bee

cs0 = H_R*sqrt(G*star_m/R_ref)*R_ref**q_index

phantom/src/setup/set_disc.f90

Line 270 in f6d5bee

polyk = cs0**2

So it looks to me that in phantom the current setup is correct only if R_ref = R_0 = 1

I can open a PR if you want & confirm this