function [f]=rho(S,T,p)
format long

% 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 

% The bulk secant Modulus has already been filled out for you.  
%  You may wish to 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 = ...