% Function to calculate density through non-linear formulation of Millero and Poisson (1981)
% relative density (rho-rho0= AS+BS^(3/2)+CS^2:
% A=A1+A2*t+A3*t^2+A4*t^3+A5*t^4
% B=B1+B2*t+B3*t^2
% C=C1

function [drho]=density(temp,salt)

persistent dens_con

if isempty(dens_con)
    dens_con=[273.15;8.24493*10^-1;-4.0899*10^-3;7.6438*10^-5;-8.2467*10^-7;5.3875*10^-9;-5.72466*10^-3;1.0227*10^-4;-1.6546*10^-6;4.8314*10^-4];
end

temp=temp-dens_con(1);

A=dens_con(2)+dens_con(3)*temp+dens_con(4)*temp^2+dens_con(5)*temp^3+dens_con(6)*temp^4;
B=dens_con(7)+dens_con(8)*temp+dens_con(9)*temp^2;

drho=A*salt+B*(max(0,salt))^(3/2)+dens_con(10)*salt^2;

end
