% This function will plot contours of the density of water vs. temperature
%   and salt content.
% theta is the temperature (in degrees C) of the solution.
% wt is the fraction weight of salt over the total weight of the solution.
% Maximum error is 5*10^(-5) for 10<theta<60, wt<0.16 and error is <5*10^(-4)
%   for 0<theta<60, wt<0.2.
% Equation taken from Ruddick and Shirtcliffe, 1979.

close all;
clear;

theta = [0:0.1:60];
wt = [0:0.005:0.2];
theta1 = ones(length(wt),1);
wt1 = ones(length(theta),1);

theta = theta1*theta;
wt = (wt1*wt)';

T = theta - 25;

t1 = -2.539*10^(-4);
t2 = -4.968*10^(-6);
t3 = 2.7*10^(-8);

w1 = 0.69976;
w2 = 0.14042;
w3 = 0.337;

w1t1 = 1.6803*10^(-3);
w1t2 = -3.551*10^(-5);
w1t3 = 3.52*10^(-7);

w2t1 = 2.714*10^(-3);
w2t2 = -8.11*10^(-5);
w2t3 = 9.0*10^(-7);

rho0 = 0.99708;

rho = rho0 + t1*T + t2*T.^2 + t3*T.^3 + w1*wt + w2*wt.^2 + w3*wt.^3 - (w1t1*T + w1t2*T.^2 + w1t3*T.^3).*wt + ...
    (w2t1*T + w2t2*T.^2 + w2t3*T.^3).*wt.^2;

contour(wt,theta,rho,450);
colorbar('vert');
xlabel('Salt weight fraction');
ylabel('Temperature (^oC)');
title('Density (10^{-3} kg/m^3)');

