function OutVar = H2Oprop(query, Var1, Var2, Var3)

%
% H2Oprop.m v1.01
%
% Developed by Bryan Ruddy (ruddy AT mit.edu), 9/2006 - 5/2009
%
% The most recent version of this function can be found at http://web.mit.edu/ruddman/Public/Thermodynamics/H2Oprop.m
%
% This work is licensed under the Creative Commons Attribution 3.0 United 
% States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by/3.0/us/ or send a letter to Creative
% Commons, 171 Second Street, Suite 300, San Francisco, California, 94105,
% USA.
%
% This script computes the thermodynamic and transport properties of water
% using the latest (December 2008) IAPWS scientific formulations, with 
% rigorous bounds-checking.  Information on the range of validity and 
% accuracy of these properties can be found in the IAPWS releases at 
% www.iapws.org.  A total of 70 different combinations of input and output
% parameters are supported, including 45 basic thermodynamic relationships,
% 14 saturation relationships, and 11 other thermodynamic and transport properties.
%
% All quantities are represented in SI units.
%
% Usage:
% The function should be called with a case-insensitive text string (query)
% and up to three scalar numerical parameters.  The identities of these 
% parameters depend on the choice of query, as described below.  The
% function outputs a single scalar value, whose meaning also depends on the
% choice of query.  Note that queries in terms of density and temperature
% are the most computationally efficient; others use root-finding and/or
% minimization routines.
%
% Query Options:
%
% t_dp  Temperature from density and pressure (gives high-temperature point near density anomaly)
%
% t_dh  Temperature from density and specific enthalpy
%
% t_ph  Temperature from pressure and specific enthalpy
%
% t_du  Temperature from density and specific internal energy
%
% t_pu  Temperature from pressure and specific internal energy
%
% t_ds  Temperature from density and specific entropy
%
% t_ps  Temperature from pressure and specific entropy
%
% t_hs  Temperature from specific enthalpy and specific entropy
%
% t_us  Temperature from specific internal energy and specific entropy
%
% p_dt  Pressure from density and temperature
%
% p_dh  Pressure from density and specific enthalpy
%
% p_du  Pressure from density and specific internal energy
%
% p_ds  Pressure from density and specific entropy
%
% p_ts  Pressure from temperature and specific entropy
%
% p_hs  Pressure from specific enthalpy and specific entropy
%
% p_us  Pressure from specific internal energy and specific entropy
%
% d_pt  Density from pressure and temperature (gives liquid value if saturated)
%
% d_ph  Density from pressure and specific enthalpy
%
% d_pu  Density from pressure and specific internal energy
%
% d_ps  Density from pressure and specific entropy
%
% d_ts  Density from temperature and specific entropy
%
% d_hs  Density from specific enthalpy and specific entropy
%
% d_us  Density from specific internal energy and specific entropy
%
% h_dt  Specific enthalpy from density and temperature
%
% h_pt  Specific enthalpy from pressure and temperature (gives liquid value if saturated)
%
% h_du  Specific enthalpy from density and specific internal energy
%
% h_ds  Specific enthalpy from density and specific entropy
%
% h_pu  Specific enthalpy from pressure and specific internal energy
%
% h_ps  Specific enthalpy from pressure and specific entropy
%
% h_ts  Specific enthalpy from temperature and specific entropy
%
% h_us  Specific enthalpy from specific internal energy and specific entropy
%
% u_dt  Specific internal energy from density and temperature
%
% u_dh  Specific internal energy from density and specific enthalpy
%
% u_ds  Specific internal energy from density and specific entropy
%
% u_pt  Specific internal energy from pressure and temperature (gives liquid value if saturated)
%
% u_ph  Specific internal energy from pressure and specific enthalpy
%
% u_ps  Specific internal energy from pressure and specific entropy
%
% u_ts  Specific internal energy from temperature and specific entropy
%
% u_hs  Specific internal energy from specific enthalpy and specific entropy
%
% s_dt  Specific entropy from density and temperature
%
% s_dh  Specific entropy from density and specific enthalpy
%
% s_du  Specific entropy from density and specific internal energy
%
% s_pt  Specific entropy from pressure and temperature (gives liquid value if saturated)
%
% s_ph  Specific entropy from pressure and specific enthalpy
%
% s_pu  Specific entropy from pressure and specific internal energy
%
% psat_t   Saturation pressure from temperature
%
% tsat_p   Saturation temperature from pressure
%
% tsat_dl  Saturation temperature from liquid density (gives high-temperature point near density anomaly)
%
% tsat_dv  Saturation temperature from vapor density
%
% dlsat_t  Saturated liquid density from temperature
%
% dlsat_p  Saturated liquid density from pressure
%
% dvsat_t  Saturated vapor density from temperature
%
% dvsat_p  Saturated vapor density from pressure
%
% hlsat_t  Saturated liquid specific enthalpy from temperature
%
% hlsat_p  Saturated liquid specific enthalpy from pressure
%
% hvsat_t  Saturated vapor specific enthalpy from temperature
%
% hvsat_p  Saturated vapor specific enthalpy from pressure
%
% ulsat_t  Saturated liquid specific internal energy from temperature
%
% ulsat_p  Saturated liquid specific internal energy from pressure
%
% uvsat_t  Saturated vapor specific internal energy from temperature
%
% uvsat_p  Saturated vapor specific internal energy from pressure
%
% slsat_t  Saturated liquid specific entropy from temperature
%
% slsat_p  Saturated liquid specific entropy from pressure
%
% svsat_t  Saturated vapor specific entropy from temperature
%
% svsat_p  Saturated vapor specific entropy from pressure
%            
% x_dt        Steam quality from density and temperature
%
% cp_dt       Constant pressure specific heat from density and temperature
%
% cv_dt       Constant volume specific heat from density and temperature
%
% betat_dt    Thermal expansion coefficient from density and temperature
%
% kappat_dt   Isothermal compressibility from density and temperature
%
% w_dt        Sound velocity from density and temperature
%
% mu_dt       Viscosity from density and temperature
%
% k_dt        Thermal conductivity from density and temperature
%
% sigma_t     Surface tension from temperature
%
% epsilon_dt  Dielectric constant from density and temperature
%
% n_dtl       Refractive index from density, temperature, and wavelength

% Programming Notes:
%
% Code known to be problematic, for efficiency and/or accuracy reasons, is
% marked with a *** in the comments.  This code should be fixed in future
% versions.
%
% In general, a major objective in future revisions should be to reduce the
% degree to which the code is recursive.  This will increase performance
% and reliability by eliminating unneccessary bounds checks and reducing
% the use of nested fzero calls.
%
% Change Log:
%
% v1.01: Replaced recursive function references with references to this
%        function's current name, rather than the development name.

%---------------------------------------------------------------------------------------  
% Equation of State Constants

% 1995 IAPWS Scientific Formulation

global Tc Tt Pc rhoc R n0 gamma0 c d t a b B C D A alpha beta gamma epsilon n

Tc   = 647.096;					%Critical Temperature, K
rhoc = 322;						%Critical Density, kg/m^3
R    = 461.51805;               %Water Gas Constant, J/kg-K

n0 = [-8.32044648374976; 6.68321052759325; 3.00632; 0.012436; 0.97315; 1.27950; ...
      0.96956; 0.24873];
gamma0 = [0; 0; 0; 1.28728967; 3.53734222; 7.74073708; 9.24437796; 27.5075105];

c = [zeros(7,1); ones(15,1); 2*ones(20,1); 3; 3; 3; 3; 4; 6; 6; 6; 6; zeros(5,1)];
d = [1; 1; 1; 2; 2; 3; 4; 1; 1; 1; 2; 2; 3; 4; 4; 5; 7; 9; 10; 11; 13; 15; 1; 2;...
     2; 2; 3; 4; 4; 4; 5; 6; 6; 7; 9*ones(5,1); 10; 10; 12; 3; 4; 4; 5; 14; 3; 6;...
     6; 6; 3; 3; 3; 0; 0];
t = [-0.5; 0.875; 1; 0.5; 0.75; 0.375; 1; 4; 6; 12; 1; 5; 4; 2; 13; 9; 3; 4; 11;...
     4; 13; 1; 7; 1; 9; 10; 10; 3; 7; 10; 10; 6; 10; 10; 1; 2; 3; 4; 8; 6; 9; 8;...
     16; 22; 23; 23; 10; 50; 44; 46; 50; 0; 1; 4; 0; 0];
a = [zeros(54,1); 3.5; 3.5];
b = [zeros(54,1); 0.85; 0.95];
B = [zeros(54,1); 0.2; 0.2];
C = [zeros(54,1); 28; 32];
D = [zeros(54,1); 700; 800];
A = [zeros(54,1); 0.32; 0.32];
alpha   = [zeros(51,1); 20; 20; 20; 0; 0];
beta    = [ones(51,1); 150; 150; 250; 0.3; 0.3];
gamma   = [zeros(51,1); 1.21; 1.21; 1.25; 0; 0];
epsilon = [zeros(51,1); 1; 1; 1; 0; 0];
n = ...
[0.12533547935523e-1   0.78957634722828e1  -0.87803203303561e1   0.31802509345418 ...
-0.26145533859358     -0.78199751687981e-2  0.88089493102134e-2 -0.66856572307965 ...
 0.20433810950965     -0.66212605039687e-4 -0.19232721156002    -0.25709043003438 ...
 0.16074868486251     -0.40092828925807e-1  0.39343422603254e-6 -0.75941377088144e-5 ...
 0.56250979351888e-3  -0.15608652257135e-4  0.11537996422951e-8  0.36582165144204e-6 ...
-0.13251180074668e-11 -0.62639586912454e-9 -0.10793600908932     0.17611491008752e-1 ...
 0.22132295167546     -0.40247669763528     0.58083399985759     0.49969146990806e-2 ...
-0.31358700712549e-1  -0.74315929710341     0.47807329915480     0.20527940895948e-1 ...
-0.13636435110343      0.14180634400617e-1  0.83326504880713e-2 -0.29052336009585e-1 ...
 0.38615085574206e-1  -0.20393486513704e-1 -0.16554050063734e-2  0.19955571979541e-2 ...
 0.15870308324157e-3  -0.16388568342530e-4  0.43613615723811e-1  0.34994005463765e-1 ...
-0.76788197844621e-1   0.22446277332006e-1 -0.62689710414685e-4 -0.55711118565645e-9 ...
-0.19905718354408      0.31777497330738    -0.11841182425981    -0.31306260323435e2 ...
 0.31546140237781e2   -0.25213154341695e4  -0.14874640856724     0.31806110878444]';

%---------------------------------------------------------------------------------------  
% Phase Boundary Constants

% Used in a variety of bounds-checking, root-finding, and minimization
% tasks throughout the code.

Pc = 22.064e6;                  % Critical Pressure (IAPWS95)
Pt = 611.654771007894;          % Triple Point Pressure (IAPWS95)
%Pmax = 1e9;                     % Maximum Pressure for Experimental Verification (IAPWS95)
Picemax = 20617812820.4493;     % Maximum Pressure on Melting Curve (715 K) (IAPWS08)
Pmin = 1.9349587189864e-040;    % Minimum Pressure on Sublimation curve (50 K) (IAPWS08)

Tt = 273.16;                    % Triple Point Temperature (IAPWS95)
Tmax = 1273;                    % Maximum Temperature for Experimental Verification (IAPWS95)
Tmin = 251.165;                 % Minimum Liquid Temperature (Ice-Ice-Liquid Triple Point) (IAPWS08)
Trhomax = 277.148195578572;     % Highest Temperature of Density Maximum Anomaly
Tabsmin = 50;                   % Minimum Temperature for Experimental Verification (IAPWS08)
Ticemax = 715;                  % Maximum Freezing Temperature (IAPWS08)
Tsmin = 694.274446246169;       % Temperature of Local Minimum Freezing Entropy
Tsmax = 603.867446392581;       % Temperature of Maximum Freezing Entropy

rhovt   = 0.00485457572477878;  % Triple Point Vapor Density
rholt   = 999.792520031622;     % Triple Point Liquid Density
rhoamax = 999.925131084086;     % Saturation Density Anomaly Point 
rhomax = 1237.39112310848;      % Maximum Density (1GPa on melting curve)
rholTmin = 1091.68582714549;    % Liquid Density at Tmin
rhomin = 3.29347680703935e-46;  % Minimum Density (Pmin at Tmax)
rhoicemax = 1930.89627377429;   % Maximum Known Liquid Density
rhosatmin = 9.86827543571457e-045;  % Minimum Phase Boundary Density (50 K)

hlt = 0.611781702320426;        % Saturated Liquid Enthalpy at Triple Point
hvt = 2500915.19146571;         % Saturated Vapor Enthalpy at Triple Point
ha  = 16805.6061030613;         % Saturated Liquid Enthalpy at Density Anomaly
hvsatmin = 2088366.77643276;    % Saturated Vapor Enthalpy at 50 K

uvsatmin = 2065290.87393276;    % Saturated Vapor Energy at 50 K 
ult = -0.000236902622937808;    % Saturated Liquid Energy at Triple Point
uvt = 2374919.6757225;          % Saturated Vapor Energy at Triple Point
ua = 16804.792597053;           % Saturated Liquid Energy at Density Anomaly

svsatmin = 51104.228263276876532; % Saturated Vapor Entropy at 50 K
slt = 2.72849374306335e-012;    % Saturated Liquid Entropy at Triple Point
svt = 9155.49340929858;         % Saturated Vapor Entropy at Triple Point
sa = 61.0753924234768;          % Saturated Liquid Entropy at Density Anomaly
spimax = 1300.31975929576;      % Highest Liquid Entropy on Pressurized Ice Line
spitmax = 1246.430489334724371; % Liquid Entropy at High End of Pressurized Ice Line
spimin = -337.232036172874;     % Lowest Liquid Entropy (at Tmin)
spitmin = 1213.949641489883561; % Liquid Entropy Minimum near Top of Hot Ice Line
scrit = 4407.116565421005362;   % Entropy at Critical Point
sdtmax = 3267.889828921342541;  % Entropy at Tmax and rhoicemax

%---------------------------------------------------------------------------------------   
% Function Body

query = lower(query);   % Strip out case from input parameter
switch query
    case 'p_dt'         % Pressure from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Tc)                                 % Supercritical
            P = Water_P(rho, T);                        % Compute directly from correlation
        elseif (T >= Tt)                            % Liquid - Vapor Regime
            [Psat rholsat rhovsat] = GoodSatProp(T);    % Find accurate saturation properties
            if (rho > rhovsat) && (rho < rholsat)           % Two-Phase
                P = Psat;                                       % Clearly the fluid is saturated...
            else                                            % Liquid or Vapor
                P = Water_P(rho, T);                            % We aren't saturated, so directly compute pressure
                if P > MeltPress(T)                                 % High-Pressure Ice
                    P = NaN;                                            % There's no correlation for this stuff...
                end
            end
        elseif (T >= Tmin)                          % Ice - Liquid - Vapor Regime
            [Pm Ps Pf] = MeltPress(T);                  % Find melting, freezing, and sublimation pressures
            rhovsub = fzero(@(x) Water_P(x, T) - Ps, Ps/R/T); % Find vapor density at sublimation point
            rholmin = fzero(@(x) Water_P(x, T) - Pm, rhomax); % Find liquid density at ice Ih melting point
            rholmax = fzero(@(x) Water_P(x, T) - Pf, rhomax); % Find liquid density at high-pressure ice melting point
            if (rho <= rhovsub) || ((rho >= rholmin) && (rho <= rholmax))   %Liquid or Gas
                P = Water_P(rho, T);                            % Single-phase, so directly compute pressure
            else                                            % Ice
                P = NaN;                                        % Would need the ice correlation to compute this...
            end
        else                                        % Ice - Vapor Regime
            Ps = MeltPress(T);                          % Find sublimation pressure
            rhovsat = fzero(@(x) Water_P(x, T) - Ps, Ps/R/T); % Find saturated vapor density
            if (rho <= rhovsat)                             % Vapor
                P = Water_P(rho, T);                            % Single-phase, so directly compute pressure
            else                                            % Ice
                P = NaN;                                        % Would need the ice correlation to compute this...
            end
        end
        OutVar = P;         % Output the pressure
        
%---------------------------------------------------------------------------------------
        
    case 'd_pt'         % Density from pressure, temperature
        P = Var1;           % Input pressure
        T = Var2;           % Input temperature
        Tm = MeltTemp(P);   % Compute solid-fluid phase transition temperature
        if (T < (Tm - 2*eps(Tm))) || isnan(Tm)                          % Ice or Absurd Pressure
            rho = NaN;                                                      % Would need the ice correlation to compute this...
        elseif T > Tc                                                   % Supercritical
            rho = fzero(@(x) Water_P(x, T) - P, [rhomin*.999 rhoicemax*1.01]); % Use root-finding method to compute density
        elseif T > Tt                                                   % Liquid - Vapor Regime
            [Psat rholsat rhovsat] = GoodSatProp(T);                        % Find accurate saturation properties
            if (P >= Psat)                                                      % Liquid
                rho = fzero(@(x) Water_P(x, T) - P, [rholsat*.999 rhoicemax]);      % Use root-finding method, and give liquid density if saturated
            else                                                                % Vapor
                rho = fzero(@(x) Water_P(x, T) - P, [rhomin rhovsat]);              % Use root-finding method to compute density
            end
        elseif T >= Tabsmin                                             % Liquid - Ice - Vapor Regime
            if (P < Pt)                                                         % Thin Vapor
                rho = fzero(@(x) Water_P(exp(x), T) - P, log(P/R/T));               % Use root-finding method to compute density
                rho = exp(rho);                                                     % Due to the extreme pressure range, examine the logarithm of density
            else                                                                % Cold Liquid
                rho = fzero(@(x) Water_P(x, T) - P, rhomax);                        % Use root-finding method to compute density
            end
        else                                                            % Too Cold
            rho = NaN;                                                      % Vapor is effectively nonexistent below Tabsmin...
        end
        OutVar = rho;       % Output density
                
%---------------------------------------------------------------------------------------   
        
    case 't_dp'         % Temperature from density, pressure, only valid above density anomaly
        rho = Var1;         % Input density
        P   = Var2;         % Input pressure
        Tm = MeltTemp(P);   % Melting temperature
        if (rho >= rholt) && (P < 3e7)                                  % Density anomaly region
            [Trhomax1, rhomax1] = fminbnd(@(x) 1/H2Oprop('d_pt',P,x), ...
                                  Tm, Trhomax*1.001, optimset('TolX',1e-15)); % Find density anomaly
            rhomax1 = 1/rhomax1;                                            % Correct the maximum density
            if (rho > rhomax1)                                          % Non-Physical
                T = NaN;                                                    % At a given pressure, ice at any temperature floats...
            else                                                        % Liquid
                Tsat = H2Oprop('tsat_p',P);                              % Find the saturation temperature
                T = fzero(@(x) Water_P(rho,x) - P, [Trhomax1 Tsat]);        % Use root-finding method to compute temperature
            end
        else                                                            % No density anomaly, skip the minimization
            Tsat = H2Oprop('tsat_p',P);                                  % Find the saturation temperature from pressure
            [junk rhol rhov] = GoodSatProp(Tsat);                           % Find the saturation densities
            rhom = H2Oprop('d_pt',P,Tm);
            if (rho >= rhom)                                            % Ice or Non-Physical
                T = NaN;
            elseif (rho >= rhol)                                        % Liquid
                Tmin1 = H2Oprop('tsat_dl',rho);                          % Find the saturation temperature from density
                if isnan(Tmin1)                                             % Check to see if that density is actually on the saturation line
                    Tmin1 = Tm;                                             % If not, use the melting temperature as the lower limit instead
                end
                T = fzero(@(x) Water_P(rho,x) - P, [Tmin1 Tsat]);           % Use root-finding method to compute temperature
            elseif (rho <= rhov)                                        % Vapor
                T = fzero(@(x) Water_P(rho,x) - P, [Tsat Tmax]);            % Use root-finding method to compute temperature
            elseif (P > Pc)                                             % Supercritical fluid
                if rho > rhoc                                               % Liquid-like fluid
                    Tmin1 = H2Oprop('tsat_dl',rho);                      % Find the saturation temperature from density
                else                                                        % Vapor-like fluid
                    Tmin1 = H2Oprop('tsat_dv',rho);                      % Find the saturation temperature from density
                end
                if isnan(Tmin1)                                             % Find the saturation temperature from density
                    Tmin1 = Tm;                                             % If not, use the melting temperature as the lower limit instead
                end
                T = fzero(@(x) Water_P(rho,x) - P, [Tmin1 Tmax]);           % Use root-finding method to compute temperature
            else                                                        % Two-phase fluid
                T = Tsat;                                                   % Fluid is at saturation temperature
            end
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   

    case 'u_dt'         % Specific energy from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T < Tmin)           % Ice or Vapor
            Pm = MeltPress(T);      % Find the sublimation pressure from temperature
            rhos = H2Oprop('d_pt',Pm,T); % Compute the vapor density at the sublimation point
            if (rho > rhos)             % Ice
                u = NaN;                    % Would need the ice correlation to compute this...
            else                        % Vapor
                u = Water_u(rho, T);        % Directly compute the specific energy
            end
        elseif (T < Tt)         % Ice, Liquid, or Vapor
            [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the Ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density at the high-pressure ice melting point
            if (rho < rhos)             % Vapor
                u = Water_u(rho, T);        % Directly compute the specific energy
            elseif (rho < rhom)         % Ice Ih
                u = NaN;                    % Would need the ice correlation to compute this...
            elseif (rho < rhof)         % Liquid
                u = Water_u(rho, T);        % Directly compute the specific energy
            else                        % High-Pressure Ice
                u = NaN;                    % There's no correlation for this stuff...
            end
        else                    % Liquid or Vapor
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhovsat) && (rho < rholsat) % Two-phase Mixture
                v  = 1/rho;             % Convert density to specific volume
                vl = 1/rholsat;         % Convert saturated liquid density to specific volume
                vv = 1/rhovsat;         % Convert saturated vapor density to specific volume
                X  = (v - vl)/(vv - vl); % Compute the quality from specific volume
                ul = Water_u(rholsat, T); % Directly compute the saturated liquid specific energy
                uv = Water_u(rhovsat, T); % Directly compute the saturated vapor specific energy
                u = X*uv + (1-X)*ul;    % Find the mixture specific energy
            else                    %Single-phase
                u = Water_u(rho, T);    % Directly compute the specific energy
            end
        end
        OutVar = u;         % Output specific internal energy
        
%---------------------------------------------------------------------------------------   

    case 's_dt'         % Specific entropy from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T < Tmin)           % Ice or Vapor
            rhos = H2Oprop('dvsat_t',T); % Compute the vapor density at the sublimation point
            if (rho > rhos)             % Ice
                s = NaN;                    % Would need the ice correlation to compute this...
            else                        % Vapor
                s = Water_s(rho, T);        % Directly compute the specific entropy
            end
        elseif (T < Tt)         % Ice, Liquid, or Vapor
            [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the Ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density at the high-pressure ice melting point
            if (rho < (rhos+10*eps(rhos))) % Vapor
                s = Water_s(rho, T);        % Directly compute the specific entropy
            elseif (rho < rhom)         % Ice Ih
                s = NaN;                    % Would need the ice correlation to compute this...
            elseif (rho <= rhof)         % Liquid
                s = Water_s(rho, T);        % Directly compute the specific entropy
            else                        % High-Pressure Ice
                s = NaN;                    % There's no correlation for this stuff...
            end
        else                    % Liquid or Vapor
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhovsat) && (rho < rholsat) % Two-phase Mixture
                v  = 1/rho;             % Convert density to specific volume
                vl = 1/rholsat;         % Convert saturated liquid density to specific volume
                vv = 1/rhovsat;         % Convert saturated vapor density to specific volume
                X  = (v - vl)/(vv - vl); % Compute the quality from specific volume
                sl = Water_s(rholsat, T); % Directly compute the saturated liquid specific entropy
                sv = Water_s(rhovsat, T); % Directly compute the saturated vapor specific entropy
                s = X*sv + (1-X)*sl;    % Find the mixture specific entropy
            else                    %Single-phase
                s = Water_s(rho, T);    % Directly compute the specific entropy
            end
        end
        OutVar = s;         % Output specific entropy
        
%---------------------------------------------------------------------------------------   

    case 'h_dt'         % Specific enthalpy from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T < Tmin)           % Ice or Vapor
            Pm = MeltPress(T);      % Find the sublimation pressure from temperature
            rhos = H2Oprop('d_pt',Pm,T); % Compute the vapor density at the sublimation point
            if (rho > (rhos+10*eps(rhos))) % Ice
                h = NaN;                    % Would need the ice correlation to compute this...
            else                        % Vapor
                h = Water_h(rho, T);        % Directly compute the specific enthalpy
            end
        elseif (T < Tt)         % Ice, Liquid, or Vapor
            [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the Ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density at the high-pressure ice melting point
            if (rho < rhos)             % Vapor
                h = Water_h(rho, T);        % Directly compute the specific enthalpy
            elseif (rho < rhom)         % Ice Ih
                h = NaN;                    % Would need the ice correlation to compute this...
            elseif (rho <= rhof)         % Liquid
                h = Water_h(rho, T);        % Directly compute the specific enthalpy
            else                        % High-Pressure Ice
                h = NaN;                    % There's no correlation for this stuff...
            end
        else                    % Liquid or Vapor
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhovsat) && (rho < rholsat) % Two-phase Mixture
                v  = 1/rho;             % Convert density to specific volume
                vl = 1/rholsat;         % Convert saturated liquid density to specific volume
                vv = 1/rhovsat;         % Convert saturated vapor density to specific volume
                X  = (v - vl)/(vv - vl); % Compute the quality from specific volume
                hl = Water_h(rholsat, T); % Directly compute the saturated liquid specific enthalpy
                hv = Water_h(rhovsat, T); % Directly compute the saturated vapor specific enthalpy
                h = X*hv + (1-X)*hl;    % Find the mixture specific enthalpy
            else                    %Single-phase
                h = Water_h(rho, T);    % Directly compute the specific enthalpy
            end
        end
        OutVar = h;         % Output specific enthalpy
        
%---------------------------------------------------------------------------------------   
        
    case 't_dh'         % Temperature from density, enthalpy
        rho = Var1;         % Input density
        h   = Var2;         % Input specific enthalpy
        if (rho < rhomin)                                               % Too rarefied
            T = NaN;                                                        % Densities this low violate continuum approximations anyway...
        elseif (rho < rhosatmin)                                        % Extremely Thin Vapor (less dense than 50 K saturated vapor)
            if (h >= hvsatmin)                                              % Vapor
                T = fzero(@(x) Water_h(rho, x) - h, [Tabsmin Tmax]);            % Use root-finding method to compute temperature
            else                                                            % Non-Physical
                T = NaN;                                                        % It's not ice at this density...
            end
        elseif (rho < rhovt)                                            % Thin Vapor (less dense than triple point vapor)
            Tsat = H2Oprop('tsat_dv',rho);                               % Find saturation temperature from density
            hsat = Water_h(rho, Tsat);                                      % Find saturated vapor enthalpy at this density
            if (h >= hsat)                                                      % Vapor
                T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]);                   % Use root-finding method to compute temperature
            else                                                                % Two-Phase with Ice
                T = NaN;                                                            % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoc)                                             % Vapor or Two-Phase with Liquid
            Tsat = H2Oprop('tsat_dv',rho);                               % Find saturation temperature from density
            hsat = Water_h(rho, Tsat);                                      % Find saturated vapor enthalpy at this density
            Xmin = H2Oprop('x_dt',rho,Tt);                               % Find quality for this density at the triple point
            hmin = Xmin*hvt + (1-Xmin)*hlt;                                 % Find enthalpy for two-phase mixture at triple point
            if (h >= hsat)                                                  % Vapor
                T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]);               % Use root-finding method to compute temperature
            elseif (h >= hmin)                                              % Two-Phase with Liquid
                T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tt Tsat]);        % Use root-finding method to compute temperature
            else                                                            % Multi-Phase with Ice
                T = NaN;                                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rholt)                                            % Liquid or Two-Phase with Vapor (No Density Anomaly)
            Tsat = H2Oprop('tsat_dl',rho);                               % Find saturation temperature from density
            hsat = Water_h(rho, Tsat);                                      % Find saturated liquid enthalpy at this density
            Xmin = H2Oprop('x_dt',rho,Tt);                               % Find quality for this density at the triple point
            hmin = Xmin*hvt + (1-Xmin)*hlt;                                 % Find enthalpy for two-phase mixture at triple point
            if (h >= hsat)                                                  % Liquid
                T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]);               % Use root-finding method to compute temperature
            elseif (h >= hmin)                                              % Two-Phase with Vapor
                T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tt Tsat]);        % Use root-finding method to compute temperature
            else                                                            % Ice
                T = NaN;                                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoamax)                                          % Liquid or Two-Phase with Vapor (Density Anomaly)
            if h >= ha                                                      % Warm Liquid or Two-Phase with Vapor
                Tsat = H2Oprop('tsat_dl',rho);                               % Find saturation temperature from density
                hsat = Water_h(rho, Tsat);                                      % Find saturated liquid enthalpy at this density
                if (h >= hsat)                                                  % Liquid
                    T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]);               % Use root-finding method to compute temperature
                else                                                            % Two-Phase with Vapor
                    T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tt Tsat]);        % Use root-finding method to compute temperature
                end
            else                                                        % Cold Liquid or Two-Phase with Vapor
                Tsat = fzero(@(x) GoodSatProp(x,'rholsat') - rho, [Tt Trhomax]);% Find saturation temperature from density
                hsat = Water_h(rho, Tsat);                                      % Find saturated liquid enthalpy at this density
                Tm   = fzero(@(x) Water_P(rho, x) - MeltPress(x), [Tmin Tt]);   % Find melting temperature at this density
                hm   = Water_h(rho, Tm);                                        % Find enthalpy of liquid at melting point with this density
                if (h >= hsat)                                                  % Two-Phase with Vapor
                    T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tsat Tmax]);      % Use root-finding method to compute temperature
                elseif (h >= hm)                                                % Liquid
                    T = fzero(@(x) Water_h(rho, x) - h, [Tm Tsat]);                 % Use root-finding method to compute temperature
                else                                                            % Ice
                    T = NaN;                                                        % Would need the ice correlation to compute this...
                end
            end
        elseif (rho < rholTmin)                                         % Pressurized Liquid (up to density of liquid at liquid-Ice Ih-Ice III triple point)
            Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x),x) - rho, [Tmin Tt]);% Find melting temperature from density
            hm = Water_h(rho,Tm);                                           % Find enthalpy of liquid at melting point with this density
            if (h >= hm)                                                    % Liquid
                T = fzero(@(x) Water_h(rho, x) - h, [Tm Tmax]);                 % Use root-finding method to compute temperature
            else                                                            % Ice
                T = NaN;                                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoicemax)                                        % Very Pressurized Liquid (Over 200 MPa)
            Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x,'hp'),x) - rho, ...
                 [Tmin Ticemax]);                                           % Find melting temperature from density
            hm = Water_h(rho,Tm);                                           % Find enthalpy of liquid at melting point with this density
            if (h >= hm)                                                    % Liquid
                T = fzero(@(x) Water_h(rho, x) - h, [Tm Tmax]);                 % Use root-finding method to compute temperature
            else                                                            % High-Pressure Ice
                T = NaN;                                                        % There's no correlation for this stuff...
            end
        else                                                            % Too Dense (High-Pressure Ice or Over 20 GPa)
            T = NaN;                                                        % There's no correlation for this stuff...
        end
        OutVar = T;         % Output Temperature
        
%---------------------------------------------------------------------------------------   
        
    case 'd_ph'         % Density from pressure, enthalpy
        P = Var1;           % Input pressure
        h = Var2;           % Input specific enthalpy
        if (P < Pmin)                           % Too rarefied
            rho = NaN;                              % Pressures this low likely violate continuum approximations...
        elseif (P < Pt)                         % Thin vapor (below triple point pressure)
            Tsat = MeltTemp(P);                     % Compute the sublimation temperature at this pressure
            rhosat = H2Oprop('d_pt',P,Tsat);     % Compute the density of saturated vapor at this pressure
            hsat = Water_h(rhosat, Tsat);           % Compute the enthalpy of saturated vapor at this pressure
            if (h >= hsat)                              % Vapor
                T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            else                                        % Ice
                rho = NaN;                                  % Would need the ice correlation to compute this...
            end
        elseif (P < Pc)                         % Liquid or vapor
            Tm = MeltTemp(P);                       % Compute the melting temperature at this pressure
            rhom = H2Oprop('d_pt',P,Tm);         % Compute the density of freezing liquid at this pressure
            hm = Water_h(rhom,Tm);                  % Compute the enthalpy of freezing liquid at this pressure
            Tsat = H2Oprop('tsat_p',P);          % Compute the saturation temperature at this pressure
            [junk rholsat rhovsat] = GoodSatProp(Tsat); % Compute the saturated liquid and vapor densities at this pressure
            hlsat = Water_h(rholsat,Tsat);          % Compute the saturated liquid enthalpy at this pressure
            hvsat = Water_h(rhovsat,Tsat);          % Compute the saturated vapor enthalpy at this pressure
            if (h <= hm)                                % Ice
                rho = NaN;                                  % Would need the ice correlation to compute this...
            elseif (h <= hlsat)                         % Liquid
                T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tm Tsat]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            elseif (h < hvsat)                          % Two-Phase
                rho = fzero(@(x) H2Oprop('h_dt',x,Tsat) - h, [rhovsat rholsat]); % Use root-finding method to compute density
            else                                        % Vapor
                T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            end
        elseif (P < Picemax)                    % Supercritical fluid
            Tm = MeltTemp(P);                       % Compute the melting temperature at this pressure
            rhom = H2Oprop('d_pt',P,Tm);         % Compute the density of freezing fluid at this pressure
            hm = Water_h(rhom,Tm);                  % Compute the enthalpy of freezing fluid at this pressure
            if (h >= hm)                                % Fluid
                T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tm Tmax]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            else                                        % Ice
                rho = NaN;                                  % There's no correlation for this stuff...
            end
        else                                    % Too pressurized (cannot verify fluidity)
            rho = NaN;                              % There's no correlation for this stuff...
        end
        OutVar = rho;       % Output density
        
%---------------------------------------------------------------------------------------   
        
    case 't_ph'         % Temperature from pressure, enthalpy 
        P = Var1;           % Input pressure
        h = Var2;           % Input specific enthalpy
        rho = H2Oprop('d_ph',P,h); % Compute density from given pressure and enthalpy
        if isnan(rho)           % Do the given pressure and enthalpy lie in a valid region?
            T = NaN;                % If not, propagate the NaN
        else                    % If so, we can compute temperature
            T = H2Oprop('t_dh',rho,h); % Compute temperature from density and enthalpy
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   

    case 't_du'         % Temperature from density, energy
        rho = Var1;         % Input density
        u   = Var2;         % Input specific internal energy
        if (rho < rhomin)                           % Too rarefied
            T = NaN;                                    % Densities this low likely violate continuum approximations...
        elseif (rho < rhosatmin)                    % Extremely Thin Vapor (less dense than 50 K saturated vapor)
            if (u >= uvsatmin)                              % Vapor
                T = fzero(@(x) Water_u(rho, x) - u, [Tabsmin Tmax]); % Use root-finding method to compute temperature
            else                                            % Non-Physical
                T = NaN;                                        % It's not ice at this density...
            end
        elseif (rho < rhovt)                        % Thin Vapor (less dense than triple point vapor)
            Tsat = H2Oprop('tsat_dv',rho);           % Find saturation temperature from density
            usat = Water_u(rho, Tsat);                  % Find saturated vapor energy at this density
            if (u >= usat)                                  % Vapor
                T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature
            else                                            % Ice
                T = NaN;                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoc)                         % Vapor or Two-Phase with Liquid
            Tsat = H2Oprop('tsat_dv',rho);           % Find saturation temperature from density
            usat = Water_u(rho, Tsat);                  % Find saturated vapor energy at this density
            Xmin = H2Oprop('x_dt',rho,Tt);           % Find quality for this density at the triple point
            umin = Xmin*uvt + (1-Xmin)*ult;             % Find energy for two-phase mixture at triple point
            if (u >= usat)                                  % Vapor
                T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature
            elseif (u >= umin)                              % Two-Phase with Liquid
                T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tt Tsat]); % Use root-finding method to compute temperature
            else                                            % Multi-Phase with Ice
                T = NaN;                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rholt)                        % Liquid or Two-Phase with Vapor (No Density Anomaly)
            Tsat = H2Oprop('tsat_dl',rho);           % Find saturation temperature from density
            usat = Water_u(rho, Tsat);                  % Find saturated liquid energy at this density
            Xmin = H2Oprop('x_dt',rho,Tt);           % Find quality for this density at the triple point
            umin = Xmin*uvt + (1-Xmin)*ult;             % Find energy for two-phase mixture at triple point
            if (u >= usat)                                  % Liquid
                T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature
            elseif (u >= umin)                              % Two-Phase with Vapor
                T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tt Tsat]); % Use root-finding method to compute temperature
            else                                            % Ice
                T = NaN;                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoamax)                      % Liquid or Two-Phase with Vapor (Density Anomaly)
            if u >= ua                                      % Warm Liquid or Two-Phase with Vapor
                Tsat = H2Oprop('tsat_dl',rho);               % Find saturation temperature from density (assume above Trhomax)
                usat = Water_u(rho, Tsat);                      % Find saturated liquid energy at this density above Trhomax
                if (u >= usat)                                      % Liquid
                    T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]);   % Use root-finding method to compute temperature
                else                                                % Two-Phase with Vapor
                    T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tt Tsat]); % Use root-finding method to compute temperature
                end
            else                                    % Cold Liquid or Two-Phase with Vapor
                Tsat = fzero(@(x) GoodSatProp(x,'rholsat') - rho, [Tt Trhomax]); % Find saturation temperature from density (assume below Trhomax)
                usat = Water_u(rho, Tsat);                      % Find saturated liquid energy at this density below Trhomax
                Tm   = fzero(@(x) Water_P(rho, x) - MeltPress(x), [Tmin Tt]); % Find melting temperature at this density
                um   = Water_u(rho, Tm);                        % Find energy of liquid at melting point with this density
                if (u >= usat)                                      % Two-Phase with Vapor
                    T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature
                elseif (u >= um)                                    % Liquid
                    T = fzero(@(x) Water_u(rho, x) - u, [Tm Tsat]);     % Use root-finding method to compute temperature
                else                                                % Ice
                    T = NaN;                                            % Would need the ice correlation to compute this...
                end
            end
        elseif (rho < rholTmin)                     % Pressurized Liquid (up to density of liquid at liquid-Ice Ih-Ice III triple point)
            Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x),x) - rho, [Tmin Tt]); % Find melting temperature from density
            um = Water_u(rho,Tm);                       % Find energy of liquid at melting point with this density
            if (u >= um)                                    % Liquid
                T = fzero(@(x) Water_u(rho, x) - u, [Tm Tmax]); % Use root-finding method to compute temperature
            else                                            % Ice
                T = NaN;                                        % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoicemax)                    % Very Pressurized Liquid (over 200 MPa)
            Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x,'hp'),x) - rho, ...
                 [Tmin Ticemax]);                       % Find melting temperature from density
            um = Water_u(rho,Tm);                       % Find energy of liquid at melting point with this density
            if (u >= um)                                    % Liquid
                T = fzero(@(x) Water_u(rho, x) - u, [Tm Tmax]); % Use root-finding method to compute temperature
            else                                            % High-Pressure Ice
                T = NaN;                                        % There's no correlation for this stuff...
            end
        else                                        % Too Dense (high-pressure ice or over 20 GPa)
            T = NaN;                                    % There's no correlation for this stuff, and we are well outside the region of validity!
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   
        
    case 'd_pu'         % Density from pressure, energy
        P = Var1;               % Input pressure
        u = Var2;               % Input specific internal energy
        if (P < Pmin)                       % Too rarefied
            rho = NaN;                          % Pressures this low likely violate continuum approximations...
        elseif (P < Pt)                     % Thin vapor (below triple point pressure)
            Tsat = MeltTemp(P);                 % Compute the sublimation temperature at this pressure
            rhosat = H2Oprop('d_pt',P,Tsat); % Compute the density of saturated vapor at this pressure
            usat = Water_u(rhosat, Tsat);       % Compute the energy of saturated vapor at this pressure
            if (u >= usat)                          % Vapor
                T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            else                                    % Ice
                rho = NaN;                              % Would need the ice correlation to compute this...
            end
        elseif (P < Pc)                     % Liquid, vapor, or two-phase
            Tm = MeltTemp(P);                   % Compute the melting temperature at this pressure
            rhom = H2Oprop('d_pt',P,Tm);     % Compute the density of freezing liquid at this pressure
            um = Water_u(rhom,Tm);              % Compute the energy of freezing liquid at this pressure
            Tsat = H2Oprop('tsat_p',P);      % Compute the saturation temperature at this pressure
            [junk rholsat rhovsat] = GoodSatProp(Tsat); % Compute the saturated liquid and vapor densities at this pressure
            ulsat = Water_u(rholsat,Tsat);      % Compute the saturated liquid energy at this pressure
            uvsat = Water_u(rhovsat,Tsat);      % Compute the saturated vapor energy at this pressure
            if (u <= um)                            % Ice
                rho = NaN;                              % Would need the ice correlation to compute this...
            elseif (u <= ulsat)                     % Liquid
                T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tm Tsat]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            elseif (u < uvsat)                      % Two-Phase
                rho = fzero(@(x) H2Oprop('u_dt',x,Tsat) - u, [rhovsat rholsat]); % Use root-finding method to compute density
            else                                    % Vapor
                T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);               % Use the temperature and given pressure to find density
            end
        elseif (P < Picemax)                % Supercritical fluid
            Tm = MeltTemp(P);                   % Compute the melting temperature at this pressure
            rhom = H2Oprop('d_pt',P,Tm);     % Compute the density of freezing fluid at this pressure
            um = Water_u(rhom,Tm);              % Compute the energy of freezing fluid at this pressure
            if (u >= um)                            % Fluid
                T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tm Tmax]); % Use root-finding method to compute temperature
                rho = H2Oprop('d_pt',P,T);           % Use the temperature and given pressure to find density
            else                                    % Ice
                rho = NaN;                              % There's no correlation for this stuff...
            end
        else                                % Too pressurized (cannot verify fluidity)
            rho = NaN;                          % There's no correlation for this stuff...
        end
        OutVar = rho;           % Output density
        
%---------------------------------------------------------------------------------------   
        
    case 't_pu'         % Temperature from pressure, energy
        P = Var1;           % Input pressure
        u = Var2;           % Input specific internal energy
        rho = H2Oprop('d_pu',P,u); % Find density from pressure and energy
        if isnan(rho)           % Do the given pressure and energy lie in a valid region?
            T = NaN;                % If not, propagate the NaN
        else                    % If so, we can compute temperature
            T = H2Oprop('t_du',rho,u); % Compute temperature from density and enthalpy
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   

    case 'd_ts'         % Density from temperature, entropy
        T = Var1;           % Input temperature
        s = Var2;           % Input specific entropy
        if (T < Tmin)               % Cold Vapor
            Psat = MeltPress(T);        % Find sublimation pressure from temperature
            rhosat = H2Oprop('d_pt',Psat,T); % Find vapor density at sublimation point
            ssat = Water_s(rhosat,T);   % Find vapor entropy at sublimation point
            if (s >= ssat*.999)             % Vapor
                rho = fzero(@(x) s - Water_s(x,T),[rhomin*.999 rhosat*1.001]); % Use root-finding method to compute density
            else                            % Ice
                rho = NaN;                      % Would need the ice correlation to compute this...
            end
        elseif (T < Tt)             % Cold Vapor or Pressurized Liquid
            [Pm Ps Pf] = MeltPress(T);  % Find all phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find vapor density at sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find liquid density at Ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find liquid density at high-pressure ice freezing point
            ss = Water_s(rhos,T);       % Find vapor entropy at sublimation point
            sm = Water_s(rhom,T);       % Find liquid entropy at Ice Ih melting point
            sf = Water_s(rhof,T);       % Find liquid entropy at high-pressure ice freezing point
            if (s < sf)                     % High-pressure Ice
                rho = NaN;                      % There's no correlation for this stuff...
            elseif (s <= sm)                % Pressurized Liquid
                rho = fzero(@(x) s - Water_s(x,T),[rhom rhof]); % Use root-finding method to compute density
            elseif (s < ss)                 % Ice Ih
                rho = NaN;                      % Would need the ice correlation to compute this...
            else                            % Vapor
                rho = fzero(@(x) s - Water_s(x,T),[rhomin rhos]); % Use root-finding method to compute density
            end
        elseif (T < Tc)             % Vapor, Liquid, or Two-Phase
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find saturation properties at this temperature
            sl = Water_s(rholsat,T);    % Find entropy of saturated liquid
            sv = Water_s(rhovsat,T);    % Find entropy of saturated vapor
            Pm = MeltPress(T);          % Find melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find density of liquid at melting point
            sm = Water_s(rhom,T);       % Find entropy of freezing liquid
            if (s < sm*.999)                % High-Pressure Ice
                rho = NaN;                      % There's no correlation for this stuff...
            elseif (s <= sl)                % Liquid
                rho = fzero(@(x) s - Water_s(x,T),[rholsat rhom*1.001]); % Use root-finding method to compute density
            elseif (s < sv)                 % Two-Phase Mixture
                rho = fzero(@(x) s - H2Oprop('s_dt',x,T),[rhovsat rholsat]); % Use root-finding method to compute density
            else                            % Vapor
                rho = fzero(@(x) s - Water_s(x,T),[rhomin rhovsat]); % Use root-finding method to compute density
            end
        elseif (T < Ticemax)        % Supercritical Fluid
            Pm = MeltPress(T);          % Find melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find density of liquid at melting point
            sm = Water_s(rhom,T);       % Find entropy of freezing liquid
            if (s < (sm - 10*eps(sm)))      % High-Pressure Ice
                rho = NaN;                      % There's no correlation for this stuff...
            else                            % Supercritical Fluid
                rho = fzero(@(x) s - Water_s(x,T),[rhomin rhom+5*eps(rhom)]); % Use root-finding method to compute density
            end
        else                        % Hot Supercritical Fluid
            try                         % Use error handling in case the true density is too large
                rho = fzero(@(x) s - Water_s(x,T),[rhomin rhoicemax]); % Use root-finding method to compute density
            catch                       % If that didn't work...
                [errmsg,errcode] = lasterr; % Fetch the error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error was the expected failure of fzero...
                    rho = NaN;                      % ...fail gracefully with a NaN
                else
                    rethrow(lasterror);             % Otherwise, pass the error through
                end
            end
        end
        OutVar = rho;       % Output density
        
%---------------------------------------------------------------------------------------   

    case 't_ds'         % Temperature from density, entropy
        rho = Var1;         % Input density
        s   = Var2;         % Input specific entropy
        if (rho < rhomin)               % Too rarefied
            T = NaN;                        % Pressures this low likely violate continuum approximations... 
        elseif (rho < rhosatmin)        % Extremely Thin Vapor (less dense than 50 K saturated vapor)
            if (s >= svsatmin)                  % Vapor
                T = fzero(@(x) Water_s(rho, x) - s, [Tabsmin Tmax]); % Use root-finding method to compute temperature
            else                                % Non-Physical
                T = NaN;                            % It's not ice at this density...
            end
        elseif (rho < rhovt)            % Thin Vapor (less dense than triple point vapor)
            Tsat = H2Oprop('tsat_dv',rho); % Find sublimation temperature from density
            ssat = Water_s(rho, Tsat);      % Find saturated vapor entropy
            if (s >= ssat)                      % Vapor
                T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature
            else                                % Ice
                T = NaN;                            % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoc)             % Vapor or Two-Phase with Liquid
            Tsat = H2Oprop('tsat_dv',rho); % Find saturation temperature from density
            ssat = Water_s(rho, Tsat);      % Find saturated vapor entropy
            Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point
            smin = Xmin*svt + (1-Xmin)*slt; % Find entropy for two-phase mixture at the triple point
            if (s >= ssat)                      % Vapor
                T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature
            elseif (s >= smin)                  % Two-Phase with Liquid
                T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tt Tsat]); % Use root-finding method to compute temperature
            else                                % Ice
                T = NaN;                            % Would need the ice correlation to compute this...
            end
        elseif (rho < rholt)            % Liquid or Two-Phase with Vapor (No Density Anomaly)
            Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density
            ssat = Water_s(rho, Tsat);      % Find saturated liquid entropy
            Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point
            smin = Xmin*svt + (1-Xmin)*slt; % Find entropy for two-phase mixture at the triple point
            if (s >= ssat)                      % Liquid
                T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature
            elseif (s >= smin)                  % Two-Phase
                T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tt Tsat]); % Use root-finding method to compute temperature
            else                                % Ice
                T = NaN;                            % Would need the ice correlation to compute this...
            end
        elseif (rho < rhoamax)          % Liquid or Two-Phase with Vapor (Density Anomaly)
            if s >= sa                          % Warm Liquid or Two-Phase with Vapor
                Tsat = H2Oprop('tsat_dl',rho);   % Find saturation temperature from density (assume above Trhomax)
                ssat = Water_s(rho, Tsat);          % Find saturated liquid entropy
                if (s >= ssat)                          % Liquid
                    T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature
                else                                    % Two-Phase with Vapor
                    T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tt Tsat]); % Use root-finding method to compute temperature
                end
            else                                % Cold Liquid or Two-Phase Mixture
                Tsat = fzero(@(x) GoodSatProp(x,'rholsat') - rho, [Tt Trhomax]); % Find saturation temperature from density (assume below Trhomax)
                ssat = Water_s(rho, Tsat);          % Find saturated liquid entropy
                Tm   = fzero(@(x) Water_P(rho, x) - MeltPress(x), [Tmin Tt]); % Find melting temperature at this density
                sm   = Water_s(rho, Tm);            % Find entropy of liquid at freezing point
                if (s >= ssat)                          % Two-Phase with Vapor
                    T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature
                elseif (s >= sm)                        % Liquid
                    T = fzero(@(x) Water_s(rho, x) - s, [Tm Tsat]); % Use root-finding method to compute temperature
                else                                    % Ice
                    T = NaN;                                % Would need the ice correlation to compute this...
                end
            end
        elseif (rho < rholTmin)                 % Pressurized Liquid (up to density of liquid at liquid-Ice Ih-Ice III triple point)
            Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x),x) - rho, [Tmin Tt]); % Find melting temperature from density
            sm = Water_s(rho,Tm);                   % Find entropy of liquid at melting point
            if (s >= sm)                                % Liquid
                T = fzero(@(x) Water_s(rho, x) - s, [Tm Tmax]); % Use root-finding method to compute temperature
            else                                        % Ice Ih
                T = NaN;                                    % Would need the ice correlation to compute this...
            end
        elseif (rho <= rhoicemax)               % Very Pressurized Liquid (over 200 MPa)
            Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x,'hp'),x) - rho, ...
                 [Tmin Ticemax]);                   % Find melting temperature from density
            sm = Water_s(rho,Tm);                   % Find entropy of liquid at melting point
            if (s >= sm)                                % Liquid
                T = fzero(@(x) Water_s(rho, x) - s, [Tm Tmax]); % Use root-finding method to compute temperature
            else                                        % High-Pressure Ice
                T = NaN;                                    % There's no correlation for this stuff...
            end
        else                                    % Too Dense (high-pressure ice or over 20 GPa)
            T = NaN;                                % There's no correlation for this stuff, and we are well outside the region of validity!
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   

    case 'p_du'         % Pressure from density, energy
        rho = Var1;         % Input density
        u = Var2;           % Input specific internal energy
        T = H2Oprop('t_du',rho,u); % Calculate temperature from density and energy
        if isnan(T)             % Was the temperature calculation successful?
            P = NaN;                % If not, propagate the NaN
        else
            P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure
        end
        OutVar = P;         % Output pressure

%---------------------------------------------------------------------------------------   

    case 'p_dh'         % Pressure from density, enthalpy
        rho = Var1;         % Input density
        h = Var2;           % Input specific enthalpy
        T = H2Oprop('t_dh',rho,h); % Calculate temperature from density and enthalpy
        if isnan(T)             % Was the temperature calculation successful?
            P = NaN;                % If not, propagate the NaN
        else
            P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure
        end
        OutVar = P;         % Output pressure

%---------------------------------------------------------------------------------------   

    case 'p_ds'         % Pressure from density, entropy
        rho = Var1;         % Input density
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_ds',rho,s); % Calculate temperature from density and entropy
        if isnan(T)             % Was the temperature calculation successful?
            P = NaN;                % If not, propagate the NaN
        else
            P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure
        end
        OutVar = P;         % Output pressure

%---------------------------------------------------------------------------------------   

    case 'p_ts'         % Pressure from temperature, entropy
        T = Var1;           % Input temperature
        s = Var2;           % Input specific entropy
        rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy
        if isnan(rho)           % Was the density calculation successful?
            P = NaN;                % If not, propagate the NaN
        else
            P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure
        end
        OutVar = P;         % Output pressure
%---------------------------------------------------------------------------------------   

    case 'u_dh'         % Specific energy from density, enthalpy
        rho = Var1;         % Input density
        h = Var2;           % Input specific enthalpy
        T = H2Oprop('t_dh',rho,h); % Calculate temperature from density and enthalpy
        if isnan(T)             % Was the temperature calculation successful?
            u = NaN;                % If not, propagate the NaN
        else
            u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy
        end
        OutVar = u;         % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 'u_ds'         % Specific energy from density, entropy
        rho = Var1;         % Input density
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_ds',rho,s); % Calculate temperature from density and entropy
        if isnan(T)             % Was the temperature calculation successful?
            u = NaN;                % If not, propagate the NaN
        else
            u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy
        end
        OutVar = u;         % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 'u_pt'         % Specific energy from pressure, temperature
        P = Var1;           % Input pressure
        T = Var2;           % Input temperature
        rho = H2Oprop('d_pt',P,T); % Calculate density from pressure and temperature
        if isnan(rho)           % Was the density calculation successful?
            u = NaN;                % If not, propagate the NaN
        else
            u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy
        end
        OutVar = u;         % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 'u_ph'         % Specific energy from pressure, enthalpy
        P = Var1;           % Input pressure
        h = Var2;           % Input specific enthalpy
        T = H2Oprop('t_ph',P,h); % Calculate temperature from pressure and enthalpy
        rho = H2Oprop('d_ph',P,h); % Calculate density from pressure and enthalpy
        if isnan(T)||isnan(rho) % Were the temperature and density calculations successful?
            u = NaN;                % If not, propagate the NaN
        else
            u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy
        end
        OutVar = u;         % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 'u_ts'         % Specific energy from entropy, temperature
        T = Var1;           % Input temperature
        s = Var2;           % Input specific entropy
        rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy
        if isnan(rho)           % Was the density calculation successful?
            u = NaN;                % If not, propagate the NaN
        else
            u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy
        end
        OutVar = u;         % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 'h_du'         % Specific enthalpy from density, energy
        rho = Var1;         % Input density
        u = Var2;           % Input specific internal energy
        T = H2Oprop('t_du',rho,u); % Calculate temperature from density and internal energy
        if isnan(T)             % Was the temperature calculation successful?
            h = NaN;                % If not, propagate the NaN
        else
            h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 'h_ds'         % Specific enthalpy from density, entropy
        rho = Var1;         % Input density
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_ds',rho,s); % Calculate temperature from density and entropy
        if isnan(T)             % Was the temperature calculation successful?
            h = NaN;                % If not, propagate the NaN
        else
            h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 'h_pt'         % Specific enthalpy from pressure, temperature
        P = Var1;           % Input pressure
        T = Var2;           % Input temperature
        rho = H2Oprop('d_pt',P,T); % Calculate density from pressure and temperature
        if isnan(rho)           % Was the density calculation successful?
            h = NaN;                % If not, propagate the NaN
        else
            h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 'h_pu'         % Specific enthalpy from pressure, energy
        P = Var1;           % Input pressure
        u = Var2;           % Input specific internal energy
        T = H2Oprop('t_pu',P,u); % Calculate temperature from pressure and energy
        rho = H2Oprop('d_pu',P,u); % Calculate density from pressure and energy
        if isnan(T)||isnan(rho) % Were the density and temperature calculations successful?
            h = NaN;                % If not, propagate the NaN
        else
            h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 'h_ts'         % Specific enthalpy from entropy, temperature
        T = Var1;           % Input temperature
        s = Var2;           % Input specific entropy
        rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy
        if isnan(rho)           % Was the density calculation successful?
            h = NaN;                % If not, propagate the NaN
        else
            h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 's_dh'         % Specific entropy from density, enthalpy
        rho = Var1;         % Input density
        h = Var2;           % Input specific enthalpy
        T = H2Oprop('t_dh',rho,h); % Calculate temperature from density and enthalpy
        if isnan(T)             % Was the temperature calculation successful?
            s = NaN;                % If not, propagate the NaN
        else
            s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy
        end
        OutVar = s;         % Output specific entropy

%---------------------------------------------------------------------------------------   

    case 's_du'         % Specific entropy from density, energy
        rho = Var1;         % Input density
        u = Var2;           % Input specific internal energy
        T = H2Oprop('t_du',rho,u); % Calculate temperature from density and energy
        if isnan(T)             % Was the temperature calculation successful?
            s = NaN;                % If not, propagate the NaN
        else
            s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy
        end
        OutVar = s;         % Output specific entropy

%---------------------------------------------------------------------------------------   

    case 's_pt'         % Specific entropy from pressure, temperature
        P = Var1;           % Input pressure
        T = Var2;           % Input temperature
        rho = H2Oprop('d_pt',P,T); % Calculate density from pressure and temperature
        if isnan(rho)           % Was the density calculation successful?
            s = NaN;                % If not, propagate the NaN
        else
            s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy
        end
        OutVar = s;         % Output specific entropy

%---------------------------------------------------------------------------------------   

    case 's_ph'         % Specific entropy from pressure, enthalpy
        P = Var1;           % Input pressure
        h = Var2;           % Input specific enthalpy
        T = H2Oprop('t_ph',P,h); % Calculate temperature from pressure and enthalpy
        rho = H2Oprop('d_ph',P,h); % Calculate density from pressure and enthalpy
        if isnan(T)||isnan(rho) % Were the density and temperature calculations successful?
            s = NaN;                % If not, propagate the NaN
        else
            s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy
        end
        OutVar = s;         % Output specific entropy

%---------------------------------------------------------------------------------------   

    case 's_pu'         % Specific entropy from pressure, energy
        P = Var1;           % Input pressure
        u = Var2;           % Input specific internal energy
        T = H2Oprop('t_pu',P,u); % Calculate temperature from pressure and energy
        rho = H2Oprop('d_pu',P,u); % Calculate density from pressure and energy
        if isnan(T)||isnan(rho) % Were the density and temperature calculations successful?
            s = NaN;                % If not, propagate the NaN
        else
            s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy
        end
        OutVar = s;         % Output specific entropy

%---------------------------------------------------------------------------------------   

    case 't_ps'         % Temperature from pressure, entropy
        P = Var1;               % Input pressure
        s = Var2;               % Input specific entropy
        if (P < Pmin)               % Too rarefied
            T = NaN;                    % Pressures this low likely violate continuum approximations...
        elseif (P < Pt)             % Thin vapor (below triple point pressure)
            Tsat = MeltTemp(P);         % Find sublimation temperature at this pressure
            rhosat = H2Oprop('d_pt',P,Tsat); % Find vapor density at sublimation point
            ssat = Water_s(rhosat, Tsat); % Find vapor entropy at sublimation point
            if (s >= ssat)                  % Vapor
                T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tsat Tmax]); % Use root-finding method to find temperature
            else                            % Ice
                T = NaN;                        % Would need the ice correlation to calculate this...
            end
        elseif (P < Pc)             % Liquid or vapor
            Tm = MeltTemp(P);           % Find Ice Ih melting point at this pressure
            rhom = H2Oprop('d_pt',P,Tm); % Find liquid density at melting point
            sm = Water_s(rhom,Tm);      % Find liquid entropy at melting point
            Tsat = H2Oprop('tsat_p',P); % Find saturation temperature at this pressure
            [junk rholsat rhovsat] = GoodSatProp(Tsat); % Calculate saturated liquid and vapor densities at this pressure
            slsat = Water_s(rholsat,Tsat); % Find entropy of saturated liquid
            svsat = Water_s(rhovsat,Tsat); % Find entropy of saturated vapor
            if (s <= sm)                    % Ice
                T = NaN;                        % Would need the ice correlation to calculate this...
            elseif (s <= slsat)             % Liquid
                T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tm Tsat]); % Use root-finding method to find temperature
            elseif (s < svsat)              % Two-Phase Mixture
                T = Tsat;                       % Return the saturation temperature
            else                            % Vapor
                T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tsat Tmax]); % Use root-finding method to find temperature
            end
        elseif (P < Picemax)        % Supercritical Fluid
            Tm = MeltTemp(P);           % Find the melting point at this pressure
            rhom = H2Oprop('d_pt',P,Tm); % Find liquid density at the melting point
            sm = Water_s(rhom,Tm);      % Find liquid entropy at the melting point
            if (s >= sm)                    % Fluid
                T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tm Tmax]); % Use root-finding method to find temperature
            else                            % Ice
                T = NaN;                        % Would need the ice correlation to calculate this...
            end
        else                        % Too pressurized (cannot verify fluidity)
            T = NaN;                    % There's no correlation for this stuff...
        end
        OutVar = T;             % Output temperature
        
%---------------------------------------------------------------------------------------   

    case 'd_ps'         % Density from pressure, entropy
        P = Var1;           % Input pressure
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_ps',P,s); % Calculate temperature from pressure and entropy
        if isnan(T)             % Did the temperature calculation work?
            rho = NaN;              % If not, propagate the NaN
        else
            rho = H2Oprop('d_ts',T,s); % If so, use the temperature and entropy to find the density
        end
        OutVar = rho;       % Output density

%---------------------------------------------------------------------------------------   

    case 'u_ps'         % Specific energy from pressure, entropy
        P = Var1;           % Input pressure
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_ps',P,s); % Calculate temperature from pressure and entropy
        rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy
        if isnan(T)             % Did the temperature calculation work?
            u = NaN;                % If not, propagate the NaN
        else
            u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find the internal energy
        end
        OutVar = u;         % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 'h_ps'         % Specific enthalpy from pressure, entropy
        P = Var1;           % Input pressure
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_ps',P,s); % Calculate temperature from pressure and entropy
        rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy
        if isnan(T)             % Did the temperature calculation work?
            h = NaN;                % If not, propagate the NaN
        else
            h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find the enthalpy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 't_hs'         % Temperature from enthalpy, entropy
        
    % Note that there are local maxima and minima in entropy along the high-pressure melting curve.
    % These are outside the official pressure range, and likely represent
    % model inaccuaracies rather than a physical phenomenon.  This function
    % does some gymnastics around these local extrema for numerical
    % reasons, but one should be cautious in interpreting results for fluid
    % states above 8.9 GPa and 600 K.
        
        h = Var1;           % Input specific enthalpy
        s = Var2;           % Input specific entropy
        if (s < spimin)         % Ice
            T = NaN;                % Would need the ice correlation to compute this...
        elseif (s < slt)        % Ice or Liquid
            Tm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x),x), ...
                       [Tmin+100*eps(Tmin) Tt]); % Use root-finding method to find Ice Ih melting point from entropy
            Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                       [Tmin+100*eps(Tmin) Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy
            try
                T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tm Tf]); % Use root-finding method to find temperature from enthalpy and entropy
            catch               % The root-finding would fail if enthalpy were outside the physical range
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                    T = NaN;                % Ice or Non-Physical, so return NaN
                else                    % Other kinds of errors indicate real problems...
                    rethrow(lasterror);     % ...so kick the error upstairs
                end
            end
        elseif (s < spitmin)    % Ice, Liquid, or Two-Phase with Vapor (below anomalies near maximum of melting curve)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy
            hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid
            if (h < hl)                 % Two-Phase or Ice
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                           % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr;     % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                        % Liquid or High-Pressure Ice
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tf]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                           % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr;     % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= spitmax)   % Ice, Liquid, or Two-Phase with Vapor (between entropy minimum and entropy at top of melting curve)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy
            Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tsmax Tsmin]); % Use root-finding method to find Ice VII melting point between entropy maximum and entropy minimum from entropy
            Tfm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tsmin Ticemax]); % Use root-finding method to find Ice VII melting point above entropy minimum from entropy
            hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid
            hil = H2Oprop('h_pt',MeltPress(Tfl,'hp'),Tfl); % Find enthalpy of freezing liquid below entropy maximum
            hih = H2Oprop('h_pt',MeltPress(Tfh,'hp'),Tfh); % Find enthalpy of freezing liquid between entropy extrema
            if (h < hl)             % Two-Phase with Vapor or Ice
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            elseif (h < hil)        % Liquid (below anomalies)
                T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from enthalpy and entropy
            elseif (h < hih)        % Ice VII (between anomalies)
                T = NaN;                % There's no correlation for this stuff...
            else                    % Hot Liquid (above anomalies)
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= spimax)    % Ice, Liquid, or Two-Phase with Vapor (between entropy at top of melting curve and entropy maximum)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy
            Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tsmax Ticemax]); % Use root-finding method to find high-pressure ice melting point above entropy maximum from entropy
            hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid
            hil = H2Oprop('h_pt',MeltPress(Tfl,'hp'),Tfl); % Find enthalpy of freezing liquid below entropy maximum
            hih = H2Oprop('h_pt',MeltPress(Tfh,'hp'),Tfh); % Find enthalpy of freezing liquid between entropy extrema
            if (h < hl)             % Two-Phase with Vapor or Ice
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            elseif (h < hil)    % Liquid (below anomalies)
                T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from enthalpy and entropy
            elseif (h < hih)    % Ice VII (between anomalies)
                T = NaN;            % There's no correlation for this stuff...
            else                % Hot Liquid (above anomalies)
                Tfm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range, or if the state has density > rhoicemax
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    elseif isnan(Tfm)                   % If this h,s point is at a density above rhoicemax, Tfm would be NaN and fzero would also have failed
                        T = NaN;                            % Density too great, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= sdtmax)    % Liquid or Two-Phase with Vapor (below entropy at maximum density and temperature)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid
            if (h > hl)             % Liquid
                Tm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tm-10*eps(Tm)]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range, or if the state has density > rhoicemax
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    elseif isnan(Tm)                   % If this h,s point is at a density above rhoicemax, Tm would be NaN and fzero would also have failed
                        T = NaN;                            % Density too great, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                    % Two-Phase with Vapor
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= scrit)     % Liquid or Two-Phase with Vapor (above entropy at maximum density and temperature)
            try
                Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            catch                   % The root-finding would fail if the actual entropy were very close to the critical entropy
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc...
                    Tv = Tc;                    % Just set the boiling point to equal the critical point
                else                        % Other kinds of errors indicate real problems...
                    rethrow(lasterror);         % ...so kick the error upstairs
                end
            end
            hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid
            if (h > hl)                 % Liquid
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tmax*.999]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                        % Two-Phase with Vapor
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= svt)       % Vapor or Two-Phase with Liquid
            try
                Tv = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            catch                   % The root-finding would fail if the actual entropy were very close to the critical entropy
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc...
                    Tv = Tc;                    % Just set the boiling point to equal the critical point
                else                        % Other kinds of errors indicate real problems...
                    rethrow(lasterror);         % ...so kick the error upstairs
                end
            end
            hv = H2Oprop('hvsat_t',Tv); % Find entropy of saturated vapor
            if (h < hv)                 % Two-Phase with Liquid
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                        % Vapor
                try
                    T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv+1e5*eps(Tv) Tmax]); % Use root-finding method to find temperature from enthalpy and entropy
                catch                       % The root-finding would fail if enthalpy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= svsatmin)  % Vapor or Ice
            Ts = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tt]); % Use root-finding method to find sublimation point from entropy
            try
                T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Ts+1e5*eps(Ts) Tmax]); % Use root-finding method to find temperature from enthalpy and entropy
            catch                       % The root-finding would fail if enthalpy were outside the physical range
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                    T = NaN;                            % Ice or Non-Physical, so return NaN
                else                                % Other kinds of errors indicate real problems...
                    rethrow(lasterror);                 % ...so kick the error upstairs
                end
            end
        else                    % Only vapor allowed
            try
                T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tabsmin Tmax]); % Use root-finding method to find temperature from enthalpy and entropy
            catch                       % The root-finding would fail if enthalpy were outside the physical range
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                    T = NaN;                            % Ice or Non-Physical, so return NaN
                else                                % Other kinds of errors indicate real problems...
                    rethrow(lasterror);                 % ...so kick the error upstairs
                end
            end
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   

    case 'd_hs'         % Density from enthalpy, entropy
        h = Var1;           % Input specific enthalpy
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_hs',h,s); % Calculate temperature from enthalpy and entropy
        if isnan(T)             % If that didn't work...
            rho = NaN;              % ...propagate the NaN
        else
            rho = H2Oprop('d_ts',T,s); % Otherwise, calculate density from temperature and entropy
        end
        OutVar = rho;       % Output density

%---------------------------------------------------------------------------------------   

    case 'p_hs'         % Pressure from enthalpy, entropy
        h = Var1;           % Input specific enthalpy
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_hs',h,s); % Calculate temperature from enthalpy and entropy
        if isnan(T)             % If that didn't work...
            P = NaN;              % ...propagate the NaN
        else
            P = H2Oprop('p_ts',T,s); % Otherwise, calculate pressure from temperature and entropy
        end
        OutVar = P;       % Output density

%---------------------------------------------------------------------------------------   

    case 'u_hs'         % Specific Energy from enthalpy, entropy
        h = Var1;           % Input specific enthalpy
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_hs',h,s); % Calculate temperature from enthalpy and entropy
        if isnan(T)             % If that didn't work...
            u = NaN;              % ...propagate the NaN
        else
            u = H2Oprop('u_ts',T,s); % Otherwise, calculate specific internal energy from temperature and entropy
        end
        OutVar = u;       % Output specific internal energy

%---------------------------------------------------------------------------------------   

    case 't_us'         % Temperature from energy, entropy
        u = Var1;           % Input specific internal energy
        s = Var2;           % Input specific entropy
        if (s < spimin)         % Ice
            T = NaN;                % Would need the ice correlation to compute this...
        elseif (s < slt)        % Ice or Liquid
            Tm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x),x), ...
                       [Tmin+100*eps(Tmin) Tt]); % Use root-finding method to find Ice Ih melting point from entropy
            Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                       [Tmin+100*eps(Tmin) Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy
            try
                T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tm Tf]); % Use root-finding method to find temperature from energy and entropy
            catch               % The root-finding would fail if energy were outside the physical range
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                    T = NaN;                % Ice or Non-Physical, so return NaN
                else                    % Other kinds of errors indicate real problems...
                    rethrow(lasterror);     % ...so kick the error upstairs
                end
            end
        elseif (s < spitmin)    % Ice, Liquid, or Two-Phase with Vapor (below anomalies near maximum of melting curve)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy
            ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid
            if (u < ul)                 % Two-Phase or Ice
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy
                catch                           % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr;     % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                        % Liquid or High-Pressure Ice
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tf]); % Use root-finding method to find temperature from energy and entropy
                catch                           % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr;     % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= spitmax)   % Ice, Liquid, or Two-Phase with Vapor (between entropy minimum and entropy at top of melting curve)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy
            Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tsmax Tsmin]); % Use root-finding method to find Ice VII melting point between entropy maximum and entropy minimum from entropy
            Tfm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tsmin Ticemax]); % Use root-finding method to find Ice VII melting point above entropy minimum from entropy
            ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid
            uil = H2Oprop('u_pt',MeltPress(Tfl,'hp'),Tfl); % Find energy of freezing liquid below entropy maximum
            uih = H2Oprop('u_pt',MeltPress(Tfh,'hp'),Tfh); % Find energy of freezing liquid between entropy extrema
            if (u < ul)             % Two-Phase with Vapor or Ice
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            elseif (u < uil)        % Liquid (below anomalies)
                T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from energy and entropy
            elseif (u < uih)        % Ice VII (between anomalies)
                T = NaN;                % There's no correlation for this stuff...
            else                    % Hot Liquid (above anomalies)
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= spimax)    % Ice, Liquid, or Two-Phase with Vapor (between entropy at top of melting curve and entropy maximum)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy
            Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ...
                  [Tsmax Ticemax]); % Use root-finding method to find high-pressure ice melting point above entropy maximum from entropy
            ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid
            uil = H2Oprop('u_pt',MeltPress(Tfl,'hp'),Tfl); % Find energy of freezing liquid below entropy maximum
            uih = H2Oprop('u_pt',MeltPress(Tfh,'hp'),Tfh); % Find energy of freezing liquid between entropy extrema
            if (u < ul)             % Two-Phase with Vapor or Ice
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            elseif (u < uil)    % Liquid (below anomalies)
                T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from energy and entropy
            elseif (u < uih)    % Ice VII (between anomalies)
                T = NaN;            % There's no correlation for this stuff...
            else                % Hot Liquid (above anomalies)
                Tfm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range, or if the state has density > rhoicemax
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    elseif isnan(Tfm)                   % If this u,s point is at a density above rhoicemax, Tfm would be NaN and fzero would also have failed
                        T = NaN;                            % Density too great, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= sdtmax)    % Liquid or Two-Phase with Vapor (below entropy at maximum density and temperature)
            Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid
            if (u > ul)             % Liquid
                Tm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tm-10*eps(Tm)]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range, or if the state has density > rhoicemax
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    elseif isnan(Tm)                   % If this u,s point is at a density above rhoicemax, Tm would be NaN and fzero would also have failed
                        T = NaN;                            % Density too great, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                    % Two-Phase with Vapor
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= scrit)     % Liquid or Two-Phase with Vapor (above entropy at maximum density and temperature)
            try
                Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            catch                   % The root-finding would fail if the actual entropy were very close to the critical entropy
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc...
                    Tv = Tc;                    % Just set the boiling point to equal the critical point
                else                        % Other kinds of errors indicate real problems...
                    rethrow(lasterror);         % ...so kick the error upstairs
                end
            end
            ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid
            if (u > ul)                 % Liquid
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tmax*.999]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                        % Two-Phase with Vapor
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= svt)       % Vapor or Two-Phase with Liquid
            try
                Tv = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy
            catch                   % The root-finding would fail if the actual entropy were very close to the critical entropy
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc...
                    Tv = Tc;                    % Just set the boiling point to equal the critical point
                else                        % Other kinds of errors indicate real problems...
                    rethrow(lasterror);         % ...so kick the error upstairs
                end
            end
            uv = H2Oprop('uvsat_t',Tv); % Find energy of saturated vapor
            if (u < uv)                 % Two-Phase with Liquid
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            else                        % Vapor
                try
                    T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv+1e5*eps(Tv) Tmax]); % Use root-finding method to find temperature from energy and entropy
                catch                       % The root-finding would fail if energy were outside the physical range
                    [errmsg,errcode] = lasterr; % Retrieve error code
                    if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                        T = NaN;                            % Ice or Non-Physical, so return NaN
                    else                                % Other kinds of errors indicate real problems...
                        rethrow(lasterror);                 % ...so kick the error upstairs
                    end
                end
            end
        elseif (s <= svsatmin)  % Vapor or Ice
            Ts = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tt]); % Use root-finding method to find sublimation point from entropy
            try
                T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Ts+1e5*eps(Ts) Tmax]); % Use root-finding method to find temperature from energy and entropy
            catch                       % The root-finding would fail if energy were outside the physical range
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                    T = NaN;                            % Ice or Non-Physical, so return NaN
                else                                % Other kinds of errors indicate real problems...
                    rethrow(lasterror);                 % ...so kick the error upstairs
                end
            end
        else                    % Only vapor allowed
            try
                T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tabsmin Tmax]); % Use root-finding method to find temperature from energy and entropy
            catch                       % The root-finding would fail if energy were outside the physical range
                [errmsg,errcode] = lasterr; % Retrieve error code
                if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range...
                    T = NaN;                            % Ice or Non-Physical, so return NaN
                else                                % Other kinds of errors indicate real problems...
                    rethrow(lasterror);                 % ...so kick the error upstairs
                end
            end
        end
        OutVar = T;
        
%---------------------------------------------------------------------------------------   

    case 'd_us'         % Density from energy, entropy
        u = Var1;           % Input specific internal energy
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_us',u,s); % Calculate temperature from energy and entropy
        if isnan(T)             % If that didn't work...
            rho = NaN;              % ...propagate the NaN
        else
            rho = H2Oprop('d_ts',T,s); % Otherwise, calculate density from temperature and entropy
        end
        OutVar = rho;       % Output density

%---------------------------------------------------------------------------------------   

    case 'p_us'         % Pressure from energy, entropy
        u = Var1;           % Input specific internal energy
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_us',u,s); % Calculate temperature from energy and entropy
        if isnan(T)             % If that didn't work...
            P = NaN;                % ...propagate the NaN
        else
            P = H2Oprop('p_ts',T,s); % Otherwise, calculate pressure from temperature and entropy
        end
        OutVar = P;         % Output pressure

%---------------------------------------------------------------------------------------   

    case 'h_us'         % Specific Enthalpy from energy, entropy
        u = Var1;           % Input specific internal energy
        s = Var2;           % Input specific entropy
        T = H2Oprop('t_us',u,s); % Calculate temperature from energy and entropy
        if isnan(T)             % If that didn't work...
            h = NaN;                % ...propagate the NaN
        else
            h = H2Oprop('h_ts',T,s); % Otherwise, calculate enthalpy from temperature and entropy
        end
        OutVar = h;         % Output specific enthalpy

%---------------------------------------------------------------------------------------   

    case 'x_dt'         % Quality from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
        if isnan(Psat)          % Did the given temperature correspond to a point on the saturation line?
            X = NaN;                % If not, return NaN
        elseif rho < rhovsat    % Is the given state a pure vapor?
            X = 1;                  % If so, the quality is one
        elseif rho > rholsat    % Is the given state a pure liquid?
            X = 0;                  % If so, the quality is zero
        else                    % We must actually have a two-phase mixture, then
            v  = 1/rho;             % Convert given density to specific volume
            vl = 1/rholsat;         % Convert saturated liquid density to specific volume
            vv = 1/rhovsat;         % Convert saturated vapor density to specific volume
            X  = (v - vl)/(vv - vl); % Compute the quality from the specific volumes
        end
        OutVar = X;         % Output quality
        
%---------------------------------------------------------------------------------------   
        
    case 'cp_dt'        % Constant Pressure Specific Heat from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            cp = Water_cp(rho,T);   % Directly compute specific heat
        elseif (T > Tt)         % Regular Fluid
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            if (rho > rhom)             % High-Pressure Ice
                cp = NaN;                   % There's no correlation for this stuff...
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase Mixture
                cp = Inf;                   % By definition, the constant-pressure specific heat at saturation is infinite
            else                        % Single Fluid Phase
                cp = Water_cp(rho,T);       % Directly compute specific heat
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                cp = Water_cp(rho,T);       % Directly compute specific heat
            else                        % Ice
                cp = NaN;                   % Would need an appropriate correlation...
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                cp = Water_cp(rho,T);       % Directly compute specific heat
            else                        % Ice Ih
                cp = NaN;                   % Would need the ice correlation to compute this...
            end
        end
        OutVar = cp;        % Output constant-pressure specific heat
            
%---------------------------------------------------------------------------------------   
        
    case 'cv_dt'        % Constant Volume Specific Heat from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            cv = Water_cv(rho,T);   % Directly compute specific heat
        elseif (T > Tt)         % Regular Fluid
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find saturation properties at this temperature
            Pm = MeltPress(T);      % Find melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find liquid density at melting point
            if (rho > rhom)             % Ice
                cv = NaN;                   % Would need an appropriate correlation...
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                cv = (H2Oprop('u_dt',rho,T+1e6*eps(T)) - ...
                      H2Oprop('u_dt',rho,T))/1e6/eps(T); % Use an approximate numerical derivative to find the specific heat
            else                        % Single Phase
                cv = Water_cv(rho,T);       % Directly compute specific heat
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                cv = Water_cv(rho,T);       % Directly compute specific heat
            else                        % Ice
                cv = NaN;                   % Would need an appropriate correlation...
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                cv = Water_cv(rho,T);       % Directly compute specific heat
            else                        % Ice Ih
                cv = NaN;                   % Would need the ice correlation to compute this...
            end
        end
        OutVar = cv;        % Output constant volume specific heat
        
%---------------------------------------------------------------------------------------   

    case 'betat_dt'     % Thermal Expansion Coefficient from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % Ice
                betat = NaN;                % There's no correlation for this stuff...
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                betat = NaN;                % Thermal expansion is poorly defined for a saturated two-phase mixture
            else                        % Single Phase
                betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient
            else                    % Ice
                betat = NaN;            % Would need an appropriate correlation...
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient
            else                        % Ice Ih
                betat = NaN;                % Would need the ice correlation to compute this...
            end
        end
        OutVar = betat;     % Output thermal expansion coefficient
        
%---------------------------------------------------------------------------------------   

    case 'kappat_dt'    % Isothermal Compressibility from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            kappat = Water_kappat(rho,T); % Directly compute compressibility
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % High-Pressure Ice
                kappat = NaN;               % There's no correlation for this stuff...
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                kappat = NaN;               % Compressibility is poorly defined for a saturated mixture
            else                        % Single Phase
                kappat = Water_kappat(rho,T); % Directly compute compressibility
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                kappat = Water_kappat(rho,T); % Directly compute compressibility
            else                        % Ice
                kappat = NaN;               % Would need an appropriate correlation...
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                kappat = Water_kappat(rho,T); % Directly compute compressibility
            else                        % Ice Ih
                kappat = NaN;               % Would need the ice correlation to compute this...
            end
        end
        OutVar = kappat;    % Output isothermal compressibility
        
%---------------------------------------------------------------------------------------   

    case 'w_dt'     % Sound Velocity from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            w = Water_w(rho,T);     % Directly compute sound velocity
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % High-Pressure Ice
                w = NaN;                    % There's no correlation for this stuff...
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                w = NaN;                    % The liquid and vapor components have independent sound velocities
            else                        % Single Phase Fluid
                w = Water_w(rho,T);         % Directly compute sound velocity
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                w = Water_w(rho,T);         % Directly compute sound velocity
            else                        % Ice
                w = NaN;                    % Would need an appropriate correlation...
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                w = Water_w(rho,T);         % Directly compute sound velocity
            else                        % Ice Ih
                w = NaN;                    % Would need the ice correlation to compute this...
            end
        end
        OutVar = w;         % Output sound velocity
        
%---------------------------------------------------------------------------------------   

    case 'mu_dt'     % Viscosity from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            mu = Water_mu(rho,T);   % Directly compute viscosity
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % High-Pressure Ice
                mu = Inf;                   % Solid
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                mu = NaN;                   % The liquid and vapor components have independent viscosities
            else                        % Single Phase
                mu = Water_mu(rho,T);       % Directly compute viscosity
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                mu = Water_mu(rho,T);   % Directly compute viscosity
            else                    % Ice
                mu = Inf;               % Solid
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                mu = Water_mu(rho,T);       % Directly compute viscosity
            else                        % Ice Ih
                mu = Inf;                   % Solid
            end
        end
        OutVar = mu;        % Output viscosity
        
%---------------------------------------------------------------------------------------   

    case 'k_dt'     % Thermal Conductivity from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            k = Water_k(rho,T);     % Directly compute thermal conductivity
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % High-Pressure Ice
                k = NaN;                    % There's no correlation for this stuff
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                k = NaN;                    % The liquid and vapor components have independent thermal conductivities
            else                        % Single Phase
                k = Water_k(rho,T);         % Directly compute thermal conductivity
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                k = Water_k(rho,T);         % Directly compute thermal conductivity
            else                        % Ice
                k = NaN;                    % Would need an appropriate correlation...
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                k = Water_k(rho,T);         % Directly compute thermal conductivity
            else                        % Ice Ih
                k = NaN;                    % Would need an appropriate correlation...
            end
        end
        OutVar = k;         % Output thermal conductivity
        
%---------------------------------------------------------------------------------------   

    case 'sigma_t'      % Surface tension from temperature
        T = Var1;           % Input temperature
        if (T < Tt) || (T > Tc) % Is this temperature off of the saturation line?
            sigma = NaN;            % If so, there's no surface tension
        else
            sigma = Water_sigma(T); % Otherwise, directly compute surface tension of the saturated liquid
        end
        OutVar = sigma;     % Output surface tension

%---------------------------------------------------------------------------------------   

    case 'epsilon_dt'   % Dielectric Constant from density, temperature
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        if (T > Ticemax)        % Supercritical Fluid
            epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % High-Pressure Ice
                epsilon = NaN;              % There's no correlation for this stuff
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                epsilon = NaN;              % The liquid and vapor components have independent dielctric constants
            else                        % Single Phase
                epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant
            else                        % Ice
                epsilon = NaN;              % There's no correlation for this stuff
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant
            else                        % Ice Ih
                epsilon = NaN;              % There's no correlation for this stuff
            end
        end
        OutVar = epsilon;   % Output dielectric constant
        
%---------------------------------------------------------------------------------------   

    case 'n_dtl'     % Refractive Index from density, temperature, wavelength
        rho = Var1;         % Input density
        T   = Var2;         % Input temperature
        lambda = Var3;      % Input wavelength
        if (T > Ticemax)        % Supercritical Fluid
            n = Water_n(rho, T, lambda); % Directly compute refractive index
        elseif (T > Tt)         % Ordinary Fluid
            Pm = MeltPress(T);      % Find the melting pressure at this temperature
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point
            [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature
            if (rho > rhom)             % High-Pressure Ice
                n = NaN;                    % There's no correlation for this stuff
            elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase
                n = NaN;                    % The liquid and vapor components have independent refractive indices
            else                        % Single Phase
                n = Water_n(rho, T, lambda); % Directly compute refractive index
            end
        elseif (T > Tmin)       % Cold Fluid
            [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point
            rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point
            if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase
                n = Water_n(rho, T, lambda); % Directly compute refractive index
            else                        % Ice
                n = NaN;                    % There's no correlation for this stuff
            end
        elseif (T > Tabsmin)    % Cold Vapor
            Ps = MeltPress(T);      % Find the sublimation pressure at this temperature
            rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point
            if (rho <= rhos)            % Vapor
                n = Water_n(rho, T, lambda); % Directly compute refractive index
            else                        % Ice Ih
                n = NaN;                    % There's no correlation for this stuff
            end
        end
        OutVar = n;         % Output refractive index
        
%---------------------------------------------------------------------------------------   
        
    case 'tsat_p'       % Saturation temperature from pressure
        P = Var1;           % Input pressure
        if (P < Pt)             % Are we on the sublimation curve?
            T = MeltTemp(P);        % If so, find the sublimation temperature
        elseif (P < Pc)         % Are we on the saturation curve?
            T = fzero(@(x) GoodSatProp(x) - P, [Tt Tc - .01]); % If so, use root finding to find the saturation temperature
        else                    % Otherwise, there's no vapor phase boundary...
            T = NaN;                % ...so we return NaN
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   
        
    case 'tsat_dl'      % Saturation temperature from liquid density, valid above 277.134 K density anomaly
        rhol = Var1;        % Input saturated liquid density
        if (rhol > rhoc) && (rhol < rhoamax) % Are we on the saturation curve?
            T = fzero(@(x) GoodSatProp(x,'rholsat') - rhol, [Trhomax Tc - .01]); % If so, use root finding to find the saturation temperature
        else                    % Otherwise, the input was not a valid saturated liquid density...
            T = NaN;                % ...so we return NaN
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   
        
    case 'tsat_dv'      % Saturation temperature from vapor density
        rhov = Var1;        % Input saturated vapor density
        if (rhov < rhovt)       % Are we on the sublimation curve?
            T = fzero(@(x) Water_P(rhov, x) - MeltPress(x,'sub'), Tt); % If so, use root finding to find the sublimation temperature
        elseif (rhov < rhoc)    % Are we on the saturation curve?
            T = fzero(@(x) GoodSatProp(x,'rhovsat') - rhov, [Tt Tc]); % If so, use root finding to find the saturation temperature
        else                    % Otherwise, the input was not a valid saturated vapor density...
            T = NaN;                % ...so we return NaN
        end
        OutVar = T;         % Output temperature
        
%---------------------------------------------------------------------------------------   
        
    case 'psat_t'       % Saturation pressure from temperature
        T = Var1;           % Input temperature
        if (T < Tt)             % Are we on the sublimation curve?
            P = MeltPress(T,'sub'); % If so, find the sublimation pressure
        elseif (T < Tc)         % Are we on the saturation curve?
            P = GoodSatProp(T);     % if so, find the saturation pressure
        else                    % Otherwise, there's no vapor phase boundary...
            P = NaN;                % ...so we return NaN
        end
        OutVar = P;         % Output pressure
        
%---------------------------------------------------------------------------------------   
        
    case 'dlsat_t'      % Saturated liquid density from temperature
        T = Var1;           % Input temperature
        dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density
        OutVar = dl;        % Output saturated liquid density
        
%---------------------------------------------------------------------------------------   
        
    case 'dvsat_t'      % Saturated vapor density from temperature
        T = Var1;           % Input temperature
        if (T > Tt)             % Are we on or above the saturation curve?
            dv = GoodSatProp(T,'rhovsat'); % If so, find the saturated vapor density
        else                    % Otherwise, we must be on the sublimation curve
            Ps = MeltPress(T,'sub'); % Find the sublimation pressure
            dv = fzero(@(x) Water_P(x, T) - Ps, rhomin); % Use a root-finding method to find the vapor density at the sublimation point
        end
        OutVar = dv;        % Output saturated vapor density
        
%---------------------------------------------------------------------------------------   
        
    case 'dlsat_p'      % Saturated liquid density from pressure
        P = Var1;           % Input pressure
        if (P < Pt)             % Is the pressure too low for liquid to exist?
            P = NaN;                % If so, set it to NaN and let it propagate
        end
        T = H2Oprop('tsat_p',P); % Find the saturation temperature
        if isnan(T)             % Is the pressure supercritical, or was it NaN?
            dl = NaN;               % If so, return NaN
        else
            dl = GoodSatProp(T,'rholsat'); % Otherwise, calculate the saturated liquid density
        end
        OutVar = dl;        % Output saturated liquid density
        
%---------------------------------------------------------------------------------------   
        
    case 'dvsat_p'      % Saturated vapor density from pressure
        P = Var1;           % Input pressure
        T = H2Oprop('tsat_p',P); % Find the saturation or sublimation temperature
        if isnan(T)             % Does a saturation or sublimation temperature exist?
            dv = NaN;               % If not, return NaN
        elseif (P < Pt)         % Is the process of interest sublimation?
            dv = fzero(@(x) Water_P(x, T) - P, P/R/T); % If so, use a root-finding method to find the vapor density
        else                    % Otherwise, we must be looking at the saturation line
            dv = GoodSatProp(T,'rhovsat'); % And so, calculate the saturated vapor density
        end
        OutVar = dv;        % Output saturated vapor density
        
%---------------------------------------------------------------------------------------   
        
    case 'ulsat_t'      % Saturated liquid specific energy from temperature
        T = Var1;           % Input temperature
        dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density
        if isnan(dl)            % Does saturated liquid exist?
            ul = NaN;               % If not, return NaN
        else
            ul = Water_u(dl,T);     % Otherwise, directly compute the saturated liquid specific internal energy
        end
        OutVar = ul;        % Output saturated liquid specific internal energy
        
%---------------------------------------------------------------------------------------   
        
    case 'ulsat_p'      % Saturated liquid specific energy from pressure
        P = Var1;           % Input pressure
        if (P < Pt)             % Is the pressure too low for liquid to exist?
            P = NaN;                % If so, set it to NaN and let it propagate
        end
        T = H2Oprop('tsat_p',P); % Find the saturation temperature
        if isnan(T)             % Is the pressure supercritical, or was it NaN?
            ul = NaN;               % If so, return NaN
        else
            dl = GoodSatProp(T,'rholsat'); % Otherwise, find the saturated liquid density
            ul = Water_u(dl,T);     % Directly compute the saturated liquid specific internal energy
        end
        OutVar = ul;        % Output saturated liquid specific internal energy
        
%---------------------------------------------------------------------------------------   
        
    case 'uvsat_t'      % Saturated vapor specific energy from temperature
        T = Var1;           % Input temperature
        if (T < Tt)             % Are we on the sublimation curve?
            Ps = MeltPress(T,'sub'); % If so, find the sublimation pressure
            dv = fzero(@(x) Water_P(x, T) - Ps, rhomin); % Use a root-finding method to compute the saturated vapor density
            uv = Water_u(dv,T);     % Directly compute the saturated vapor specific internal energy
        elseif (T > Tc)         % Are we supercritical?
            uv = NaN;               % If so, return NaN
        else                    
            dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density
            uv = Water_u(dv,T);     % Directly compute the saturated vapor specific internal energy
        end
        OutVar = uv;        % Output saturated vapor specific internal energy
        
%---------------------------------------------------------------------------------------   
        
    case 'uvsat_p'      % Saturated vapor specific energy from pressure
        P = Var1;           % Input pressure
        T = H2Oprop('tsat_p',P); % Find the saturation or sublimation temperature
        if isnan(T)             % Does a saturation or sublimation temperature exist?
            uv = NaN;               % If not, return NaN
        elseif (P < Pt)         % Is the process of interest sublimation?
            dv = fzero(@(x) Water_P(x, T) - P, P/R/T); % If so, use a root-finding method to find the vapor density
            uv = Water_u(dv,T);     % Directly compute the saturated vapor specific internal energy
        else
            dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density
            uv = Water_u(dv,T);     % Directly compute the saturated vapor specific internal energy
        end
        OutVar = uv;        % Output saturated vapor specific internal energy
        
%---------------------------------------------------------------------------------------   
        
    case 'hlsat_t'      % Saturated liquid specific enthalpy from temperature
        T = Var1;           % Input temperature
        dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density
        if isnan(dl)            % Does saturated liquid exist?
            hl = NaN;               % If not, return NaN
        else
            hl = Water_h(dl,T);     % Otherwise, directly compute the saturated liquid specific enthalpy
        end
        OutVar = hl;        % Output saturated liquid specific enthalpy
        
%---------------------------------------------------------------------------------------   
        
    case 'hlsat_p'      % Saturated liquid specific enthalpy from pressure
        P = Var1;           % Input pressure
        if (P < Pt)             % Is the pressure too low for liquid to exist?
            P = NaN;                % If so, set it to NaN and let it propagate
        end
        T = H2Oprop('tsat_p',P); % Find the saturation temperature
        if isnan(T)             % Is the pressure supercritical, or was it NaN?
            hl = NaN;               % If so, return NaN
        else
            dl = GoodSatProp(T,'rholsat'); % Otherwise, find the saturated liquid density
            hl = Water_h(dl,T);     % Directly compute the saturated liquid specific enthalpy
        end
        OutVar = hl;        % Output saturated liquid specific enthalpy
        
%---------------------------------------------------------------------------------------   
        
    case 'hvsat_t'      % Saturated vapor specific enthalpy from temperature
        T = Var1;           % Input temperature
        if (T < Tt)             % Are we on the sublimation curve?
            Ps = MeltPress(T,'sub'); % If so, find the sublimation pressure
            dv = fzero(@(x) Water_P(x, T) - Ps, rhomin); % Use a root-finding method to compute the saturated vapor density
            hv = Water_h(dv,T);     % Directly compute the saturated vapor specific internal energy
        elseif (T > Tc)         % Are we supercritical?
            hv = NaN;               % If so, return NaN
        else                    
            dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density
            hv = Water_h(dv,T);     % Directly compute the saturated vapor specific internal energy
        end
        OutVar = hv;        % Output saturated vapor specific internal energy
        
%---------------------------------------------------------------------------------------   
        
    case 'hvsat_p'      % Saturated vapor specific enthalpy from pressure
        P = Var1;           % Input pressure
        T = H2Oprop('tsat_p',P); % Find the saturation or sublimation temperature
        if isnan(T)             % Does a saturation or sublimation temperature exist?
            hv = NaN;               % If not, return NaN
        elseif (P < Pt)         % Is the process of interest sublimation?
            dv = fzero(@(x) Water_P(x, T) - P, P/R/T); % If so, use a root-finding method to find the vapor density
            hv = Water_h(dv,T);     % Directly compute the saturated vapor specific enthalpy
        else
            dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density
            hv = Water_h(dv,T);     % Directly compute the saturated vapor specific enthalpy
        end
        OutVar = hv;        % Output saturated vapor specific enthalpy
        
%---------------------------------------------------------------------------------------   
        
    case 'slsat_t'      % Saturated liquid specific entropy from temperature
        T = Var1;           % Input temperature
        dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density
        if isnan(dl)            % Does saturated liquid exist?
            sl = NaN;               % If not, return NaN
        else
            sl = Water_s(dl,T);     % Otherwise, directly compute the saturated liquid specific entropy
        end
        OutVar = sl;        % Output saturated liquid specific entropy
        
%---------------------------------------------------------------------------------------   
        
    case 'slsat_p'      % Saturated liquid specific entropy from pressure
        P = Var1;           % Input pressure
        if (P < Pt)             % Is the pressure too low for liquid to exist?
            P = NaN;                % If so, set it to NaN and let it propagate
        end
        T = H2Oprop('tsat_p',P); % Find the saturation temperature
        if isnan(T)             % Is the pressure supercritical, or was it NaN?
            sl = NaN;               % If so, return NaN
        else
            dl = GoodSatProp(T,'rholsat'); % Otherwise, find the saturated liquid density
            sl = Water_s(dl,T);     % Directly compute the saturated liquid specific entropy
        end
        OutVar = sl;        % Output saturated liquid specific entropy
        
%---------------------------------------------------------------------------------------   
        
    case 'svsat_t'      % Saturated vapor specific entropy from temperature
        T = Var1;           % Input temperature
        dv = H2Oprop('dvsat_t',T); % Find the saturated vapor density
        if (T > Tc)             % If the temperature is supercritical...
            sv = NaN;               % ...return NaN
        else
            sv = Water_s(dv,T);     % Otherwise, directly compute the saturated vapor specific entropy
        end
        OutVar = sv;        % Output saturated vapor specific entropy
        
%---------------------------------------------------------------------------------------        

    case 'svsat_p'      % Saturated vapor specific entropy from pressure
        P = Var1;           % Input pressure
        T = H2Oprop('tsat_p',P); % Find the saturation temperature
        dv = H2Oprop('dvsat_p',P); % Find the saturated vapor density
        if isnan(T)             % If a saturation or sublimation temperature could not be found...
            sv = NaN;               % ...return NaN
        else
            sv = Water_s(dv,T);     % Otherwise, directly compute the saturated vapor specific entropy
        end
        OutVar = sv;        % Output saturated vapor specific entropy
    
%---------------------------------------------------------------------------------------
        
    otherwise
        help H2Oprop     % Display help text to assist with choice of valid query
        error('Invalid input string.') % Throw an exception to stop any code that might have called this function
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Internal Helper Functions

function P = Water_P(rho, T)        % Compute pressure using the equation of state

global Tc rhoc R                        % Inherit a few basic physical constants

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phir_d  = phirx_d(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to density

P = rho*R*T*(1 + delta*phir_d);         % Calculate the pressure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = Water_u(rho, T)        % Compute specific internal energy using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0_t  = n0(2) + n0(3)/tau + sum(n0(4:8).*gamma0(4:8).*((1 -...
           exp(-gamma0(4:8)*tau)).^(-1) - 1)); % Compute the partial derivative of the ideal gas part of Helmholtz energy with respect to temperature
       
phir_t  = phirx_t(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature

u = R*T*tau*(phi0_t + phir_t);          % Calculate the specific internal energy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function s = Water_s(rho, T)        % Compute specific entropy using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0    = log(delta) + n0(1) + n0(2)*tau + n0(3)*log(tau) + sum(n0(4:8).*log(1-...
           exp(-gamma0(4:8)*tau)));     % Compute the ideal gas part of Helmholtz energy
phi0_t  = n0(2) + n0(3)/tau + sum(n0(4:8).*gamma0(4:8).*((1 -...
           exp(-gamma0(4:8)*tau)).^(-1) - 1)); % Compute the partial derivative of the ideal gas part of Helmholtz energy with respect to temperature

phir    = phirx(delta,tau);             % Compute the residual part of Helmholtz energy
phir_t  = phirx_t(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature

s = R*tau*(phi0_t + phir_t) - R*phi0 - R*phir; % Calculate the specific entropy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h = Water_h(rho, T)        % Compute specific enthalpy using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0_t  = n0(2) + n0(3)/tau + sum(n0(4:8).*gamma0(4:8).*((1 -...
           exp(-gamma0(4:8)*tau)).^(-1) - 1)); % Compute the partial derivative of the ideal gas part of Helmholtz energy with respect to temperature

phir_d  = phirx_d(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
phir_t  = phirx_t(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature

h = R*T*(1 + tau*(phi0_t + phir_t) + delta*phir_d); % Calculate the specific enthalpy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cp = Water_cp(rho, T)      % Compute the constant-pressure specific heat using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*...
          (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature

phir_d  = phirx_d(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
phir_dd = phirx_dd(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density
phir_tt = phirx_tt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature
phir_dt = phirx_dt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature

cp = -R*tau^2*(phi0_tt + phir_tt) + R*(1 + delta*phir_d - delta*tau*phir_dt)^2/...
     (1 + 2*delta*phir_d + delta^2*phir_dd); % Calculate the constant-pressure specific heat
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function betat = Water_betat(rho, T) % Compute the thermal expansion coefficient using the equation of state

global Tc rhoc                          % Inherit a few basic physical constants

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phir_d  = phirx_d(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
phir_dd = phirx_dd(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density
phir_dt = phirx_dt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature

betat   = 1/T*(1+delta*phir_d-delta*tau*phir_dt)/(1+2*delta*phir_d+delta^2*phir_dd); % Calculate the thermal expansion coefficient

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cv = Water_cv(rho, T)      % Compute the constant-volume specific heat using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*...
          (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature

phir_tt = phirx_tt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature
         
cv = -R*tau^2*(phi0_tt + phir_tt);      % Calculate the constant-volume specific heat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function kappat = Water_kappat(rho, T) % Compute the isothermal compressibility using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*...
          (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature

phir_d  = phirx_d(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
phir_dd = phirx_dd(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density
phir_tt = phirx_tt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature
phir_dt = phirx_dt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature

% Properties

cv = -R*tau^2*(phi0_tt + phir_tt);              % Isochoric Heat Capacity
cp = -R*tau^2*(phi0_tt + phir_tt) + R*(1 + delta*phir_d - delta*tau*phir_dt)^2/...
     (1 + 2*delta*phir_d + delta^2*phir_dd);    % Isobaric Heat Capacity
w  = sqrt(R*T*(1 + 2*delta*phir_d + delta^2*phir_dd - (1 + delta*phir_d - ...
     delta*tau*phir_dt)^2/tau^2/(phi0_tt + phir_tt)));  % Speed of Sound
kappat  = cp/cv/rho/w^2;                % Isothermal Compressibility

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function w = Water_w(rho, T)        % Compute the speed of sound using the equation of state

global Tc rhoc R n0 gamma0              % Inherit a few basic physical constants and equation of state parameters

tau   = Tc/T;                           % Compute the reduced temperature
delta = rho/rhoc;                       % Compute the reduced density

phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*...
          (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature

phir_d  = phirx_d(delta,tau);           % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
phir_dd = phirx_dd(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density
phir_tt = phirx_tt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature
phir_dt = phirx_dt(delta,tau);          % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature

w  = sqrt(R*T*(1 + 2*delta*phir_d + delta^2*phir_dd - (1 + delta*phir_d - ...
     delta*tau*phir_dt)^2/tau^2/(phi0_tt + phir_tt)));  % Speed of Sound

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function sigma = Water_sigma(T)     % Compute the surface tension using the IAPWS 1994 formulation

global Tc                               % Inherit a few basic physical constants

TAU = 1 - T/Tc;                         % Compute the (alternate) reduced temperature
sigma = 0.2358*TAU^1.256*(1 - 0.625*TAU); % Use the IAPWS correlation to compute surface tension

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function mu = Water_mu(rho, T)      % Compute the viscosity using the IAPWS 2008 formulation

% *** NOTE: SOME DISAGREEMENT WITH DATA NEAR CRITICAL POINT
% Probably a subtle error in isothermal compressibility, shared with
% thermal conductivity formulation

global Tc rhoc R                        % Inherit a few basic physical constants

Tstar   = 647.096;                      % Reference Temperature
rhostar = 322.0;                        % Reference Density
Pstar   = 22.064e6;                     % Reference Pressure
mustar  = 1e-6;                         % Reference Viscosity

% Viscosity correlation constants

iv = [0 1 2 3 0 1 2 3 5 0 1 2 3 4 0 1 0 3 4 3 5]';
jv = [0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 4 4 5 6 6]';
Nv = ...
[ 5.20094e-1  8.50895e-2 -1.08374    -2.89555e-1  2.22531e-1  9.99115e-1 ...
  1.88797     1.26613     1.20573e-1 -2.81378e-1 -9.06851e-1 -7.72479e-1 ...
 -4.89837e-1 -2.57040e-1  1.61913e-1  2.57399e-1 -3.25372e-2  6.98452e-2 ...
  8.72102e-3 -4.35673e-3 -5.93264e-4]';

xi0    = 0.13e-9;                       
nu     = 0.630;
gamma0 = 1.239;
Gamma0 = 0.06;
xmu    = 0.068;                         % Critical exponent for viscosity
qC     = 1/1.9e-9;
qD     = 1/1.1e-9;
xi1    = 0.3817016416e-9;

% Viscosity correlation computation

Tbar   = T/Tstar;                       % Compute the (alternative) reduced temperature
Tbar1  = 1.5;                           % Define a reference reduced temperature
T1     = Tbar1*Tstar;                   % Compute the reference absolute temperature
rhobar = rho/rhostar;                   % Compute the (alternative) reduced density
tau    = Tc/T;                          % Compute the (conventional) reduced temperature
tau1   = Tc/T1;                         % Compute the reference (conventional) reduced temperature
delta  = rho/rhoc;                      % Compute the (conventional) reduced density

if (Tbar <= 1.01) && (Tbar >= 0.996) && (rhobar <= 1.36) && (rhobar >= 0.31) % Near the critical point?  If so, do the expensive critical enhancement calculation
    phir_d   = phirx_d(delta, tau);     % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
    phir_d1  = phirx_d(delta, tau1);    % Compute the partial derivative of the residual part of Helmholtz energy with respect to density at the reference temperature
    phir_dd  = phirx_dd(delta, tau);    % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density
    phir_dd1 = phirx_dd(delta, tau1);   % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density at the reference temperature
    chibar   = rhobar*Pstar/rhostar/R/T/(1 + 2*delta*phir_d + delta^2*phir_dd); % Compute the compressibility
    chibar1  = rhobar*Tbar1/Tbar * ...
               Pstar/rhostar/R/T1/(1 + 2*delta*phir_d1 + delta^2*phir_dd1); % Compute the compressibility at the reference temperature
    dchibar  = max((chibar - chibar1), 0); % Compute the compressibility difference with the reference temperature and discard mu2 if it's negative
    xi       = xi0*(dchibar/Gamma0)^(nu/gamma0); % Compute the correlation length
    if xi <= xi1                            % Use the below equation for Y if the correlation length is very short
        Y    = .2*qC*xi*(qD*xi)^5*(1 - qC*xi + (qC*xi)^2 - 765/504*(qD*xi)^2);
    else                                    % Otherwise, make the more complex calculation below for Y
        psiD = acos((1 + (qD*xi)^2)^-.5);
        w    = abs((qC*xi - 1)/(qC*xi + 1))^.5 * tan(psiD/2);
        L    = (qC*xi <= 1) * 2*atan(abs(w)) + (qC*xi > 1) * log((1 + w)/(1 - w));
        Y    = sin(3*psiD)/12 - sin(2*psiD)/4/qC/xi + ...
               (1 - 1.25*(qC*xi)^2)/(qC*xi)^2*sin(psiD) - ...
               ((1 - 1.5*(qC*xi)^2)*psiD - abs((qC*xi)^2 - 1)^1.5*L)/(qC*xi)^3;
    end
    mu2 = exp(xmu*Y);                   % Compute the critical enhancement for viscosity
else
    mu2 = 1;                            % Neglect any critical enhancement if far from the critical point
end

mu0 = 100*sqrt(Tbar)/(1.67752 + 2.20462/Tbar + 0.6366564/Tbar^2 - 0.241605/Tbar^3); % Compute dilute gas viscosity
mu1 = exp(rhobar*sum(Nv.*(1/Tbar - 1).^iv.*(rhobar - 1).^jv)); % Compute finite density viscosity enhancement

mu = mu0*mu1*mu2*mustar; % Use the IAPWS 2008 formulation to compute viscosity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function k = Water_k(rho, T)        % Compute the thermal conductivity using the IAPWS 2008 formulation

% ***  NOTE: SOME DISAGREEMENT WITH DATA NEAR CRITICAL POINT
% Probably a subtle error in isothermal compressibility, shared with
% viscosity formulation

global Tc rhoc R                        % Inherit a few basic physical constants

Tstar   = 647.226;                      % Reference Temperature
rhostar = 317.763;                      % Reference Density
Pstar   = 22.115e6;                     % Reference Pressure
kstar  = 0.4945;                        % Reference Thermal Conductivity

% Thermal conductivity correlation constants

ik = [0 0 0 0 0 0 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 4 4]';
jk = [0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 0 1]';
Nk = ...
[ 1.3293046  -0.40452437  0.24409490  0.018660751 -0.12961068  0.044809953 ...
  1.7018363  -2.2156845   1.6511057  -0.76736002   0.37283344 -0.11203160  ...
  5.2246158 -10.124111    4.9874687  -0.27297694  -0.43083393  0.13333849  ...
  8.7127675  -9.5000611   4.3786606  -0.91783782  -1.8525999   0.93404690]';

iv = [0 1 2 3 0 1 2 3 5 0 1 2 3 4 0 1 0 3 4 3 5]';
jv = [0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 4 4 5 6 6]';
Nv = ...
[ 5.20094e-1  8.50895e-2 -1.08384    -2.89555e-1  2.22531e-1  9.99115e-1 ...
  1.88797     1.26613     1.20573e-1 -2.81378e-1 -9.06851e-1 -7.72479e-1 ...
 -4.89837e-1 -2.57040e-1  1.61913e-1  2.57399e-1 -3.25372e-2  6.98452e-2 ...
  8.72102e-3 -4.35673e-3 -5.93264e-4]';

% Thermal conductivity correlation computation

Tbar   = T/Tstar;                       % Compute the (alternative) reduced temperature
rhobar = rho/rhostar;                   % Compute the (alternative) reduced density
tau    = Tc/T;                          % Compute the (conventional) reduced temperature
delta  = rho/rhoc;                      % Compute the (conventional) reduced density

phir_d  = phirx_d(delta, tau);          % Compute the partial derivative of the residual part of Helmholtz energy with respect to density
phir_dd = phirx_dd(delta, tau);         % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density
phir_dt = phirx_dt(delta, tau);         % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature
chibar  = rhobar*Pstar/rhostar/R/T/(1 + 2*delta*phir_d + delta^2*phir_dd); % Compute the compressibility
dPbar   = rho*R*Tstar/Pstar*(1 + delta*phir_d - delta*tau*phir_dt); % Compute the isochoric pressure-temperature coefficient

mu0 = 100*sqrt(Tbar)/(1.67752 + 2.20462/Tbar + 0.6366564/Tbar^2 - 0.241605/Tbar^3); % Compute dilute gas viscosity
mu1 = exp(rhobar*sum(Nv.*(1/Tbar - 1).^iv.*(rhobar - 1).^jv)); % Compute finite density viscosity enhancement

k0 = sqrt(Tbar)/(1 + 6.978267/Tbar + 2.599096/Tbar^2 - 0.998254/Tbar^3); % Compute dilute gas thermal conductivity
k1 = exp(rhobar*sum(Nk.*(1/Tbar - 1).^ik.*(rhobar - 1).^jk)); % Compute finite density thermal conductivity enhancement
k2 = 55.071*0.0013848/mu0/mu1*(Tbar/rhobar)^2*dPbar^2*chibar^0.4678*rhobar^0.5 * ...
     exp(-18.66*(Tbar - 1)^2 - (rhobar - 1)^4); % Compute critical enhancement in thermal conductivity

k = (k0*k1 + k2)*kstar;                 % Use the IAPWS 2008 formulation to compute thermal conductivity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function epsilon = Water_epsilon(rho, T) % Compute the dielectric constant using the IAPWS 1997 formulation

global Tc rhoc                          % Inherit a few basic physical constants

% Dielectric constant correlation constants

Nh = ...
[ 0.978224486826    -0.957771379375     0.237511794148     0.714692244396 ...
 -0.298217036956    -0.108863472196     0.949327488264e-1 -0.980469816509e-2 ...
  0.165167634970e-4  0.937359795772e-4 -0.123179218720e-9  0.196096504426e-2]';
ih = [1 1 1 2 3 3 4 5 6 7 10]';
jh = [0.25 1 2.5 1.5 1.5 2.5 2 2 5 0.5 10]';
Na = 6.0221367e23;
epsilon0 = 1/(4e-7*pi*299792458^2);
alpha = 1.636e-40;
Mw = 0.018015268;
mu = 6.138e-30;
kB = 1.380658e-23;

% Dielectric constant correlation computation

tau    = Tc/T;                          % Compute the reduced temperature
delta  = rho/rhoc;                      % Compute the reduced density

g = 1 + sum(Nh(1:11).*delta.^ih.*tau.^jh) + Nh(12)*delta*(T/228 - 1)^-1.2; % Compute the Harris-Alder g-factor
Ah = Na*mu^2/epsilon0/kB * rho*g/T/Mw;
Bh = Na*alpha/epsilon0/3 * rho/Mw;

epsilon = (1 + Ah + 5*Bh + sqrt(9 + 2*Ah + 18*Bh + Ah^2 + 10*Ah*Bh + 9*Bh^2))/ ...
          (4 - 4*Bh);                   % Use the IAPWS 1997 formulation to compute the dielectric constant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function n = Water_n(rho, T, lambda) % Compute the refractive index using the IAPWS 1997 formulation

Tbar = T/273.15;                        % Calculate the reduced temperature
rhobar = rho/1000;                      % Calculate the reduced density
lambdabar = lambda/589e-9;              % Calculate the reduced wavelength

nx = rhobar*(0.244257733 + 9.74634476e-3*rhobar - 3.73234996e-3*Tbar + ...
     2.68678472e-4*lambdabar^2*Tbar + 1.58920570e-3/lambdabar^2 + ...
     2.45934259e-3/(lambdabar^2-0.2292020^2) + ...
     0.900704920/(lambdabar^2-5.432937^2) - 1.66626219e-2*rhobar^2);

n = sqrt((1 + 2*nx)/(1 - nx));          % Use the IAPWS 1997 formulation to compute the refractive index

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Phase Boundary Functions

function [Psat, rholsat, rhovsat] = GoodSatProp(T, prop) % Use the Maxwell relation to find saturation properties from the equation of state

global Tc rhoc R                        % Inherit a few basic physical constants

tau = Tc/T;                             % Compute the reduced temperature

% Use the 1992 IAPWS correlation for saturation properties to determine a 
% starting guess for the saturation properties

[Psat rholsat rhovsat] = SatProp(T);    % Compute the approximate saturation properties
if isnan(Psat) || ((Tc - T) < 1e-5)         % If the given temperature is not on the saturation line, or is very close to the critical point, keep the approximate values
    return;
end

% Refine the initial estimate by iterating the Maxwell relation

tol1 = 1e-10;                           % Tolerance on iteration-to-iteration property change
tol2 = 1e-12;                           % Tolerance on change in iteration-to-iteration property change
err  = 10;                              % Initial iteration-to-iteration property change
err1 = 20;                              % Initial previous iteration-to-iteration property change
kmax = 20;                              % Maximum iteration count
k    = 0;                               % Starting iteration

while (err > tol1) && (k < kmax) && (abs(err - err1) > tol2)
    rholsat1 = fzero(@(x) (Psat - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), rholsat); % Find best saturated liquid density given Psat
    if (rholsat1 < rhoc)                    % Found wrong zero, which is possible near the critical point
        Guess      = R*T*rhoc*(1 + phirx_d(1,tau)); % Use the directly calculated pressure for critical density and the current temperature in place of Psat
        rholsat1   = fzero(@(x) (Guess - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), ...
                     [322.001 350]);            % Use firm bounds to ensure that the correct zero is located
    end
    rhovsat1 = fzero(@(x) (Psat - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), rhovsat); % Find best saturated vapor density given Psat
    if (rhovsat1 > rhoc)                    % Found wrong zero, which is possible near the critical point
        Guess      = R*T*rhoc*(1 + phirx_d(1,tau)); % Use the directly calculated pressure for critical density and the current temperature in place of Psat
        rhovsat1   = fzero(@(x) (Guess - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), ...
                     [300 321.999]);            % Use firm bounds to ensure that the correct zero is located
    end
    Psat1    = R*T/(1/rhovsat1 - 1/rholsat1)*(log(rholsat1/rhovsat1) + ...
               phirx(rholsat1/rhoc,tau) - phirx(rhovsat1/rhoc,tau)); % Use the new estimates of the saturation densities to improve the estimate of Psat
    errmat   = ([Psat rholsat rhovsat] - [Psat1 rholsat1 rhovsat1])./...
               [Psat rholsat rhovsat];  % Construct a matrix of the relative change in saturation properties
    err1     = err;                     % Store the previous iteration-to-iteration property change
    err      = max(abs(errmat));        % Note the largest relative iteration-to-iteration property change
    Psat     = Psat1;                   % Update the saturation pressure with the new estimate
    rholsat  = rholsat1;                % Update the saturated liquid density with the new estimate
    rhovsat  = rhovsat1;                % Update the saturated vapor density with the new estimate
    k        = k + 1;                   % Increment the iteration counter
end

if (nargin > 1)                             % If a specific property was requested in the first function output, put it there
    switch prop
        case 'rholsat'
            Psat = rholsat;                     % Return the saturated liquid density first
        case 'rhovsat'
            Psat = rhovsat;                     % Return the saturated vapor density first
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Psat rholsat rhovsat] = SatProp(T, prop) % Estimate the saturation properties using the 1992 IAPWS supplementary release

global Tc Tt Pc rhoc                    % Inherit a few basic physical constants

% Saturation property correlation constants

asat = [-7.85951783 1.84408259 -11.7866497 22.6807411 -15.9618719 1.80122502]';
Asat = [1 1.5 3 3.5 4 7.5]';
bsat = [1.99274064 1.09965342 -0.510839303 -1.75493479 -45.5170352 -6.74694450e5]';
Bsat = [1/3 2/3 5/3 16/3 43/3 110/3]';
csat = [-2.03150240 -2.68302940 -5.38626492 -17.2991605 -44.7586581 -63.9201063]';
Csat = [2/6 4/6 8/6 18/6 37/6 71/6]';

% Saturation propery computation

TAU  = 1 - T/Tc;                        % Calculate the reduced temperature

if (T > Tc) || (T < Tt)                 % If we aren't on the saturation line, return NaN
    Psat = NaN;
    rholsat = NaN;
    rhovsat = NaN;
else                                    % Otherwise, use the correlations to calculate the saturation properties
    Psat    = Pc*exp(Tc/T*sum(asat.*TAU.^Asat));
    rholsat = rhoc*(1 + sum(bsat.*TAU.^Bsat));
    rhovsat = rhoc*exp(sum(csat.*TAU.^Csat));
end

if (nargin > 1)                             % If a specific property was requested in the first function output, put it there
    switch prop
        case 'rholsat'
            Psat = rholsat;                     % Return the saturated liquid density first
        case 'rhovsat'
            Psat = rhovsat;                     % Return the saturated vapor density first
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function T = MeltTemp(P)

% 2008 IAPWS Formulation

% Describes melting curve from 273.16 K to 715 K
% Describes sublimation curve from 50 K to 273.16 K

% Triple Point Constants

P1 = 611.657;       %Ice I - Liquid - Vapor 
T1 = 273.16;        %Ice I - Liquid - Vapor

P2 = 208.566566032976716757e6; %Ice I - Ice III - Liquid (Adjusted to match correlations; nominally 208.566 MPa)
T2 = 251.165;       %Ice I - Ice III - Liquid

P3 = 350.1e6;       %Ice III - Ice V - Liquid
T3 = 256.164;       %Ice III - Ice V - Liquid

P4 = 632.4e6;       %Ice V - Ice VI - Liquid
T4 = 273.31;        %Ice V - Ice VI - Liquid

P5 = 2216e6;        %Ice VI - Ice VII - Liquid
T5 = 355;           %Ice VI - Ice VII - Liquid

Tmax = 715;         %Limits of Validity
Tmin = 50;

Pmax = 20617812820.4493;        % Highest pressure on melting line
Pmin = 1.9349587189864e-040;    % Lowest pressure on melting line

% Computation

if (P < Pmin)       %Below valid range
    T = NaN;
elseif (P < P1)     %Sublimation
    fun = @(x) P1*exp(T1/x*(-0.212144006e2*(x/T1)^(1/300) + ...
                        0.273203819e2*(x/T1)^(181/150) - ...
                        0.610598130e1*(x/T1)^(511/300)));
    T = fzero(@(x) fun(x) - P, [Tmin T1]); % Use a root-finding technique to find the sublimation temperature
elseif (P < P2)     %Ice Ih
    fun = @(x) P1*(1 + 0.119539337e7*(1 - (x/T1)^3) + ...
                   0.808183159e5*(1 - (x/T1)^25.75) + ...
                   0.333826860e4*(1 - (x/T1)^103.75));
    T = fzero(@(x) fun(x) - P, [T2 T1]); % Use a root-finding technique to find the Ice Ih melting point
elseif (P < P3)     %Ice III
    T = T2*((-0.700052 + P/P2)/0.299948)^(1/60); % Invert the correlation explicitly to find the Ice III melting point
elseif (P < P4)     %Ice V
    T = T3*((0.18721 + P/P3)/1.18721)^(1/8); % Invert the correlation explicitly to find the Ice V melting point
elseif (P < P5)     %Ice VI
    T = T4*((0.07476 + P/P4)/1.07476)^(1/4.6); % Invert the correlation explicitly to find the Ice VI melting point
elseif (P < Pmax)   %Ice VII
    fun = @(x) P5*exp(1.73683*(1-T5/x) - 0.0544606*(1-(x/T5)^5) + ...
               8.06106e-8*(1-(x/T5)^22));
    T = fzero(@(x) fun(x) - P, [T5 Tmax]); % Use a root-finding technique to find the Ice VII melting point
else                %Above valid range
    T = NaN;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Pm Pm1 Pm2] = MeltPress(T, prop)

% 2008 IAPWS Formulation

% Gives sublimation curve from 130 K to 251.165 K
% Gives Ice Ih curve from 251.165 K to 273.16 K (Plus other phases)
% Gives Ice V curve from 273.16 K to 273.31 K
% Gives Ice VI curve from 273.31 K to 355 K
% Gives Ice VII curve from 355 K to 715 K

% From 251.165 K to 273.16 K, Pm is the Ice Ih pressure, Pm1 is the
% sublimation pressure, and Pm2 is the Ice III or Ice V pressure.  This can
% be overridden by supplying the 'prop' argument, wthich can have one of
% two values: 'sub' returns sublimation pressure in Pm, and 'hp' returns
% the high-pressure ice melting mressure in Pm.  The 'prop' argument is
% useful in some minimization and root-finding problems, where only the
% first output parameter can be used.

% Triple Point Constants

P1 = 611.657;       % Ice I - Liquid - Vapor 
T1 = 273.16;        % Ice I - Liquid - Vapor

P2 = 208.566566032976716757e6; % Ice I - Ice III - Liquid (Adjusted to match correlations; nominally 208.566 MPa)
T2 = 251.165;       % Ice I - Ice III - Liquid

P3 = 350.1e6;       % Ice III - Ice V - Liquid
T3 = 256.164;       % Ice III - Ice V - Liquid

P4 = 632.4e6;       % Ice V - Ice VI - Liquid
T4 = 273.31;        % Ice V - Ice VI - Liquid

P5 = 2216e6;        % Ice VI - Ice VII - Liquid
T5 = 355;           % Ice VI - Ice VII - Liquid

Tabsmin = 50;       % Lowest Sublimation Temperature
Ticemax = 715;      % Highest Freezing Temperature

% Computation

if (T < Tabsmin)    % No Data
    Pm  = NaN;
    Pm1 = NaN;
    Pm2 = NaN;
elseif (T < T2)     % Sublimation Only
    Pm  = P1*exp(T1/T*(-0.212144006e2*(T/T1)^(1/300) + 0.273203819e2*(T/T1)^(181/150) - ...
                        0.610598130e1*(T/T1)^(511/300))); % Report sublimation pressure in Pm
    Pm1 = NaN;
    Pm2 = NaN;
elseif (T < T3)         % Ice Ih + Sublimation + Ice III
    Pm  = P1*(1 + 0.119539337e7*(1 - (T/T1)^3) + 0.808183159e5*(1 - (T/T1)^25.75) + ...
              0.333826860e4*(1 - (T/T1)^103.75)); % Report Ice Ih pressure in Pm
    Pm1 = P1*exp(T1/T*(-0.212144006e2*(T/T1)^(1/300) + ...
                        0.273203819e2*(T/T1)^(181/150) - ...
                        0.610598130e1*(T/T1)^(511/300))); % Report sublimation pressure in Pm1
    Pm2 = P2*(1 - 0.299948*(1 - (T/T2)^60)); % Report Ice III pressure in Pm2
elseif (T <= T1)        % Ice Ih + Sublimation + Ice V
    Pm  = P1*(1 + 0.119539337e7*(1 - (T/T1)^3) + 0.808183159e5*(1 - (T/T1)^25.75) + ...
              0.333826860e4*(1 - (T/T1)^103.75)); % Report Ice Ih pressure in Pm
    Pm1 = P1*exp(T1/T*(-0.212144006e2*(T/T1)^(1/300) + ...
                        0.273203819e2*(T/T1)^(181/150) - ...
                        0.610598130e1*(T/T1)^(511/300))); % Report sublimation pressure in Pm1
    Pm2 = P3*(1 - 1.18721*(1 - (T/T3)^8)); % Report Ice V pressure in Pm2
elseif (T < T4)         % Ice V
    Pm  = P3*(1 - 1.18721*(1 - (T/T3)^8)); % Report Ice V pressure in Pm
    Pm1 = NaN;
    Pm2 = NaN;
elseif (T < T5)         % Ice VI
    Pm  = P4*(1 - 1.07476*(1 - (T/T4)^4.6)); % Report Ice VI pressure in Pm
    Pm1 = NaN;
    Pm2 = NaN;
elseif (T <= Ticemax)   % Ice VII
    Pm  = P5*exp(1.73683*(1 - T5/T) - 0.0544606*(1 - (T/T5)^5) + ...
          8.06106e-8*(1 - (T/T5)^22)); % Report Ice VII pressure in Pm
    Pm1 = NaN;
    Pm2 = NaN;
else                    % Temperature above experimental data; melting line unknown
    Pm  = NaN;
    Pm1 = NaN;
    Pm2 = NaN;
end

if (nargin > 1) && ~isnan(Pm1)  % Only listen to 'prop' parameter if it's relevant
    switch prop
        case 'sub'      % Report sublimation pressure in Pm
            Pm = Pm1;
        case 'hp'       % Report high-pressure ice melting pressure in Pm
            Pm = Pm2;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Low-level Equation of State Functions

% Consult the IAPWS 1995 release for details of these functions

function Q = phirx(delta, tau)      % Compute the residual part of Helmholtz energy

global c d t a b B C D A alpha beta gamma epsilon n

psi    = exp(-C*(delta - 1)^2 - D*(tau - 1)^2);
theta  = 1 - tau + A.*((delta - 1)^2).^(1/2./beta);
DELTA  = theta.^2 + B.*((delta - 1)^2).^a;
DELTAb = DELTA.^b;

Q      = sum(n(1:7).*delta.^d(1:7).*tau.^t(1:7)) + ...
         sum(n(8:51).*delta.^d(8:51).*tau.^t(8:51).*exp(-delta.^c(8:51))) + ...
         sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ...
       	   epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2)) + ...
         sum(n(55:56).*DELTAb(55:56).*delta.*psi(55:56));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Q = phirx_d(delta, tau)    % Compute the partial derivative of the residual part of Helmholtz energy with respect to density

global c d t a b B C D A alpha beta gamma epsilon n

psi      = exp(-C*(delta - 1)^2 - D*(tau - 1)^2);
psi_d    = -2.*C*(delta - 1).*psi;
theta    = 1 - tau + A.*((delta - 1)^2).^(1/2./beta);
DELTA    = theta.^2 + B.*((delta - 1)^2).^a;
DELTAb   = DELTA.^b;
DELTA_d  = (delta - 1)*(2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ...
           2*B.*a.*((delta - 1)^2).^(a-1));
DELTAb_d = b.*DELTA.^(b - 1).*DELTA_d;

Q        = sum(n(1:7).*d(1:7).*delta.^(d(1:7) - 1).*tau.^t(1:7)) + ...
           sum(n(8:51).*exp(-delta.^c(8:51)).*delta.^(d(8:51) - 1).*tau.^t(8:51).*...
             (d(8:51)-c(8:51).*delta.^c(8:51))) + ...
           sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ...
			 epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*(d(52:54)/...
             delta - 2*alpha(52:54).*(delta - epsilon(52:54)))) + ...
           sum(n(55:56).*(DELTAb(55:56).*(psi(55:56) + delta*psi_d(55:56)) + ...
             delta*DELTAb_d(55:56).*psi(55:56)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Q = phirx_dd(delta, tau)   % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density

global c d t a b B C D A alpha beta gamma epsilon n

psi       = exp(-C*(delta - 1)^2 - D*(tau - 1)^2);
psi_d     = -2.*C*(delta - 1).*psi;
psi_dd    = 2*C.*psi.*(2*C*(delta - 1)^2 - 1);
theta     = 1 - tau + A.*((delta - 1)^2).^(1/2./beta);
DELTA     = theta.^2 + B.*((delta - 1)^2).^a;
DELTAb    = DELTA.^b;
DELTA_d   = (delta - 1)*(2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ...
            2*B.*a.*((delta - 1)^2).^(a-1));
DELTAb_d  = b.*DELTA.^(b - 1).*DELTA_d;
DELTA_dd  = 2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ...
            2*B.*a.*((delta - 1)^2).^(a-1) + (delta - 1)^2*(4*B.*a.*(a - 1).*...
		      ((delta - 1)^2).^(a - 2) + 2*A.^2.*beta.^-2.*(((delta - 1)^2).^...
            (1/2./beta - 1)).^2) + 4*A.*theta./beta.*(2./beta - 1).*...
            ((delta - 1)^2).^(1/2./beta);
DELTAb_dd = b.*(DELTA.^(b - 1).*DELTA_dd + (b - 1).*DELTA.^(b - 2).*DELTA_d.^2);

Q       = sum(n(1:7).*d(1:7).*(d(1:7) - 1).*delta.^(d(1:7) - 2).*tau.^t(1:7)) + ...
          sum(n(8:51).*exp(-delta.^c(8:51)).*delta.^(d(8:51) - 2).*tau.^t(8:51).*...
            ((d(8:51) - c(8:51).*delta.^c(8:51)).*(d(8:51) - 1 - ...
            c(8:51).*delta.^c(8:51)) - c(8:51).^2.*delta.^c(8:51))) + ...
          sum(n(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ...
            epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*...
            (-2*alpha(52:54).*delta.^d(52:54) + 4*alpha(52:54).^2.*delta.^d(52:54).*...
            (delta - epsilon(52:54)).^2 - 4*d(52:54).*alpha(52:54).*...
            delta.^(d(52:54) - 1).*(delta - epsilon(52:54)) + d(52:54).*(d(52:54) - ...
            1).*delta.^(d(52:54) - 2))) + ...
          sum(n(55:56).*(DELTAb(55:56).*(2*psi_d(55:56) + ...
            delta*psi_dd(55:56)) + 2*DELTAb_d(55:56).*(psi(55:56) + ...
            delta*psi_d(55:56)) + delta*DELTAb_dd(55:56).*psi(55:56)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Q = phirx_t(delta, tau)    % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature

global c d t a b B C D A alpha beta gamma epsilon n

psi      = exp(-C*(delta - 1)^2 - D*(tau - 1)^2);
psi_t    = -2*D.*psi*(tau - 1);
theta    = 1 - tau + A.*((delta - 1)^2).^(1/2./beta);
DELTA    = theta.^2 + B.*((delta - 1)^2).^a;
DELTAb   = DELTA.^b;
DELTAb_t = -2*theta.*b.*DELTA.^(b - 1);

Q      = sum(n(1:7).*t(1:7).*delta.^d(1:7).*tau.^(t(1:7) - 1)) + ...
         sum(n(8:51).*t(8:51).*delta.^d(8:51).*tau.^(t(8:51) - 1).*...
           exp(-delta.^c(8:51))) + ...
         sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ...
           epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*...
           (t(52:54)/tau - 2*beta(52:54).*(tau - gamma(52:54)))) + ...
         sum(n(55:56)*delta.*(DELTAb_t(55:56).*psi(55:56) + DELTAb(55:56).*...
           psi_t(55:56)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Q = phirx_tt(delta, tau)   % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature

global c d t a b B C D A alpha beta gamma epsilon n

psi       = exp(-C*(delta - 1)^2 - D*(tau - 1)^2);
psi_t     = -2*D.*psi*(tau - 1);
psi_tt    = 2*D.*psi.*(2*D*(tau - 1)^2 - 1);
theta     = 1 - tau + A.*((delta - 1)^2).^(1/2./beta);
DELTA     = theta.^2 + B.*((delta - 1)^2).^a;
DELTAb    = DELTA.^b;
DELTAb_t  = -2*theta.*b.*DELTA.^(b - 1);
DELTAb_tt = 2*b.*DELTA.^(b - 1) + 4*theta.^2.*b.*(b - 1).*DELTA.^(b - 2);

Q       = sum(n(1:7).*t(1:7).*(t(1:7) - 1).*delta.^d(1:7).*tau.^(t(1:7) - 2)) + ...
          sum(n(8:51).*t(8:51).*(t(8:51) - 1).*delta.^d(8:51).*tau.^(t(8:51) - 2).*...
            exp(-delta.^c(8:51))) + ...
          sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ...
            epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*...
            ((t(52:54)/tau - 2*beta(52:54).*(tau - gamma(52:54))).^2 - ...
            t(52:54)/tau^2 - 2*beta(52:54))) + ...
          sum(n(55:56)*delta.*(DELTAb_tt(55:56).*psi(55:56) + 2*DELTAb_t(55:56).*...
            psi_t(55:56) + DELTAb(55:56).*psi_tt(55:56)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Q = phirx_dt(delta, tau)   % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature

global c d t a b B C D A alpha beta gamma epsilon n

psi       = exp(-C*(delta - 1)^2 - D*(tau - 1)^2);
psi_d     = -2.*C*(delta - 1).*psi;
psi_t     = -2*D.*psi*(tau - 1);
psi_dt    = 4*C.*D.*psi*(delta - 1)*(tau - 1);
theta     = 1 - tau + A.*((delta - 1)^2).^(1/2./beta);
DELTA     = theta.^2 + B.*((delta - 1)^2).^a;
DELTA_d   = (delta - 1)*(2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ...
            2*B.*a.*((delta - 1)^2).^(a-1));
DELTAb    = DELTA.^b;
DELTAb_d  = b.*DELTA.^(b - 1).*DELTA_d;
DELTAb_t  = -2*theta.*b.*DELTA.^(b - 1);
DELTAb_dt = -2*(delta - 1)*A.*b./beta.*DELTA.^(b - 1).*((delta - 1)^2).^...
			   (1/2./beta - 1) - 2*theta.*b.*(b-1).*DELTA.^(b - 2).*DELTA_d;

Q       = sum(n(1:7).*d(1:7).*t(1:7).*delta.^(d(1:7) - 1).*tau.^(t(1:7) - 1)) + ...
          sum(n(8:51).*t(8:51).*delta.^(d(8:51) - 1).*tau.^(t(8:51) - 1).*(d(8:51) - ...
            c(8:51).*delta.^c(8:51)).*exp(-delta.^c(8:51))) + ...
          sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ...
            epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*...
            (d(52:54)/delta - 2*alpha(52:54).*(delta - epsilon(52:54))).*...
            (t(52:54)/tau - 2*beta(52:54).*(tau - gamma(52:54)))) + ...
          sum(n(55:56).*(DELTAb(55:56).*(psi_t(55:56) + delta*psi_dt(55:56)) + ...
            delta*DELTAb_d(55:56).*psi_t(55:56) + DELTAb_t(55:56).*(psi(55:56) + ...
            delta*psi_d(55:56)) + delta*DELTAb_dt(55:56).*psi(55:56)));

