function [rho] = density(p, K, pcav, rhocav)
%Compute fluid density given its Bulk Modulus (K) and pressure. If pressure is
%below cavitation limit density is set to cavitation density
if p > pcav 
    rho = rhocav * exp( (p - pcav) / K);
else, rho = rhocav;
end
end