function f=rho(S,T,P)
% Lab 2: rho.m
% Purpose: Creates a function to calculate density
%

% Variables:
% S = Salinity
% T = Temperature
% P = Pressure
% m = meters
%
%--------------------------------------------------
% display more than four digits after decimal
format long
%
% Check input!!!
% input: rho(S,T,P)
% OR
% input: rho(S, T, m)
% MUST change function line (p to m) and
% need to convert meters to bars (i.e. uncomment next line)
% p=m/10;

% The coefficeints of the interantional equation of state have
% been provided here

r=[9.99842594E+02
6.79395200E-02
-9.09529000E-03
1.00168500E-04
-1.12008300E-06
6.53633200E-09
8.24493000E-01
-4.08990000E-03
7.64380000E-05
-8.24670000E-07
5.38750000E-09
-5.72466000E-03
1.02270000E-04
-1.65460000E-06
4.83140000E-04];

c=[1.965221E+04
1.484206E+02
-2.327105E+00
1.360477E-02
-5.155288E-05
3.239908E+00
1.437130E-03
1.160920E-04
-5.779050E-07
8.509350E-05
-6.122930E-06
5.278700E-08
5.467460E+01
-6.034590E-01
1.099870E-02
-6.167000E-05
7.944000E-02
1.648300E-02
-5.300900E-04
2.283800E-03
-1.098100E-05
-1.607800E-06
1.910750E-04
-9.934800E-07
2.081600E-08
9.169700E-10];

% input rho_o given the above coefficients, r(i), above 

rho_o = r(1)+r(2)*T ...
+r(3)*T^2+r(4)*T^3 ...
+r(5)*T^4+r(6)*T^5 ...
+r(7)*S+r(8)*T*S ...
+r(9)*T^2*S+r(10)*T^3*S ...
+r(11)*T^4*S+r(12)*S^(3/2) ...
+r(13)*T*S^(3/2)+r(14)*(T^2)*(S^(3/2)) ...
+r(15)*S^2;

% The bulk secant Modulus has already been filled out for you.  
% You may use it as a guide to input rho_o

K = c(1) ...
+c(2)*T + c(3)*T^2 ...
+c(4)*T^3 + c(5)*T^4 ...
+c(6)*p + c(7)*T*p ...
+c(8)*T^2*p + c(9)*T^3*p ...
+c(10)*p^2 + c(11)*T*p^2 ...
+c(12)*T^2*p^2 ...
+c(13)*S+c(14)*T*S ...
+c(15)*T^2*S+c(16)*T^3*S ...
+c(17)*S^(3/2)+c(18)*T*S^(3/2) ...
+c(19)*T^2.*S^(3/2)+c(20)*p*S ...
+c(21)*T*p*S+c(22)*T^2*p*S ...
+c(23)*p*S^(3/2)+c(24)*p^2*S ...
+c(25)*T*p^2*S+c(26)*T^2*p^2*S;

% Given rho_o and K, define the density field here

f=rho_o/(1-P/K);