%Brian Polagye
%January 28, 2019

%Description: Get water properties for experimental work

%Inputs:
%   T: Temperature [C]
%   D: Depth [m]
%   S: Salinity [ppt]
%   P: Pressure [kPa]

%Outputs:
%   rho: Density [kg/m3]
%   nu: Kinematic viscosity [m2/s]
%   c: sound speed (m/s)

function [rho,nu,c] = getWaterProps(T, D, S, P)

%initialize default values
if ~exist('D','var')
    D = 0.3;    %water depth of interest [m]
end

if ~exist('S','var')   
    S = 0;      %salinity [ppt]
end

if ~exist('T','var')
    T = 22;   %temperature [C]
end

if ~exist('P','var')
   P = 101.325 + 1000*9.81*D/1000; %pressure [kPa] (using nominal values)
end

% Salt water - use MIT salt water functions + MacKenzie for sound speed
if S > 0

    %get density
    rho = SW_Density(T,'C',S,'ppt',P,'kPa');
    
    %get kinematic viscosity
    nu = SW_Kviscosity(T,'C',S,'ppt');
    
    %nu = SW_Viscosity(T,'C',S,'ppt')/rho; %same result, because this is all the above function does
    
    %get sound speed
    c = salt_water_c(T,D,S);

% Fresh water - use IAPWS 97 property implementations with addenda for viscosity
else
    
    %get density
    rho = XSteam('rho_pT',P/100,T);    %input T [C] and P [bar]
    
    %get sound speed
    c = XSteam('w_pT',P/100,T);
    
    %get kinematic viscosity
    nu = XSteam('my_pT',P/100,T)/rho;
    
end

end