function rho = rhov(T, P, d)
rho = P./(d.R*T);
