%%%%%%%%%%%%%
% density.m
%%%%%%%%%%%%%%
function d=density(s,t,p);
% function d=density(s,t,[p]);
%
% calculate density, rho(s,t,p)-1000kg/m^3 %
% s	salinity (PSU)
% t	temperature (C)
% p	pressure (bars - 10m depth = 1 bar; since p=0 for 1 atm) optional
%
% follows steps outlined in Gill appendix 3, pages 599-600 
% equations are labeled with numbers used in Gill's text. 
%

if (nargin==3)
if ((length(s)~=length(t))|(length(s)~=length(p))) 
	error('times series are of unequal length') 
end
else
if (length(s)~=length(t))
error('times series are of unequal length') 
end
end

% calculate powers for polynomials
s32=s.*sqrt(s);
t2=t.*t;t3=t.*t2;t4=t.*t3;

%
% equation for density for pure water (A3.1) 
rhow= 999.842594 + 6.793952e-2*t - 9.095290e-3.*t2 + 1.001685e-4*t3... 
- 1.120083e-6*t4 + 6.536332e-9*(t.*t4);
%
% density at 1 atm (A3.2)
rhost0=rhow + s.*(0.824493 - 4.0899e-3*t + 7.6438e-5*t2... 
- 8.2467e-7*t3 + 5.3875e-9*t4)...
+ s32.*(-5.72466e-3 + 1.0227e-4*t...
- 1.6546e-6*t2) + 4.8314e-4*s.*s;

if (nargin==3)
%
% pure water secant bulk modulus (A3.4)
Kw=19652.21 + 148.4206*t - 2.327105*t2 + 1.360477e-2*t3... 
- 5.155288e-5*t4;

%
% secant bulk modulus at 1 atm (A3.5)
Kst0=Kw + s.*(54.6746 - .603459*t + 1.09987e-2*t2... 
- 6.1670e-5*t3) + s32.*(7.944e-2 + 1.6483e-2*t... 
- 5.3009e-4*t2); 
clear Kw

p2=p.*p;
%
% secant bulk modulus at pressure p (A3.6) 
K=Kst0 + p.*(3.239908 + 1.43713e-3*t...
+ 1.16092e-4*t2 - 5.77905e-7*t3) + (p.*s).*(2.2838e-3... 
- 1.0981e-5*t - 1.6078e-6*t2) + 1.91075e-4*(p.*s32)... 
+ p2.*(8.50935e-5 - 6.12293e-6*t + 5.2787e-8*t2)... 
+ (p2.*s).*(-9.9348e-7 + 2.0816e-8*t + 9.1697e-10*t2); 

% equation for density at pressure p (A3.3) 
d=rhost0./(1-p./K)-1000;
else
d=rhost0-1000;
end
return





