MATLAB/Engineering thermodynamics

< MATLAB

This page is devoted to XSteam.m, which is a Matlab code available available on the internet.^[1]

There seems to be a flaw in T(h,s) in MATLAB steam table code found at http://www.mathworks.com/matlabcentral/fileexchange/9817-x-steam--thermodynamic-properties-of-water-and-steam

Matlab code to generate steam tables

The code that follows was written by Magnus Holmgren, who generously released it under conditions that permit it to be posted here. The following link should be used to retrieve the code because it contains valuable instructions on using it:

http://www.mathworks.com/matlabcentral/fileexchange/9817-x-steam--thermodynamic-properties-of-water-and-steam

This website also has a US units version, as well as detailed instructions on how to call the function. Once you have this code in a MATLAB folder with the title

 XSteam.m

you can call it from the command window with a variety of arguments. The following two examples are taken from the pdf file that accompanies the download:

 XSteam('h_pt',1,20)

Returns the enthalpy of water at 1 bar and 20 C.

 XSteam('TSat_p',1)

Returns the saturation temperature of water at 1 bar.

Backup copy of XSteam.m

Click to view code, but if possible, download it the aforementioned website

%h_prho behöver T_prho för samtliga regioner!!!!
 
%***********************************************************************************************************
%* Water and steam properties according to IAPWS IF-97                                                     *
%* By Magnus Holmgren, www.x-eng.com                                                                       *
%* The steam tables are free and provided as is.                                                           *
%* We take no responsibilities for any errors in the code or damage thereby.                               *
%* You are free to use, modify and distribute the code as long as authorship is properly acknowledged.     *
%* Please notify me at magnus@x-eng.com if the code is used in commercial applications                     *
%***********************************************************************************************************
%
% XSteam provides accurate steam and water properties from 0 - 1000 bar and from 0 - 2000 deg C according to
% the standard IAPWS IF-97. For accuracy of the functions in different regions see IF-97 (www.iapws.org)
%
% *** Using XSteam *****************************************************************************************
%XSteam take 2 or 3 arguments. The first argument must always be the steam table function you want to use.
%The other arguments are the inputs to that function.
%Example: XSteam('h_pt',1,20)  Returns the enthalpy of water at 1 bar and 20 degC
%Example: XSteam('TSat_p',1)  Returns the saturation temperature of water at 1 bar.
%For a list of valid Steam Table functions se bellow or the XSteam macros for MS Excel.
%
%*** Nomenclature ******************************************************************************************
% First the wanted property then a _ then the wanted input properties. 
% Example. T_ph is temperature as a function of pressure and enthalpy. 
% For a list of valid functions se bellow or XSteam for MS Excel.
% T     Temperature (deg C)
% p	    Pressure    (bar)
% h	    Enthalpy    (kJ/kg)
% v	    Specific volume (m3/kg)
% rho	Density
% s	    Specific entropy 
% u	    Specific internal energy 
% Cp	Specific isobaric heat capacity 
% Cv	Specific isochoric heat capacity 
% w	    Speed of sound 
% my	Viscosity
% tc	Thermal Conductivity
% st	Surface Tension
% x	    Vapour fraction
% vx	Vapour Volume Fraction
%
%*** Valid Steam table functions. ****************************************************************************
%
%Temperature	
%Tsat_p	Saturation temperature
%T_ph	Temperture as a function of pressure and enthalpy
%T_ps	Temperture as a function of pressure and entropy
%T_hs	Temperture as a function of enthalpy and entropy
%
%Pressure	
%psat_T	Saturation pressure
%p_hs	Pressure as a function of h and s. 
%p_hrho Pressure as a function of h and rho. Very unaccurate for solid water region since it's almost incompressible!
%
%Enthalpy	
%hV_p	Saturated vapour enthalpy
%hL_p	Saturated liquid enthalpy
%hV_T	Saturated vapour enthalpy
%hL_T	Saturated liquid enthalpy
%h_pT	Entalpy as a function of pressure and temperature.
%h_ps	Entalpy as a function of pressure and entropy.
%h_px	Entalpy as a function of pressure and vapour fraction
%h_prho	Entalpy as a function of pressure and density. Observe for low temperatures (liquid) this equation has 2 solutions. 
%h_Tx	Entalpy as a function of temperature and vapour fraction
%
%Specific volume	
%vV_p	Saturated vapour volume
%vL_p	Saturated liquid volume
%vV_T	Saturated vapour volume
%vL_T	Saturated liquid volume
%v_pT	Specific volume as a function of pressure and temperature.
%v_ph	Specific volume as a function of pressure and enthalpy
%v_ps	Specific volume as a function of pressure and entropy.
%
%Density	
%rhoV_p	Saturated vapour density
%rhoL_p	Saturated liquid density
%rhoV_T	Saturated vapour density
%rhoL_T	Saturated liquid density
%rho_pT	Density as a function of pressure and temperature.
%rho_ph	Density as a function of pressure and enthalpy
%rho_ps	Density as a function of pressure and entropy.
%
%Specific entropy 	
%sV_p	Saturated vapour entropy
%sL_p	Saturated liquid entropy
%sV_T	Saturated vapour entropy
%sL_T	Saturated liquid entropy
%s_pT	Specific entropy as a function of pressure and temperature (Returns saturated vapour entalpy if mixture.)
%s_ph	Specific entropy as a function of pressure and enthalpy
%
%Specific internal energy 	
%uV_p	Saturated vapour internal energy
%uL_p	Saturated liquid internal energy
%uV_T	Saturated vapour internal energy
%uL_T	Saturated liquid internal energy
%u_pT	Specific internal energy as a function of pressure and temperature.
%u_ph	Specific internal energy as a function of pressure and enthalpy
%u_ps	Specific internal energy as a function of pressure and entropy.
%
%Specific isobaric heat capacity 	
%CpV_p	Saturated vapour heat capacity 
%CpL_p	Saturated liquid heat capacity 
%CpV_T	Saturated vapour heat capacity 
%CpL_T	Saturated liquid heat capacity 
%Cp_pT	Specific isobaric heat capacity as a function of pressure and temperature.
%Cp_ph	Specific isobaric heat capacity as a function of pressure and enthalpy
%Cp_ps	Specific isobaric heat capacity as a function of pressure and entropy.
%
%Specific isochoric heat capacity 	
%CvV_p	Saturated vapour isochoric heat capacity
%CvL_p	Saturated liquid isochoric heat capacity
%CvV_T	Saturated vapour isochoric heat capacity
%CvL_T	Saturated liquid isochoric heat capacity
%Cv_pT	Specific isochoric heat capacity as a function of pressure and temperature.
%Cv_ph	Specific isochoric heat capacity as a function of pressure and enthalpy
%Cv_ps	Specific isochoric heat capacity as a function of pressure and entropy.
%
%Speed of sound 	
%wV_p	Saturated vapour speed of sound
%wL_p	Saturated liquid speed of sound
%wV_T	Saturated vapour speed of sound
%wL_T	Saturated liquid speed of sound
%w_pT	Speed of sound as a function of pressure and temperature.
%w_ph	Speed of sound as a function of pressure and enthalpy
%w_ps	Speed of sound as a function of pressure and entropy.
%
%Viscosity	
%Viscosity is not part of IAPWS Steam IF97. Equations from 
%"Revised Release on the IAPWS Formulation 1985 for the Viscosity of Ordinary Water Substance", 2003 are used.	
%Viscosity in the mixed region (4) is interpolated according to the density. This is not true since it will be two fases.	
%my_pT	Viscosity as a function of pressure and temperature.
%my_ph	Viscosity as a function of pressure and enthalpy
%my_ps	Viscosity as a function of pressure and entropy.
%
%Thermal Conductivity	
%Revised release on the IAPS Formulation 1985 for the Thermal Conductivity of ordinary water substance (IAPWS 1998)	
%tcL_p	Saturated vapour thermal conductivity
%tcV_p	Saturated liquid thermal conductivity
%tcL_T	Saturated vapour thermal conductivity
%tcV_T	Saturated liquid thermal conductivity
%tc_pT	Thermal conductivity as a function of pressure and temperature.
%tc_ph	Thermal conductivity as a function of pressure and enthalpy
%tc_hs	Thermal conductivity as a function of enthalpy and entropy
%
%Surface tension
%st_T	Surface tension for two phase water/steam as a function of T
%st_p	Surface tension for two phase water/steam as a function of T
%Vapour fraction
%x_ph	Vapour fraction as a function of pressure and enthalpy
%x_ps	Vapour fraction as a function of pressure and entropy.
%
%Vapour volume fraction
%vx_ph	Vapour volume fraction as a function of pressure and enthalpy
%vx_ps	Vapour volume fraction as a function of pressure and entropy.
 
 
function Out=XSteam(fun,In1,In2)
%*Contents.
%*1 Calling functions
%*1.1
%*1.2 Temperature (T)
%*1.3 Pressure (p)
%*1.4 Enthalpy (h)
%*1.5 Specific Volume (v)
%*1.6 Density (rho)
%*1.7 Specific entropy (s)
%*1.8 Specific internal energy (u)
%*1.9 Specific isobaric heat capacity (Cp)
%*1.10 Specific isochoric heat capacity (Cv)
%*1.11 Speed of sound
%*1.12 Viscosity
%*1.13 Prandtl
%*1.14 Kappa
%*1.15 Surface tension
%*1.16 Heat conductivity
%*1.17 Vapour fraction
%*1.18 Vapour Volume Fraction
%
%*2 IAPWS IF 97 Calling functions
%*2.1 Functions for region 1
%*2.2 Functions for region 2
%*2.3 Functions for region 3
%*2.4 Functions for region 4
%*2.5 Functions for region 5
%
%*3 Region Selection
%*3.1 Regions as a function of pT
%*3.2 Regions as a function of ph
%*3.3 Regions as a function of ps
%*3.4 Regions as a function of hs
%
%4 Region Borders
%4.1 Boundary between region 1 and 3.
%4.2 Region 3. pSat_h and pSat_s
%4.3 Region boundary 1to3 and 3to2 as a functions of s
%
%5 Transport properties
%5.1 Viscosity (IAPWS formulation 1985)
%5.2 Thermal Conductivity (IAPWS formulation 1985)
%5.3 Surface Tension
%
%6 Units
%
%7 Verification
%7.1 Verifiy region 1
%7.2 Verifiy region 2
%7.3 Verifiy region 3
%7.4 Verifiy region 4
%7.5 Verifiy region 5
 
%***********************************************************************************************************
%*1 Calling functions                                                                                      *
%***********************************************************************************************************
 
%***********************************************************************************************************
%*1.1
 
 
fun=lower(fun);
 
switch fun
    %***********************************************************************************************************
    %*1.2 Temperature
case 'tsat_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        Out = fromSIunit_T(T4_p(p));
    else
        Out = NaN;
    end
case 'tsat_s'
    s = toSIunit_s(In1);
    if s > -0.0001545495919  && s < 9.155759395 
        ps = p4_s(s);
        Out = fromSIunit_T(T4_p(ps));
    else
        Out = NaN;
    end
case 't_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        Out = fromSIunit_T(T1_ph(p, h));
    case 2
        Out = fromSIunit_T(T2_ph(p, h));
    case 3
        Out = fromSIunit_T(T3_ph(p, h));
    case 4
        Out = fromSIunit_T(T4_p(p));
    case 5
        Out = fromSIunit_T(T5_ph(p, h));
    otherwise
        Out = NaN;
    end
case 't_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Out = fromSIunit_T(T1_ps(p, s));
    case 2
        Out = fromSIunit_T(T2_ps(p, s));
    case 3
        Out = fromSIunit_T(T3_ps(p, s));
    case 4
        Out = fromSIunit_T(T4_p(p));
    case 5
        Out = fromSIunit_T(T5_ps(p, s));
    otherwise
        Out = NaN;
    end 
 
case 't_hs'
    h = toSIunit_h(In1);
    s = toSIunit_s(In2);
    Region = region_hs(h, s);
    switch Region
    case 1
        p1 = p1_hs(h, s);
        Out = fromSIunit_T(T1_ph(p1, h));
    case 2
        p2 = p2_hs(h, s);
        Out = fromSIunit_T(T2_ph(p2, h));
    case 3
        p3 = p3_hs(h, s);
        Out = fromSIunit_T(T3_ph(p3, h));
    case 4
        Out = fromSIunit_T(T4_hs(h, s));
    case 5
        error('functions of hs is not avlaible in region 5');
    otherwise
        Out = NaN;
    end
    %***********************************************************************************************************
    %*1.3 Pressure (p)
case 'psat_t'
    T = toSIunit_T(In1);
    if T < 647.096  && T > 273.15 
        Out = fromSIunit_p(p4_T(T));
    else
        Out = NaN;
    end
 
case 'psat_s'
    s = toSIunit_s(In1);
    if s > -0.0001545495919  && s < 9.155759395 
        Out = fromSIunit_p(p4_s(s));
    else
        Out = NaN;
    end
 
case 'p_hs'
    h = toSIunit_h(In1);
    s = toSIunit_s(In2);
    Region = region_hs(h, s);
    switch Region
    case 1
        Out = fromSIunit_p(p1_hs(h, s));
    case 2
        Out = fromSIunit_p(p2_hs(h, s));
    case 3
        Out = fromSIunit_p(p3_hs(h, s));
    case 4
        tSat = T4_hs(h, s);
        Out = fromSIunit_p(p4_T(tSat));
    case 5
        error('functions of hs is not avlaible in region 5');
    otherwise
        Out = NaN;
    end
 
    case 'p_hrho'
    h=In1;
    rho=In2;
%Not valid for water or sumpercritical since water rho does not change very much with p.
%Uses iteration to find p.
  High_Bound = fromSIunit_p(100);
  Low_Bound = fromSIunit_p(0.000611657);
  ps = fromSIunit_p(10);
  rhos = 1 / XSteam('v_ph',ps, h);
  while abs(rho - rhos) > 0.0000001
    rhos = 1 / XSteam('v_ph',ps, h);
    if rhos >= rho
      High_Bound = ps;
    else
      Low_Bound = ps;
    end
    ps = (Low_Bound + High_Bound) / 2;
  end
    Out = ps;
 
    %***********************************************************************************************************
    %*1.4 Enthalpy (h)
case 'hv_p'
    p = toSIunit_p(In1);
    if p > 0.000611657 && p < 22.06395 
        Out = fromSIunit_h(h4V_p(p));
    else
        Out = NaN;
    end
 
case 'hl_p'
    p = toSIunit_p(In1);
    if p > 0.000611657 && p < 22.06395 
        Out = fromSIunit_h(h4L_p(p));
    else
        Out = NaN;
    end
 
case 'hv_t'
    T = toSIunit_T(In1);
    if T > 273.15 && T < 647.096 
        p = p4_T(T);
        Out = fromSIunit_h(h4V_p(p));
    else
        Out = NaN;
    end
 
case 'hl_t'
    T = toSIunit_T(In1);
    if T > 273.15 && T < 647.096 
        p = p4_T(T);
        Out = fromSIunit_h(h4L_p(p));
    else
        Out = NaN;
    end
 
 
case 'h_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_h(h1_pT(p, T));
    case 2
        Out = fromSIunit_h(h2_pT(p, T));
    case 3
        Out = fromSIunit_h(h3_pT(p, T));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_h(h5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
case 'h_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Out = fromSIunit_h(h1_pT(p, T1_ps(p, s)));
    case 2
        Out = fromSIunit_h(h2_pT(p, T2_ps(p, s)));
    case 3
        Out = fromSIunit_h(h3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)));
    case 4
        xs = x4_ps(p, s);
        Out = fromSIunit_h(xs * h4V_p(p) + (1 - xs) * h4L_p(p));
    case 5
        Out = fromSIunit_h(h5_pT(p, T5_ps(p, s)));
    otherwise
        Out = NaN;
    end
 
case 'h_px'
    p = toSIunit_p(In1);
    x = toSIunit_x(In2);
    if x > 1 || x < 0 || p >= 22.064 
        Out = NaN;
        return
    end
    hL = h4L_p(p);
    hV = h4V_p(p);
    Out = hL + x * (hV - hL);
 
case 'h_prho'
    p = toSIunit_p(In1);
    rho = 1 / toSIunit_v(1 / In2);
    Region = Region_prho(p, rho);
    switch Region
        case 1
            Out = fromSIunit_h(h1_pT(p, T1_prho(p, rho)));
        case 2
            Out = fromSIunit_h(h2_pT(p, T2_prho(p, rho)));
        case 3
            Out = fromSIunit_h(h3_rhoT(rho, T3_prho(p, rho)));
        case 4
            if p < 16.529
                vV = v2_pT(p, T4_p(p));
                vL = v1_pT(p, T4_p(p));
            else
                vV = v3_ph(p, h4V_p(p));
                vL = v3_ph(p, h4L_p(p));
            end
            hV = h4V_p(p);
            hL = h4L_p(p);
            x = (1 / rho - vL) / (vV - vL);
            Out = fromSIunit_h((1 - x) * hL + x * hV);
        case 5
            Out = fromSIunit_h(h5_pT(p, T5_prho(p, rho)));
        otherwise
            Out = NaN;
    end
 
case 'h_tx'
    T = toSIunit_T(In1);
    x = toSIunit_x(In2);
    if x > 1 || x < 0 || T >= 647.096 
        Out = NaN;
        return
    end
    p = p4_T(T);
    hL = h4L_p(p);
    hV = h4V_p(p);
    Out = hL + x * (hV - hL);
 
    %***********************************************************************************************************
    %*1.5 Specific Volume (v)
case {'vv_p','rhov_p'}
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_v(v2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_v(v3_ph(p, h4V_p(p)));
        end
    else
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
 
case {'vl_p','rhol_p'}
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_v(v1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_v(v3_ph(p, h4L_p(p)));
        end
    else
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
case {'vv_t','rhov_t'}
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_v(v2_pT(p4_T(T), T));
        else
            Out = fromSIunit_v(v3_ph(p4_T(T), h4V_p(p4_T(T))));
        end
    else
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
case {'vl_t','rhol_t'}
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_v(v1_pT(p4_T(T), T));
        else
            Out = fromSIunit_v(v3_ph(p4_T(T), h4L_p(p4_T(T))));
        end
    else
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
case {'v_pt','rho_pt'}
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_v(v1_pT(p, T));
    case 2
        Out = fromSIunit_v(v2_pT(p, T));
    case 3
        Out = fromSIunit_v(v3_ph(p, h3_pT(p, T)));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_v(v5_pT(p, T));
    otherwise
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
 
case {'v_ph','rho_ph'}
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        Out = fromSIunit_v(v1_pT(p, T1_ph(p, h)));
    case 2
        Out = fromSIunit_v(v2_pT(p, T2_ph(p, h)));
    case 3
        Out = fromSIunit_v(v3_ph(p, h));
    case 4
        xs = x4_ph(p, h);
        if p < 16.529 
            v4v = v2_pT(p, T4_p(p));
            v4L = v1_pT(p, T4_p(p));
        else
            v4v = v3_ph(p, h4V_p(p));
            v4L = v3_ph(p, h4L_p(p));
        end
        Out = fromSIunit_v((xs * v4v + (1 - xs) * v4L));
    case 5
        Ts = T5_ph(p, h);
        Out = fromSIunit_v(v5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
case {'v_ps','rho_ps'}
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Out = fromSIunit_v(v1_pT(p, T1_ps(p, s)));
    case 2
        Out = fromSIunit_v(v2_pT(p, T2_ps(p, s)));
    case 3
        Out = fromSIunit_v(v3_ps(p, s));
    case 4
        xs = x4_ps(p, s);
        if p < 16.529 
            v4v = v2_pT(p, T4_p(p));
            v4L = v1_pT(p, T4_p(p));
        else
            v4v = v3_ph(p, h4V_p(p));
            v4L = v3_ph(p, h4L_p(p));
        end
        Out = fromSIunit_v((xs * v4v + (1 - xs) * v4L));
    case 5
        Ts = T5_ps(p, s);
        Out = fromSIunit_v(v5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
    if fun(1)=='r' 
        Out=1/Out;
    end
 
    %***********************************************************************************************************
    %*1.6 Density (rho)
    % Density is calculated as 1/v. Se section 1.5 Volume
    %***********************************************************************************************************
    %*1.7 Specific entropy (s)
case 'sv_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_s(s2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'sl_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_s(s1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'sv_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_s(s2_pT(p4_T(T), T));
        else
            Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'sl_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_s(s1_pT(p4_T(T), T));
        else
            Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 's_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_s(s1_pT(p, T));
    case 2
        Out = fromSIunit_s(s2_pT(p, T));
    case 3
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_s(s3_rhoT(rhos, T));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_s(s5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
case 's_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        T = T1_ph(p, h);
        Out = fromSIunit_s(s1_pT(p, T));
    case 2
        T = T2_ph(p, h);
        Out = fromSIunit_s(s2_pT(p, T));
    case 3
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_s(s3_rhoT(rhos, Ts));
    case 4
        Ts = T4_p(p);
        xs = x4_ph(p, h);
        if p < 16.529 
            s4v = s2_pT(p, Ts);
            s4L = s1_pT(p, Ts);
        else
            v4v = v3_ph(p, h4V_p(p));
            s4v = s3_rhoT(1 / v4v, Ts);
            v4L = v3_ph(p, h4L_p(p));
            s4L = s3_rhoT(1 / v4L, Ts);
        end
        Out = fromSIunit_s((xs * s4v + (1 - xs) * s4L));
    case 5
        T = T5_ph(p, h);
        Out = fromSIunit_s(s5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
 
    %***********************************************************************************************************
    %*1.8 Specific internal energy (u)
case 'uv_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_u(u2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'ul_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_u(u1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'uv_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_u(u2_pT(p4_T(T), T));
        else
            Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'ul_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_u(u1_pT(p4_T(T), T));
        else
            Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'u_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_u(u1_pT(p, T));
    case 2
        Out = fromSIunit_u(u2_pT(p, T));
    case 3
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_u(u3_rhoT(rhos, T));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_u(u5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
case 'u_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        Ts = T1_ph(p, h);
        Out = fromSIunit_u(u1_pT(p, Ts));
    case 2
        Ts = T2_ph(p, h);
        Out = fromSIunit_u(u2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_u(u3_rhoT(rhos, Ts));
    case 4
        Ts = T4_p(p);
        xs = x4_ph(p, h);
        if p < 16.529 
            u4v = u2_pT(p, Ts);
            u4L = u1_pT(p, Ts);
        else
            v4v = v3_ph(p, h4V_p(p));
            u4v = u3_rhoT(1 / v4v, Ts);
            v4L = v3_ph(p, h4L_p(p));
            u4L = u3_rhoT(1 / v4L, Ts);
        end
        Out = fromSIunit_u((xs * u4v + (1 - xs) * u4L));
    case 5
        Ts = T5_ph(p, h);
        Out = fromSIunit_u(u5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
 
case 'u_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Ts = T1_ps(p, s);
        Out = fromSIunit_u(u1_pT(p, Ts));
    case 2
        Ts = T2_ps(p, s);
        Out = fromSIunit_u(u2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_u(u3_rhoT(rhos, Ts));
    case 4
        if p < 16.529 
            uLp = u1_pT(p, T4_p(p));
            uVp = u2_pT(p, T4_p(p));
        else
            uLp = u3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p));
            uVp = u3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p));
        end
        xs = x4_ps(p, s);
        Out = fromSIunit_u((xs * uVp + (1 - xs) * uLp));
    case 5
        Ts = T5_ps(p, s);
        Out = fromSIunit_u(u5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
    %***********************************************************************************************************
    %*1.9 Specific isobaric heat capacity (Cp)
case 'cpv_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_Cp(Cp2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'cpl_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_Cp(Cp1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'cpv_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_Cp(Cp2_pT(p4_T(T), T));
        else
            Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'cpl_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_Cp(Cp1_pT(p4_T(T), T));
        else
            Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'cp_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_Cp(Cp1_pT(p, T));
    case 2
        Out = fromSIunit_Cp(Cp2_pT(p, T));
    case 3
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_Cp(Cp3_rhoT(rhos, T));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_Cp(Cp5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
case 'cp_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        Ts = T1_ph(p, h);
        Out = fromSIunit_Cp(Cp1_pT(p, Ts));
    case 2
        Ts = T2_ph(p, h);
        Out = fromSIunit_Cp(Cp2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_Cp(Cp3_rhoT(rhos, Ts));
    case 4
        Out = NaN;
    case 5
        Ts = T5_ph(p, h);
        Out = fromSIunit_Cp(Cp5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
 
case 'cp_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Ts = T1_ps(p, s);
        Out = fromSIunit_Cp(Cp1_pT(p, Ts));
    case 2
        Ts = T2_ps(p, s);
        Out = fromSIunit_Cp(Cp2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_Cp(Cp3_rhoT(rhos, Ts));
    case 4
        Out = NaN;
    case 5
        Ts = T5_ps(p, s);
        Out = fromSIunit_Cp(Cp5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
 
    %***********************************************************************************************************
    %*1.10 Specific isochoric heat capacity (Cv)
case 'cvv_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_Cv(Cv2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'cvl_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_Cv(Cv1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'cvv_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_Cv(Cv2_pT(p4_T(T), T));
        else
            Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'cvl_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_Cv(Cv1_pT(p4_T(T), T));
        else
            Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'cv_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_Cv(Cv1_pT(p, T));
    case 2
        Out = fromSIunit_Cv(Cv2_pT(p, T));
    case 3
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_Cv(Cv3_rhoT(rhos, T));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_Cv(Cv5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
case 'cv_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        Ts = T1_ph(p, h);
        Out = fromSIunit_Cv(Cv1_pT(p, Ts));
    case 2
        Ts = T2_ph(p, h);
        Out = fromSIunit_Cv(Cv2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_Cv(Cv3_rhoT(rhos, Ts));
    case 4
        Out = NaN;
    case 5
        Ts = T5_ph(p, h);
        Out = fromSIunit_Cv(Cv5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
 
 
case 'cv_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Ts = T1_ps(p, s);
        Out = fromSIunit_Cv(Cv1_pT(p, Ts));
    case 2
        Ts = T2_ps(p, s);
        Out = fromSIunit_Cv(Cv2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_Cv(Cv3_rhoT(rhos, Ts));
    case 4
        Out = NaN;  %(xs * CvVp + (1 - xs) * CvLp) / Cv_scale - Cv_offset
    case 5
        Ts = T5_ps(p, s);
        Out = fromSIunit_Cv(Cv5_pT(p, Ts));
    otherwise
        Out = CVErr(xlErrValue);
    end
 
 
    %***********************************************************************************************************
    %*1.11 Speed of sound
case 'wv_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_w(w2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'wl_p'
    p = toSIunit_p(In1);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            Out = fromSIunit_w(w1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
        end
    else
        Out = NaN;
    end
 
case 'wv_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_w(w2_pT(p4_T(T), T));
        else
            Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'wl_t'
    T = toSIunit_T(In1);
    if T > 273.15  && T < 647.096 
        if T <= 623.15 
            Out = fromSIunit_w(w1_pT(p4_T(T), T));
        else
            Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
        end
    else
        Out = NaN;
    end
 
case 'w_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 1
        Out = fromSIunit_w(w1_pT(p, T));
    case 2
        Out = fromSIunit_w(w2_pT(p, T));
    case 3
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_w(w3_rhoT(rhos, T));
    case 4
        Out = NaN;
    case 5
        Out = fromSIunit_w(w5_pT(p, T));
    otherwise
        Out = NaN;
    end
 
 
case 'w_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case 1
        Ts = T1_ph(p, h);
        Out = fromSIunit_w(w1_pT(p, Ts));
    case 2
        Ts = T2_ph(p, h);
        Out = fromSIunit_w(w2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_w(w3_rhoT(rhos, Ts));
    case 4
        Out = NaN;
    case 5
        Ts = T5_ph(p, h);
        Out = fromSIunit_w(w5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
 
case 'w_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    switch Region
    case 1
        Ts = T1_ps(p, s);
        Out = fromSIunit_w(w1_pT(p, Ts));
    case 2
        Ts = T2_ps(p, s);
        Out = fromSIunit_w(w2_pT(p, Ts));
    case 3
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_w(w3_rhoT(rhos, Ts));
    case 4
        Out = NaN; %(xs * wVp + (1 - xs) * wLp) / w_scale - w_offset
    case 5
        Ts = T5_ps(p, s);
        Out = fromSIunit_w(w5_pT(p, Ts));
    otherwise
        Out = NaN;
    end
    %***********************************************************************************************************
    %*1.12 Viscosity
case 'my_pt'
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    switch Region
    case 4
        Out = NaN;
    case {1, 2, 3, 5}
        Out = fromSIunit_my(my_AllRegions_pT(p, T));
    otherwise
        Out = NaN;
    end
 
case {'my_ph'}
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    switch Region
    case {1, 2, 3, 5}
        Out = fromSIunit_my(my_AllRegions_ph(p, h));
    case {4}
        Out = NaN;  
    otherwise
        Out = NaN;
    end
 
case 'my_ps'
    h = XSteam('h_ps',In1, In2);
    Out = XSteam('my_ph',In1, h);
 
    %***********************************************************************************************************
    %*1.13 Prandtl
    case 'pr_pt'
  Cp = toSIunit_Cp(XSteam('Cp_pT',In1, In2));
  my = toSIunit_my(XSteam('my_pT',In1,In2));
  tc = toSIunit_tc(XSteam('tc_pT',In1,In2));
  Out = Cp * 1000 * my / tc;
 
    case 'pr_ph'
  Cp = toSIunit_Cp(XSteam('Cp_ph',In1, In2));
  my = toSIunit_my(XSteam('my_ph',In1,In2));
  tc = toSIunit_tc(XSteam('tc_ph',In1,In2));
  Out = Cp * 1000 * my / tc;
 
    %***********************************************************************************************************
    %*1.14 Kappa
    %***********************************************************************************************************
    %***********************************************************************************************************
    %*1.15 Surface tension
case 'st_t'
    T = toSIunit_T(In1);
    Out = fromSIunit_st(Surface_Tension_T(T));
 
case 'st_p'
    T = XSteam('Tsat_p',In1);
    T = toSIunit_T(T);
    Out = fromSIunit_st(Surface_Tension_T(T));
 
    %***********************************************************************************************************
    %*1.16 Thermal conductivity
case 'tcl_p'
    T = XSteam('Tsat_p',In1);
    v = XSteam('vL_p',In1);
    p = toSIunit_p(In1);
    T = toSIunit_T(T);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
 
 
case 'tcv_p'
    ps = In1;
    T = XSteam('Tsat_p',ps);
    v = XSteam('vV_p',ps);
    p = toSIunit_p(In1);
    T = toSIunit_T(T);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
 
case 'tcl_t'
    Ts = In1;
    p = XSteam('psat_T',Ts);
    v = XSteam('vL_T',Ts);
    p = toSIunit_p(p);
    T = toSIunit_T(Ts);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
 
case 'tcv_t'
    Ts = In1;
    p = XSteam('psat_T',Ts);
    v = XSteam('vV_T',Ts);
    p = toSIunit_p(p);
    T = toSIunit_T(Ts);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
 
case 'tc_pt'
    Ts = In2;
    ps = In1;
    v = XSteam('v_pT',ps, Ts);
    p = toSIunit_p(ps);
    T = toSIunit_T(Ts);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
 
case 'tc_ph'
    hs = In2;
    ps = In1;
    v = XSteam('v_ph',ps, hs);
    T = XSteam('T_ph',ps, hs);
    p = toSIunit_p(ps);
    T = toSIunit_T(T);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
 
case 'tc_hs'
    hs = In1;
    p = XSteam('p_hs',hs, In2);
    ps = p;
    v = XSteam('v_ph',ps, hs);
    T = XSteam('T_ph',ps, hs);
    p = toSIunit_p(p);
    T = toSIunit_T(T);
    v = toSIunit_v(v);
    rho = 1 / v;
    Out = fromSIunit_tc(tc_ptrho(p, T, rho));
    %***********************************************************************************************************
    %*1.17 Vapour fraction
case 'x_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    if p > 0.000611657  && p < 22.06395 
        Out = fromSIunit_x(x4_ph(p, h));
    else
        Out = NaN;
    end
 
case 'x_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    if p > 0.000611657  && p < 22.06395 
        Out = fromSIunit_x(x4_ps(p, s));
    else
        Out = NaN;
    end
 
    %***********************************************************************************************************
    %*1.18 Vapour Volume Fraction
case 'vx_ph'
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            vL = v1_pT(p, T4_p(p));
            vV = v2_pT(p, T4_p(p));
        else
            vL = v3_ph(p, h4L_p(p));
            vV = v3_ph(p, h4V_p(p));
        end
        xs = x4_ph(p, h);
        Out = fromSIunit_vx((xs * vV / (xs * vV + (1 - xs) * vL)));
    else
        Out = NaN;
    end
 
case 'vx_ps'
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    if p > 0.000611657  && p < 22.06395 
        if p < 16.529 
            vL = v1_pT(p, T4_p(p));
            vV = v2_pT(p, T4_p(p));
        else
            vL = v3_ph(p, h4L_p(p));
            vV = v3_ph(p, h4V_p(p));
        end
        xs = x4_ps(p, s);
        Out = fromSIunit_vx((xs * vV / (xs * vV + (1 - xs) * vL)));
    else
        Out = NaN;
    end
 
 case 'check'
   err=check();
 
otherwise
    error(['Unknown calling function to XSteam, ',fun, ' See help XSteam for valid calling functions']);
end
%***********************************************************************************************************
%*2 IAPWS IF 97 Calling functions                                                                          *
%***********************************************************************************************************
%
%***********************************************************************************************************
%*2.1 Functions for region 1
function v1_pT = v1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_pi = 0;
for i = 1 : 34
    gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
end
v1_pT = R * T / p * Pi * gamma_der_pi / 1000;
 
function h1_pT = h1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_tau = 0;
for i = 1 : 34
    gamma_der_tau = gamma_der_tau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1));
end
h1_pT = R * T * tau * gamma_der_tau;
 
function u1_pT = u1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_tau = 0;
gamma_der_pi = 0;
for i = 1 : 34
    gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
    gamma_der_tau = gamma_der_tau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1));
end
u1_pT = R * T * (tau * gamma_der_tau - Pi * gamma_der_pi);
 
function s1_pT = s1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma = 0;
gamma_der_tau = 0;
for i = 1 : 34
    gamma_der_tau = gamma_der_tau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1));
    gamma = gamma + n1(i) * (7.1 - Pi) ^ I1(i) * (tau - 1.222) ^ J1(i);
end
s1_pT = R * tau * gamma_der_tau - R * gamma;
 
function Cp1_pT = Cp1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_tautau = 0;
for i = 1 :34
    gamma_der_tautau = gamma_der_tautau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2));
end
Cp1_pT = -R * tau ^ 2 * gamma_der_tautau;
 
function Cv1_pT = Cv1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_pi = 0;
gamma_der_pipi = 0;
gamma_der_pitau = 0;
gamma_der_tautau = 0;
for i = 1 : 34
    gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
    gamma_der_pipi = gamma_der_pipi + n1(i) * I1(i) * (I1(i) - 1) * (7.1 - Pi) ^ (I1(i) - 2) * (tau - 1.222) ^ J1(i);
    gamma_der_pitau = gamma_der_pitau - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * J1(i) * (tau - 1.222) ^ (J1(i) - 1);
    gamma_der_tautau = gamma_der_tautau + n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2);
end
Cv1_pT = R * (-tau ^ 2 * gamma_der_tautau + (gamma_der_pi - tau * gamma_der_pitau) ^ 2 / gamma_der_pipi);
 
function w1_pT = w1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_pi = 0;
gamma_der_pipi = 0;
gamma_der_pitau = 0;
gamma_der_tautau = 0;
for i = 1 : 34
    gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
    gamma_der_pipi = gamma_der_pipi + n1(i) * I1(i) * (I1(i) - 1) * (7.1 - Pi) ^ (I1(i) - 2) * (tau - 1.222) ^ J1(i);
    gamma_der_pitau = gamma_der_pitau - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * J1(i) * (tau - 1.222) ^ (J1(i) - 1);
    gamma_der_tautau = gamma_der_tautau + n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2);
end
w1_pT = (1000 * R * T * gamma_der_pi ^ 2 / ((gamma_der_pi - tau * gamma_der_pitau) ^ 2 / (tau ^ 2 * gamma_der_tautau) - gamma_der_pipi)) ^ 0.5;
 
function T1_ph = T1_ph(p, h)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.1 The Backward Equation T ( p,h )
%Eqution 11, Table 6, Page 10
I1 = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6];
J1 = [0, 1, 2, 6, 22, 32, 0, 1, 2, 3, 4, 10, 32, 10, 32, 10, 32, 32, 32, 32];
n1 = [-238.72489924521, 404.21188637945, 113.49746881718, -5.8457616048039, -1.528548241314E-04, -1.0866707695377E-06, -13.391744872602, 43.211039183559, -54.010067170506, 30.535892203916, -6.5964749423638, 9.3965400878363E-03, 1.157364750534E-07, -2.5858641282073E-05, -4.0644363084799E-09, 6.6456186191635E-08, 8.0670734103027E-11, -9.3477771213947E-13, 5.8265442020601E-15, -1.5020185953503E-17];
Pi = p / 1;
eta = h / 2500;
T = 0;
for i = 1 : 20
    T = T + n1(i) * Pi ^ I1(i) * (eta + 1) ^ J1(i);
end
T1_ph = T;
 
function T1_ps = T1_ps(p, s)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.2 The Backward Equation T ( p, s )
%Eqution 13, Table 8, Page 11
I1 = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4];
J1 = [0, 1, 2, 3, 11, 31, 0, 1, 2, 3, 12, 31, 0, 1, 2, 9, 31, 10, 32, 32];
n1 = [174.78268058307, 34.806930892873, 6.5292584978455, 0.33039981775489, -1.9281382923196E-07, -2.4909197244573E-23, -0.26107636489332, 0.22592965981586, -0.064256463395226, 7.8876289270526E-03, 3.5672110607366E-10, 1.7332496994895E-24, 5.6608900654837E-04, -3.2635483139717E-04, 4.4778286690632E-05, -5.1322156908507E-10, -4.2522657042207E-26, 2.6400441360689E-13, 7.8124600459723E-29, -3.0732199903668E-31];
Pi = p / 1;
Sigma = s / 1;
T = 0;
for i = 1 : 20
    T = T + n1(i) * Pi ^ I1(i) * (Sigma + 2) ^ J1(i);
end
T1_ps = T;
 
function p1_hs = p1_hs(h, s)
%Supplementary Release on Backward Equations for Pressure as a Function of Enthalpy and Entropy p(h,s) to the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%5 Backward Equation p(h,s) for Region 1
%Eqution 1, Table 2, Page 5
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 5];
J1 = [0, 1, 2, 4, 5, 6, 8, 14, 0, 1, 4, 6, 0, 1, 10, 4, 1, 4, 0];
n1 = [-0.691997014660582, -18.361254878756, -9.28332409297335, 65.9639569909906, -16.2060388912024, 450.620017338667, 854.68067822417, 6075.23214001162, 32.6487682621856, -26.9408844582931, -319.9478483343, -928.35430704332, 30.3634537455249, -65.0540422444146, -4309.9131651613, -747.512324096068, 730.000345529245, 1142.84032569021, -436.407041874559];
eta = h / 3400;
Sigma = s / 7.6;
p = 0;
for i = 1 : 19
    p = p + n1(i) * (eta + 0.05) ^ I1(i) * (Sigma + 0.05) ^ J1(i);
end
p1_hs = p * 100;
 
function T1_prho = T1_prho(p ,rho)
  %Solve by iteration. Observe that for low temperatures this equation has 2 solutions.
  %Solve with half interval method
  Low_Bound = 273.15;
  High_Bound = T4_p(p);
  rhos=-1000;
  while abs(rho - rhos) > 0.00001
    Ts = (Low_Bound + High_Bound) / 2;
    rhos = 1 / v1_pT(p, Ts);
    if rhos < rho
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
  T1_prho = Ts;
 
%***********************************************************************************************************
%*2.2 functions for region 2
 
function v2_pT = v2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_pi = 1 / Pi;
gr_pi = 0;
for i = 1 : 43
    gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
end
v2_pT = R * T / p * Pi * (g0_pi + gr_pi) / 1000;
 
function h2_pT = h2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tau = 0;
for i = 1 : 9
    g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1);
end
gr_tau = 0;
for i = 1 : 43
    gr_tau = gr_tau + nr(i) * Pi ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
end
h2_pT = R * T * tau * (g0_tau + gr_tau);
 
function u2_pT = u2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_pi = 1 / Pi;
g0_tau = 0;
for i = 1 : 9
    g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1);
end
gr_pi = 0;
gr_tau = 0;
for i = 1 : 43
    gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
    gr_tau = gr_tau + nr(i) * Pi ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
end
u2_pT = R * T * (tau * (g0_tau + gr_tau) - Pi * (g0_pi + gr_pi));
 
 
function s2_pT = s2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0 = log(Pi);
g0_tau = 0;
for i = 1 : 9
    g0 = g0 + n0(i) * tau ^ J0(i);
    g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1);
end
gr = 0;
gr_tau = 0;
for i = 1 : 43
    gr = gr + nr(i) * Pi ^ Ir(i) * (tau - 0.5) ^ Jr(i);
    gr_tau = gr_tau + nr(i) * Pi ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
end
s2_pT = R * (tau * (g0_tau + gr_tau) - (g0 + gr));
 
function Cp2_pT = Cp2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tautau = 0;
for i = 1 : 9
    g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2);
end
gr_tautau = 0;
for i = 1 : 43
    gr_tautau = gr_tautau + nr(i) * Pi ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2);
end
Cp2_pT = -R * tau ^ 2 * (g0_tautau + gr_tautau);
 
function Cv2_pT = Cv2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tautau = 0;
for i =  1 : 9
    g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2);
end
gr_pi = 0;
gr_pitau = 0;
gr_pipi = 0;
gr_tautau = 0;
for i = 1 : 43
    gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
    gr_pipi = gr_pipi + nr(i) * Ir(i) * (Ir(i) - 1) * Pi ^ (Ir(i) - 2) * (tau - 0.5) ^ Jr(i);
    gr_pitau = gr_pitau + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
    gr_tautau = gr_tautau + nr(i) * Pi ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2);
end
Cv2_pT = R * (-tau ^ 2 * (g0_tautau + gr_tautau) - (1 + Pi * gr_pi - tau * Pi * gr_pitau) ^ 2 / (1 - Pi ^ 2 * gr_pipi));
 
function w2_pT = w2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tautau = 0;
for i =  1 : 9
    g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2);
end
gr_pi = 0;
gr_pitau = 0;
gr_pipi = 0;
gr_tautau = 0;
for i = 1 : 43
    gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
    gr_pipi = gr_pipi + nr(i) * Ir(i) * (Ir(i) - 1) * Pi ^ (Ir(i) - 2) * (tau - 0.5) ^ Jr(i);
    gr_pitau = gr_pitau + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
    gr_tautau = gr_tautau + nr(i) * Pi ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2);
end
w2_pT = (1000 * R * T * (1 + 2 * Pi * gr_pi + Pi ^ 2 * gr_pi ^ 2) / ((1 - Pi ^ 2 * gr_pipi) + (1 + Pi * gr_pi - tau * Pi * gr_pitau) ^ 2 / (tau ^ 2 * (g0_tautau + gr_tautau)))) ^ 0.5;
 
function T2_ph = T2_ph(p, h)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2,6.3.1 The Backward Equations T( p, h ) for Subregions 2a, 2b, and 2c
if p < 4 
    sub_reg = 1;
else
    if p < (905.84278514723 - 0.67955786399241 * h + 1.2809002730136E-04 * h ^ 2) 
        sub_reg = 2;
    else
        sub_reg = 3;
    end
end
 
switch sub_reg
case 1
    %Subregion A
    %Table 20, Eq 22, page 22
    Ji = [0, 1, 2, 3, 7, 20, 0, 1, 2, 3, 7, 9, 11, 18, 44, 0, 2, 7, 36, 38, 40, 42, 44, 24, 44, 12, 32, 44, 32, 36, 42, 34, 44, 28];
    Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7];
    ni = [1089.8952318288, 849.51654495535, -107.81748091826, 33.153654801263, -7.4232016790248, 11.765048724356, 1.844574935579, -4.1792700549624, 6.2478196935812, -17.344563108114, -200.58176862096, 271.96065473796, -455.11318285818, 3091.9688604755, 252266.40357872, -6.1707422868339E-03, -0.31078046629583, 11.670873077107, 128127984.04046, -985549096.23276, 2822454697.3002, -3594897141.0703, 1722734991.3197, -13551.334240775, 12848734.66465, 1.3865724283226, 235988.32556514, -13105236.545054, 7399.9835474766, -551966.9703006, 3715408.5996233, 19127.72923966, -415351.64835634, -62.459855192507];
    Ts = 0;
    hs = h / 2000;
    for i = 1 : 34
        Ts = Ts + ni(i) * p ^ (Ii(i)) * (hs - 2.1) ^ Ji(i);
    end
    T2_ph = Ts;
case 2
    %Subregion B
    %Table 21, Eq 23, page 23
    Ji = [0, 1, 2, 12, 18, 24, 28, 40, 0, 2, 6, 12, 18, 24, 28, 40, 2, 8, 18, 40, 1, 2, 12, 24, 2, 12, 18, 24, 28, 40, 18, 24, 40, 28, 2, 28, 1, 40];
    Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 7, 7, 9, 9];
    ni = [1489.5041079516, 743.07798314034, -97.708318797837, 2.4742464705674, -0.63281320016026, 1.1385952129658, -0.47811863648625, 8.5208123431544E-03, 0.93747147377932, 3.3593118604916, 3.3809355601454, 0.16844539671904, 0.73875745236695, -0.47128737436186, 0.15020273139707, -0.002176411421975, -0.021810755324761, -0.10829784403677, -0.046333324635812, 7.1280351959551E-05, 1.1032831789999E-04, 1.8955248387902E-04, 3.0891541160537E-03, 1.3555504554949E-03, 2.8640237477456E-07, -1.0779857357512E-05, -7.6462712454814E-05, 1.4052392818316E-05, -3.1083814331434E-05, -1.0302738212103E-06, 2.821728163504E-07, 1.2704902271945E-06, 7.3803353468292E-08, -1.1030139238909E-08, -8.1456365207833E-14, -2.5180545682962E-11, -1.7565233969407E-18, 8.6934156344163E-15];
    Ts = 0;
    hs = h / 2000;
    for i = 1 : 38
        Ts = Ts + ni(i) * (p - 2) ^ (Ii(i)) * (hs - 2.6) ^ Ji(i);
    end
    T2_ph = Ts;
otherwise
    %Subregion C
    %Table 22, Eq 24, page 24
    Ji = [0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16, 20, 22];
    Ii = [-7, -7, -6, -6, -5, -5, -2, -2, -1, -1, 0, 0, 1, 1, 2, 6, 6, 6, 6, 6, 6, 6, 6];
    ni = [-3236839855524.2, 7326335090218.1, 358250899454.47, -583401318515.9, -10783068217.47, 20825544563.171, 610747.83564516, 859777.2253558, -25745.72360417, 31081.088422714, 1208.2315865936, 482.19755109255, 3.7966001272486, -10.842984880077, -0.04536417267666, 1.4559115658698E-13, 1.126159740723E-12, -1.7804982240686E-11, 1.2324579690832E-07, -1.1606921130984E-06, 2.7846367088554E-05, -5.9270038474176E-04, 1.2918582991878E-03];
    Ts = 0;
    hs = h / 2000;
    for i = 1 : 23
        Ts = Ts + ni(i) * (p + 25) ^ (Ii(i)) * (hs - 1.8) ^ Ji(i);
    end
    T2_ph = Ts;
end
 
function T2_ps = T2_ps(p, s)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2,6.3.2 The Backward Equations T( p, s ) for Subregions 2a, 2b, and 2c
%Page 26
if p < 4 
    sub_reg = 1;
else
    if s < 5.85 
        sub_reg = 3;
    else
        sub_reg = 2;
    end
end
switch sub_reg
case 1
    %Subregion A
    %Table 25, Eq 25, page 26
    Ii = [-1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.25, -1.25, -1.25, -1, -1, -1, -1, -1, -1, -0.75, -0.75, -0.5, -0.5, -0.5, -0.5, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.75, 0.75, 0.75, 0.75, 1, 1, 1.25, 1.25, 1.5, 1.5];
    Ji = [-24, -23, -19, -13, -11, -10, -19, -15, -6, -26, -21, -17, -16, -9, -8, -15, -14, -26, -13, -9, -7, -27, -25, -11, -6, 1, 4, 8, 11, 0, 1, 5, 6, 10, 14, 16, 0, 4, 9, 17, 7, 18, 3, 15, 5, 18];
    ni = [-392359.83861984, 515265.7382727, 40482.443161048, -321.93790923902, 96.961424218694, -22.867846371773, -449429.14124357, -5011.8336020166, 0.35684463560015, 44235.33584819, -13673.388811708, 421632.60207864, 22516.925837475, 474.42144865646, -149.31130797647, -197811.26320452, -23554.39947076, -19070.616302076, 55375.669883164, 3829.3691437363, -603.91860580567, 1936.3102620331, 4266.064369861, -5978.0638872718, -704.01463926862, 338.36784107553, 20.862786635187, 0.033834172656196, -4.3124428414893E-05, 166.53791356412, -139.86292055898, -0.78849547999872, 0.072132411753872, -5.9754839398283E-03, -1.2141358953904E-05, 2.3227096733871E-07, -10.538463566194, 2.0718925496502, -0.072193155260427, 2.074988708112E-07, -0.018340657911379, 2.9036272348696E-07, 0.21037527893619, 2.5681239729999E-04, -0.012799002933781, -8.2198102652018E-06];
    Pi = p;
    Sigma = s / 2;
    teta = 0;
    for i = 1 : 46
        teta = teta + ni(i) * Pi ^ Ii(i) * (Sigma - 2) ^ Ji(i);
    end
    T2_ps = teta;
case 2
    %Subregion B
    %Table 26, Eq 26, page 27
    Ii = [-6, -6, -5, -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5];
    Ji = [0, 11, 0, 11, 0, 1, 11, 0, 1, 11, 12, 0, 1, 6, 10, 0, 1, 5, 8, 9, 0, 1, 2, 4, 5, 6, 9, 0, 1, 2, 3, 7, 8, 0, 1, 5, 0, 1, 3, 0, 1, 0, 1, 2];
    ni = [316876.65083497, 20.864175881858, -398593.99803599, -21.816058518877, 223697.85194242, -2784.1703445817, 9.920743607148, -75197.512299157, 2970.8605951158, -3.4406878548526, 0.38815564249115, 17511.29508575, -1423.7112854449, 1.0943803364167, 0.89971619308495, -3375.9740098958, 471.62885818355, -1.9188241993679, 0.41078580492196, -0.33465378172097, 1387.0034777505, -406.63326195838, 41.72734715961, 2.1932549434532, -1.0320050009077, 0.35882943516703, 5.2511453726066E-03, 12.838916450705, -2.8642437219381, 0.56912683664855, -0.099962954584931, -3.2632037778459E-03, 2.3320922576723E-04, -0.1533480985745, 0.029072288239902, 3.7534702741167E-04, 1.7296691702411E-03, -3.8556050844504E-04, -3.5017712292608E-05, -1.4566393631492E-05, 5.6420857267269E-06, 4.1286150074605E-08, -2.0684671118824E-08, 1.6409393674725E-09];
    Pi = p;
    Sigma = s / 0.7853;
    teta = 0;
    for i = 1 : 44
        teta = teta + ni(i) * Pi ^ Ii(i) * (10 - Sigma) ^ Ji(i);
    end
    T2_ps = teta;
otherwise
    %Subregion C
    %Table 27, Eq 27, page 28
    Ii = [-2, -2, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 7, 7, 7, 7];
    Ji = [0, 1, 0, 0, 1, 2, 3, 0, 1, 3, 4, 0, 1, 2, 0, 1, 5, 0, 1, 4, 0, 1, 2, 0, 1, 0, 1, 3, 4, 5];
    ni = [909.68501005365, 2404.566708842, -591.6232638713, 541.45404128074, -270.98308411192, 979.76525097926, -469.66772959435, 14.399274604723, -19.104204230429, 5.3299167111971, -21.252975375934, -0.3114733441376, 0.60334840894623, -0.042764839702509, 5.8185597255259E-03, -0.014597008284753, 5.6631175631027E-03, -7.6155864584577E-05, 2.2440342919332E-04, -1.2561095013413E-05, 6.3323132660934E-07, -2.0541989675375E-06, 3.6405370390082E-08, -2.9759897789215E-09, 1.0136618529763E-08, 5.9925719692351E-12, -2.0677870105164E-11, -2.0874278181886E-11, 1.0162166825089E-10, -1.6429828281347E-10];
    Pi = p;
    Sigma = s / 2.9251;
    teta = 0;
    for i = 1 : 30
        teta = teta + ni(i) * Pi ^ Ii(i) * (2 - Sigma) ^ Ji(i);
    end
    T2_ps = teta;
end
 
function p2_hs = p2_hs(h, s)
%Supplementary Release on Backward Equations for Pressure as a function of Enthalpy and Entropy p(h,s) to the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%Chapter 6:Backward Equations p(h,s) for Region 2
if h < (-3498.98083432139 + 2575.60716905876 * s - 421.073558227969 * s ^ 2 + 27.6349063799944 * s ^ 3) 
    sub_reg = 1;
else
    if s < 5.85 
        sub_reg = 3;
    else
        sub_reg = 2;
    end
end
switch sub_reg
case 1
    %Subregion A
    %Table 6, Eq 3, page 8
    Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 5, 6, 7];
    Ji = [1, 3, 6, 16, 20, 22, 0, 1, 2, 3, 5, 6, 10, 16, 20, 22, 3, 16, 20, 0, 2, 3, 6, 16, 16, 3, 16, 3, 1];
    ni = [-1.82575361923032E-02, -0.125229548799536, 0.592290437320145, 6.04769706185122, 238.624965444474, -298.639090222922, 0.051225081304075, -0.437266515606486, 0.413336902999504, -5.16468254574773, -5.57014838445711, 12.8555037824478, 11.414410895329, -119.504225652714, -2847.7798596156, 4317.57846408006, 1.1289404080265, 1974.09186206319, 1516.12444706087, 1.41324451421235E-02, 0.585501282219601, -2.97258075863012, 5.94567314847319, -6236.56565798905, 9659.86235133332, 6.81500934948134, -6332.07286824489, -5.5891922446576, 4.00645798472063E-02];
    eta = h / 4200;
    Sigma = s / 12;
    Pi = 0;
    for i = 1 : 29
        Pi = Pi + ni(i) * (eta - 0.5) ^ Ii(i) * (Sigma - 1.2) ^ Ji(i);
    end
    p2_hs = Pi ^ 4 * 4;
case 2
    %Subregion B
    %Table 7, Eq 4, page 9
    Ii = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 8, 8, 12, 14];
    Ji = [0, 1, 2, 4, 8, 0, 1, 2, 3, 5, 12, 1, 6, 18, 0, 1, 7, 12, 1, 16, 1, 12, 1, 8, 18, 1, 16, 1, 3, 14, 18, 10, 16];
    ni = [8.01496989929495E-02, -0.543862807146111, 0.337455597421283, 8.9055545115745, 313.840736431485, 0.797367065977789, -1.2161697355624, 8.72803386937477, -16.9769781757602, -186.552827328416, 95115.9274344237, -18.9168510120494, -4334.0703719484, 543212633.012715, 0.144793408386013, 128.024559637516, -67230.9534071268, 33697238.0095287, -586.63419676272, -22140322476.9889, 1716.06668708389, -570817595.806302, -3121.09693178482, -2078413.8463301, 3056059461577.86, 3221.57004314333, 326810259797.295, -1441.04158934487, 410.694867802691, 109077066873.024, -24796465425889.3, 1888019068.65134, -123651009018773];
    eta = h / 4100;
    Sigma = s / 7.9;
    Pi = 0;
    for i = 1 : 33
        Pi = Pi + ni(i) * (eta - 0.6) ^ Ii(i) * (Sigma - 1.01) ^ Ji(i);
    end
    p2_hs = Pi ^ 4 * 100;
otherwise
    %Subregion C
    %Table 8, Eq 5, page 10
    Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 6, 6, 10, 12, 16];
    Ji = [0, 1, 2, 3, 4, 8, 0, 2, 5, 8, 14, 2, 3, 7, 10, 18, 0, 5, 8, 16, 18, 18, 1, 4, 6, 14, 8, 18, 7, 7, 10];
    ni = [0.112225607199012, -3.39005953606712, -32.0503911730094, -197.5973051049, -407.693861553446, 13294.3775222331, 1.70846839774007, 37.3694198142245, 3581.44365815434, 423014.446424664, -751071025.760063, 52.3446127607898, -228.351290812417, -960652.417056937, -80705929.2526074, 1626980172256.69, 0.772465073604171, 46392.9973837746, -13731788.5134128, 1704703926305.12, -25110462818730.8, 31774883083552, 53.8685623675312, -55308.9094625169, -1028615.22421405, 2042494187562.34, 273918446.626977, -2.63963146312685E+15, -1078908541.08088, -29649262098.0124, -1.11754907323424E+15];
    eta = h / 3500;
    Sigma = s / 5.9;
    Pi = 0;
    for i = 1 : 31
        Pi = Pi + ni(i) * (eta - 0.7) ^ Ii(i) * (Sigma - 1.1) ^ Ji(i);
    end
    p2_hs = Pi ^ 4 * 100;
end
 
function T2_prho=T2_prho(p,rho) 
  %Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  %Solve with half interval method
 
  if p < 16.5292 
    Low_Bound = T4_p(p);
  else
    Low_Bound = B23T_p(p);
  end
  High_Bound = 1073.15;
  rhos=-1000;
  while abs(rho - rhos) > 0.000001
    Ts = (Low_Bound + High_Bound) / 2;
    rhos = 1 / v2_pT(p, Ts);
    if rhos < rho
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
    T2_prho = Ts;
 
%***********************************************************************************************************
%*2.3 functions for region 3
 
function p3_rhoT = p3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fidelta = 0;
for i = 2 : 40
    fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
end
fidelta = fidelta + ni(1) / delta;
p3_rhoT = rho * R * T * delta * fidelta / 1000;
 
function u3_rhoT = u3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitau = 0;
for i = 2 : 40
    fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1);
end
u3_rhoT = R * T * (tau * fitau);
 
 
function h3_rhoT = h3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fidelta = 0;
fitau = 0;
for i = 2 : 40
    fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
    fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1);
end
fidelta = fidelta + ni(1) / delta;
h3_rhoT = R * T * (tau * fitau + delta * fidelta);
 
function s3_rhoT = s3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fi = 0;
fitau = 0;
for i = 2 : 40
    fi = fi + ni(i) * delta ^ Ii(i) * tau ^ Ji(i);
    fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1);
end
fi = fi + ni(1) * log(delta);
s3_rhoT = R * (tau * fitau - fi);
 
function Cp3_rhoT = Cp3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitautau = 0;
fidelta = 0;
fideltatau = 0;
fideltadelta = 0;
for i = 2 : 40
    fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2);
    fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
    fideltatau = fideltatau + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * Ji(i) * tau ^ (Ji(i) - 1);
    fideltadelta = fideltadelta + ni(i) * Ii(i) * (Ii(i) - 1) * delta ^ (Ii(i) - 2) * tau ^ Ji(i);
end
fidelta = fidelta + ni(1) / delta;
fideltadelta = fideltadelta - ni(1) / (delta ^ 2);
Cp3_rhoT = R * (-tau ^ 2 * fitautau + (delta * fidelta - delta * tau * fideltatau) ^ 2 / (2 * delta * fidelta + delta ^ 2 * fideltadelta));
 
function Cv3_rhoT = Cv3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitautau = 0;
for i = 1 : 40
    fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2);
end
Cv3_rhoT = R * -(tau * tau * fitautau);
 
function w3_rhoT = w3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitautau = 0;
fidelta = 0;
fideltatau = 0;
fideltadelta = 0;
for i = 2 : 40
    fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2);
    fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
    fideltatau = fideltatau + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * Ji(i) * tau ^ (Ji(i) - 1);
    fideltadelta = fideltadelta + ni(i) * Ii(i) * (Ii(i) - 1) * delta ^ (Ii(i) - 2) * tau ^ Ji(i);
end
fidelta = fidelta + ni(1) / delta;
fideltadelta = fideltadelta - ni(1) / (delta ^ 2);
w3_rhoT = (1000 * R * T * (2 * delta * fidelta + delta ^ 2 * fideltadelta - (delta * fidelta - delta * tau * fideltatau) ^ 2 / (tau ^ 2 * fitautau))) ^ 0.5;
 
function T3_ph = T3_ph(p, h)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
%Boundary equation, Eq 1 Page 5
h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3;
if h < h3ab 
    %Subregion 3a
    %Eq 2, Table 3, Page 7
    Ii = [-12, -12, -12, -12, -12, -12, -12, -12, -10, -10, -10, -8, -8, -8, -8, -5, -3, -2, -2, -2, -1, -1, 0, 0, 1, 3, 3, 4, 4, 10, 12];
    Ji = [0, 1, 2, 6, 14, 16, 20, 22, 1, 5, 12, 0, 2, 4, 10, 2, 0, 1, 3, 4, 0, 2, 0, 1, 1, 0, 1, 0, 3, 4, 5];
    ni = [-1.33645667811215E-07, 4.55912656802978E-06, -1.46294640700979E-05, 6.3934131297008E-03, 372.783927268847, -7186.54377460447, 573494.7521034, -2675693.29111439, -3.34066283302614E-05, -2.45479214069597E-02, 47.8087847764996, 7.64664131818904E-06, 1.28350627676972E-03, 1.71219081377331E-02, -8.51007304583213, -1.36513461629781E-02, -3.84460997596657E-06, 3.37423807911655E-03, -0.551624873066791, 0.72920227710747, -9.92522757376041E-03, -0.119308831407288, 0.793929190615421, 0.454270731799386, 0.20999859125991, -6.42109823904738E-03, -0.023515586860454, 2.52233108341612E-03, -7.64885133368119E-03, 1.36176427574291E-02, -1.33027883575669E-02];
    ps = p / 100;
    hs = h / 2300;
    Ts = 0;
    for i = 1 : 31
        Ts = Ts + ni(i) * (ps + 0.24) ^ Ii(i) * (hs - 0.615) ^ Ji(i);
    end
    T3_ph = Ts * 760;
else
    %Subregion 3b
    %Eq 3, Table 4, Page 7,8
    Ii = [-12, -12, -10, -10, -10, -10, -10, -8, -8, -8, -8, -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1, -1, -1, -1, -1, 0, 0, 1, 3, 5, 6, 8];
    Ji = [0, 1, 0, 1, 5, 10, 12, 0, 1, 2, 4, 10, 0, 1, 2, 0, 1, 5, 0, 4, 2, 4, 6, 10, 14, 16, 0, 2, 1, 1, 1, 1, 1];
    ni = [3.2325457364492E-05, -1.27575556587181E-04, -4.75851877356068E-04, 1.56183014181602E-03, 0.105724860113781, -85.8514221132534, 724.140095480911, 2.96475810273257E-03, -5.92721983365988E-03, -1.26305422818666E-02, -0.115716196364853, 84.9000969739595, -1.08602260086615E-02, 1.54304475328851E-02, 7.50455441524466E-02, 2.52520973612982E-02, -6.02507901232996E-02, -3.07622221350501, -5.74011959864879E-02, 5.03471360939849, -0.925081888584834, 3.91733882917546, -77.314600713019, 9493.08762098587, -1410437.19679409, 8491662.30819026, 0.861095729446704, 0.32334644281172, 0.873281936020439, -0.436653048526683, 0.286596714529479, -0.131778331276228, 6.76682064330275E-03];
    hs = h / 2800;
    ps = p / 100;
    Ts = 0;
    for i = 1 : 33
        Ts = Ts + ni(i) * (ps + 0.298) ^ Ii(i) * (hs - 0.72) ^ Ji(i);
    end
    T3_ph = Ts * 860;
end
 
function v3_ph = v3_ph(p, h)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
%Boundary equation, Eq 1 Page 5
h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3;
if h < h3ab 
    %Subregion 3a
    %Eq 4, Table 6, Page 9
    Ii = [-12, -12, -12, -12, -10, -10, -10, -8, -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1, -1, -1, 0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 8];
    Ji = [6, 8, 12, 18, 4, 7, 10, 5, 12, 3, 4, 22, 2, 3, 7, 3, 16, 0, 1, 2, 3, 0, 1, 0, 1, 2, 0, 2, 0, 2, 2, 2];
    ni = [5.29944062966028E-03, -0.170099690234461, 11.1323814312927, -2178.98123145125, -5.06061827980875E-04, 0.556495239685324, -9.43672726094016, -0.297856807561527, 93.9353943717186, 1.92944939465981E-02, 0.421740664704763, -3689141.2628233, -7.37566847600639E-03, -0.354753242424366, -1.99768169338727, 1.15456297059049, 5683.6687581596, 8.08169540124668E-03, 0.172416341519307, 1.04270175292927, -0.297691372792847, 0.560394465163593, 0.275234661176914, -0.148347894866012, -6.51142513478515E-02, -2.92468715386302, 6.64876096952665E-02, 3.52335014263844, -1.46340792313332E-02, -2.24503486668184, 1.10533464706142, -4.08757344495612E-02];
    ps = p / 100;
    hs = h / 2100;
    vs = 0;
    for i = 1 : 32
        vs = vs + ni(i) * (ps + 0.128) ^ Ii(i) * (hs - 0.727) ^ Ji(i);
    end
    v3_ph = vs * 0.0028;
else
    %Subregion 3b
    %Eq 5, Table 7, Page 9
    Ii = [-12, -12, -8, -8, -8, -8, -8, -8, -6, -6, -6, -6, -6, -6, -4, -4, -4, -3, -3, -2, -2, -1, -1, -1, -1, 0, 1, 1, 2, 2];
    Ji = [0, 1, 0, 1, 3, 6, 7, 8, 0, 1, 2, 5, 6, 10, 3, 6, 10, 0, 2, 1, 2, 0, 1, 4, 5, 0, 0, 1, 2, 6];
    ni = [-2.25196934336318E-09, 1.40674363313486E-08, 2.3378408528056E-06, -3.31833715229001E-05, 1.07956778514318E-03, -0.271382067378863, 1.07202262490333, -0.853821329075382, -2.15214194340526E-05, 7.6965608822273E-04, -4.31136580433864E-03, 0.453342167309331, -0.507749535873652, -100.475154528389, -0.219201924648793, -3.21087965668917, 607.567815637771, 5.57686450685932E-04, 0.18749904002955, 9.05368030448107E-03, 0.285417173048685, 3.29924030996098E-02, 0.239897419685483, 4.82754995951394, -11.8035753702231, 0.169490044091791, -1.79967222507787E-02, 3.71810116332674E-02, -5.36288335065096E-02, 1.6069710109252];
    ps = p / 100;
    hs = h / 2800;
    vs = 0;
    for i = 1 : 30
        vs = vs + ni(i) * (ps + 0.0661) ^ Ii(i) * (hs - 0.72) ^ Ji(i);
    end
    v3_ph = vs * 0.0088;
end
 
function T3_ps = T3_ps(p, s)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
%Boundary equation, Eq 6 Page 11
if s <= 4.41202148223476 
    %Subregion 3a
    %Eq 6, Table 10, Page 11
    Ii = [-12, -12, -10, -10, -10, -10, -8, -8, -8, -8, -6, -6, -6, -5, -5, -5, -4, -4, -4, -2, -2, -1, -1, 0, 0, 0, 1, 2, 2, 3, 8, 8, 10];
    Ji = [28, 32, 4, 10, 12, 14, 5, 7, 8, 28, 2, 6, 32, 0, 14, 32, 6, 10, 36, 1, 4, 1, 6, 0, 1, 4, 0, 0, 3, 2, 0, 1, 2];
    ni = [1500420082.63875, -159397258480.424, 5.02181140217975E-04, -67.2057767855466, 1450.58545404456, -8238.8953488889, -0.154852214233853, 11.2305046746695, -29.7000213482822, 43856513263.5495, 1.37837838635464E-03, -2.97478527157462, 9717779473494.13, -5.71527767052398E-05, 28830.794977842, -74442828926270.3, 12.8017324848921, -368.275545889071, 6.64768904779177E+15, 0.044935925195888, -4.22897836099655, -0.240614376434179, -4.74341365254924, 0.72409399912611, 0.923874349695897, 3.99043655281015, 3.84066651868009E-02, -3.59344365571848E-03, -0.735196448821653, 0.188367048396131, 1.41064266818704E-04, -2.57418501496337E-03, 1.23220024851555E-03];
    Sigma = s / 4.4;
    Pi = p / 100;
    teta = 0;
    for i = 1 : 33
        teta = teta + ni(i) * (Pi + 0.24) ^ Ii(i) * (Sigma - 0.703) ^ Ji(i);
    end
    T3_ps = teta * 760;
else
    %Subregion 3b
    %Eq 7, Table 11, Page 11
    Ii = [-12, -12, -12, -12, -8, -8, -8, -6, -6, -6, -5, -5, -5, -5, -5, -4, -3, -3, -2, 0, 2, 3, 4, 5, 6, 8, 12, 14];
    Ji = [1, 3, 4, 7, 0, 1, 3, 0, 2, 4, 0, 1, 2, 4, 6, 12, 1, 6, 2, 0, 1, 1, 0, 24, 0, 3, 1, 2];
    ni = [0.52711170160166, -40.1317830052742, 153.020073134484, -2247.99398218827, -0.193993484669048, -1.40467557893768, 42.6799878114024, 0.752810643416743, 22.6657238616417, -622.873556909932, -0.660823667935396, 0.841267087271658, -25.3717501764397, 485.708963532948, 880.531517490555, 2650155.92794626, -0.359287150025783, -656.991567673753, 2.41768149185367, 0.856873461222588, 0.655143675313458, -0.213535213206406, 5.62974957606348E-03, -316955725450471, -6.99997000152457E-04, 1.19845803210767E-02, 1.93848122022095E-05, -2.15095749182309E-05];
    Sigma = s / 5.3;
    Pi = p / 100;
    teta = 0;
    for i = 1 : 28
        teta = teta + ni(i) * (Pi + 0.76) ^ Ii(i) * (Sigma - 0.818) ^ Ji(i);
    end
    T3_ps = teta * 860;
end
 
function v3_ps = v3_ps(p, s)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
%Boundary equation, Eq 6 Page 11
if s <= 4.41202148223476 
    %Subregion 3a
    %Eq 8, Table 13, Page 14
    Ii = [-12, -12, -12, -10, -10, -10, -10, -8, -8, -8, -8, -6, -5, -4, -3, -3, -2, -2, -1, -1, 0, 0, 0, 1, 2, 4, 5, 6];
    Ji = [10, 12, 14, 4, 8, 10, 20, 5, 6, 14, 16, 28, 1, 5, 2, 4, 3, 8, 1, 2, 0, 1, 3, 0, 0, 2, 2, 0];
    ni = [79.5544074093975, -2382.6124298459, 17681.3100617787, -1.10524727080379E-03, -15.3213833655326, 297.544599376982, -35031520.6871242, 0.277513761062119, -0.523964271036888, -148011.182995403, 1600148.99374266, 1708023226634.27, 2.46866996006494E-04, 1.6532608479798, -0.118008384666987, 2.537986423559, 0.965127704669424, -28.2172420532826, 0.203224612353823, 1.10648186063513, 0.52612794845128, 0.277000018736321, 1.08153340501132, -7.44127885357893E-02, 1.64094443541384E-02, -6.80468275301065E-02, 0.025798857610164, -1.45749861944416E-04];
    Pi = p / 100;
    Sigma = s / 4.4;
    omega = 0;
    for i = 1 : 28
        omega = omega + ni(i) * (Pi + 0.187) ^ Ii(i) * (Sigma - 0.755) ^ Ji(i);
    end
    v3_ps = omega * 0.0028;
else
    %Subregion 3b
    %Eq 9, Table 14, Page 14
    Ii = [-12, -12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -5, -5, -5, -4, -4, -4, -4, -3, -2, -2, -2, -2, -2, -2, 0, 0, 0, 1, 1, 2];
    Ji = [0, 1, 2, 3, 5, 6, 0, 1, 2, 4, 0, 1, 2, 3, 0, 1, 2, 3, 1, 0, 1, 2, 3, 4, 12, 0, 1, 2, 0, 2, 2];
    ni = [5.91599780322238E-05, -1.85465997137856E-03, 1.04190510480013E-02, 5.9864730203859E-03, -0.771391189901699, 1.72549765557036, -4.67076079846526E-04, 1.34533823384439E-02, -8.08094336805495E-02, 0.508139374365767, 1.28584643361683E-03, -1.63899353915435, 5.86938199318063, -2.92466667918613, -6.14076301499537E-03, 5.76199014049172, -12.1613320606788, 1.67637540957944, -7.44135838773463, 3.78168091437659E-02, 4.01432203027688, 16.0279837479185, 3.17848779347728, -3.58362310304853, -1159952.60446827, 0.199256573577909, -0.122270624794624, -19.1449143716586, -1.50448002905284E-02, 14.6407900162154, -3.2747778718823];
    Pi = p / 100;
    Sigma = s / 5.3;
    omega = 0;
    for i = 1 : 31
        omega = omega + ni(i) * (Pi + 0.298) ^ Ii(i) * (Sigma - 0.816) ^ Ji(i);
    end
    v3_ps = omega * 0.0088;
end
 
 
function p3_hs = p3_hs(h, s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Equations as a function of h and s for the Region Boundaries, and an Equation
%( ) sat , T hs for Region 4 of the IAPWS Industrial formulation 1997 for the
%Thermodynamic Properties of Water and Steam
%2004
%Section 3 Backward functions p(h,s), T(h,s), and v(h,s) for Region 3
if s < 4.41202148223476 
    %Subregion 3a
    %Eq 1, Table 3, Page 8
    Ii = [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 6, 7, 8, 10, 10, 14, 18, 20, 22, 22, 24, 28, 28, 32, 32];
    Ji = [0, 1, 5, 0, 3, 4, 8, 14, 6, 16, 0, 2, 3, 0, 1, 4, 5, 28, 28, 24, 1, 32, 36, 22, 28, 36, 16, 28, 36, 16, 36, 10, 28];
    ni = [7.70889828326934, -26.0835009128688, 267.416218930389, 17.2221089496844, -293.54233214597, 614.135601882478, -61056.2757725674, -65127225.1118219, 73591.9313521937, -11664650591.4191, 35.5267086434461, -596.144543825955, -475.842430145708, 69.6781965359503, 335.674250377312, 25052.6809130882, 146997.380630766, 5.38069315091534E+19, 1.43619827291346E+21, 3.64985866165994E+19, -2547.41561156775, 2.40120197096563E+27, -3.93847464679496E+29, 1.47073407024852E+24, -4.26391250432059E+31, 1.94509340621077E+38, 6.66212132114896E+23, 7.06777016552858E+33, 1.75563621975576E+41, 1.08408607429124E+28, 7.30872705175151E+43, 1.5914584739887E+24, 3.77121605943324E+40];
    Sigma = s / 4.4;
    eta = h / 2300;
    Pi = 0;
    for i = 1 : 33
        Pi = Pi + ni(i) * (eta - 1.01) ^ Ii(i) * (Sigma - 0.75) ^ Ji(i);
    end
    p3_hs = Pi * 99;
else
    %Subregion 3b
    %Eq 2, Table 4, Page 8
    Ii = [-12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -8, -6, -6, -6, -6, -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -1, 0, 2, 2, 5, 6, 8, 10, 14, 14];
    Ji = [2, 10, 12, 14, 20, 2, 10, 14, 18, 2, 8, 2, 6, 7, 8, 10, 4, 5, 8, 1, 3, 5, 6, 0, 1, 0, 3, 0, 1, 0, 1, 1, 1, 3, 7];
    ni = [1.25244360717979E-13, -1.26599322553713E-02, 5.06878030140626, 31.7847171154202, -391041.161399932, -9.75733406392044E-11, -18.6312419488279, 510.973543414101, 373847.005822362, 2.99804024666572E-08, 20.0544393820342, -4.98030487662829E-06, -10.230180636003, 55.2819126990325, -206.211367510878, -7940.12232324823, 7.82248472028153, -58.6544326902468, 3550.73647696481, -1.15303107290162E-04, -1.75092403171802, 257.98168774816, -727.048374179467, 1.21644822609198E-04, 3.93137871762692E-02, 7.04181005909296E-03, -82.910820069811, -0.26517881813125, 13.7531682453991, -52.2394090753046, 2405.56298941048, -22736.1631268929, 89074.6343932567, -23923456.5822486, 5687958081.29714];
    Sigma = s / 5.3;
    eta = h / 2800;
    Pi = 0;
    for i = 1 : 35
        Pi = Pi + ni(i) * (eta - 0.681) ^ Ii(i) * (Sigma - 0.792) ^ Ji(i);
    end
    p3_hs = 16.6 / Pi;
end
 
function h3_pT = h3_pT(p, T)
%Not avalible with if 97
%Solve function T3_ph-T=0 with half interval method.
%ver2.6 Start corrected bug
if p < 22.06395    %Bellow tripple point
  Ts = T4_p(p);    %Saturation temperature
  if T <= Ts      %Liquid side
    High_Bound = h4L_p(p); %Max h är liauid h.
    Low_Bound = h1_pT(p, 623.15);
  else
    Low_Bound = h4V_p(p);  %Min h är Vapour h.
    High_Bound = h2_pT(p, B23T_p(p));
  end 
else                  %Above tripple point. R3 from R2 till R3.
  Low_Bound = h1_pT(p, 623.15);
  High_Bound = h2_pT(p, B23T_p(p));
end
%ver2.6 End corrected bug
Ts = T+1;
while abs(T - Ts) > 0.00001
    hs = (Low_Bound + High_Bound) / 2;
    Ts = T3_ph(p, hs);
    if Ts > T 
        High_Bound = hs;
    else
        Low_Bound = hs;
    end
end
h3_pT = hs;
 
function T3_prho = T3_prho(p, rho)
  %Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  %Solve with half interval method
  Low_Bound = 623.15;
  High_Bound = 1073.15;
  ps=-1000;
  while abs(p - ps) > 0.00000001
    Ts = (Low_Bound + High_Bound) / 2;
    ps = p3_rhoT(rho, Ts);
    if ps > p
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
  T3_prho = Ts;
 
 
 
%***********************************************************************************************************
%*2.4 functions for region 4
function p4_T = p4_T(T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Section 8.1 The Saturation-Pressure Equation
%Eq 30, Page 33
teta = T - 0.23855557567849 / (T - 650.17534844798);
a = teta ^ 2 + 1167.0521452767 * teta - 724213.16703206;
B = -17.073846940092 * teta ^ 2 + 12020.82470247 * teta - 3232555.0322333;
C = 14.91510861353 * teta ^ 2 - 4823.2657361591 * teta + 405113.40542057;
p4_T = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4;
 
 
function T4_p = T4_p(p)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Section 8.2 The Saturation-Temperature Equation
%Eq 31, Page 34
beta = p ^ 0.25;
E = beta ^ 2 - 17.073846940092 * beta + 14.91510861353;
f = 1167.0521452767 * beta ^ 2 + 12020.82470247 * beta - 4823.2657361591;
G = -724213.16703206 * beta ^ 2 - 3232555.0322333 * beta + 405113.40542057;
D = 2 * G / (-f - (f ^ 2 - 4 * E * G) ^ 0.5);
T4_p = (650.17534844798 + D - ((650.17534844798 + D) ^ 2 - 4 * (-0.23855557567849 + 650.17534844798 * D)) ^ 0.5) / 2;
 
function h4_s = h4_s(s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,Equations as a function of h and s for the Region Boundaries, and an Equation( ) sat , T hs for Region 4 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%4 Equations for Region Boundaries Given Enthalpy and Entropy
% Se picture page 14
if (s > -0.0001545495919 & s <= 3.77828134)==1
    %hL1_s
    %Eq 3,Table 9,Page 16
    Ii = [0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 7, 8, 12, 12, 14, 14, 16, 20, 20, 22, 24, 28, 32, 32];
    Ji = [14, 36, 3, 16, 0, 5, 4, 36, 4, 16, 24, 18, 24, 1, 4, 2, 4, 1, 22, 10, 12, 28, 8, 3, 0, 6, 8];
    ni = [0.332171191705237, 6.11217706323496E-04, -8.82092478906822, -0.45562819254325, -2.63483840850452E-05, -22.3949661148062, -4.28398660164013, -0.616679338856916, -14.682303110404, 284.523138727299, -113.398503195444, 1156.71380760859, 395.551267359325, -1.54891257229285, 19.4486637751291, -3.57915139457043, -3.35369414148819, -0.66442679633246, 32332.1885383934, 3317.66744667084, -22350.1257931087, 5739538.75852936, 173.226193407919, -3.63968822121321E-02, 8.34596332878346E-07, 5.03611916682674, 65.5444787064505];
    Sigma = s / 3.8;
    eta = 0;
    for i = 1 : 27
        eta = eta + ni(i) * (Sigma - 1.09) ^ Ii(i) * (Sigma + 0.0000366) ^ Ji(i);
    end
    h4_s = eta * 1700;
elseif (s > 3.77828134 & s <= 4.41202148223476 )==1
    %hL3_s
    %Eq 4,Table 10,Page 16
    Ii = [0, 0, 0, 0, 2, 3, 4, 4, 5, 5, 6, 7, 7, 7, 10, 10, 10, 32, 32];
    Ji = [1, 4, 10, 16, 1, 36, 3, 16, 20, 36, 4, 2, 28, 32, 14, 32, 36, 0, 6];
    ni = [0.822673364673336, 0.181977213534479, -0.011200026031362, -7.46778287048033E-04, -0.179046263257381, 4.24220110836657E-02, -0.341355823438768, -2.09881740853565, -8.22477343323596, -4.99684082076008, 0.191413958471069, 5.81062241093136E-02, -1655.05498701029, 1588.70443421201, -85.0623535172818, -31771.4386511207, -94589.0406632871, -1.3927384708869E-06, 0.63105253224098];
    Sigma = s / 3.8;
    eta = 0;
    for i = 1 : 19
        eta = eta + ni(i) * (Sigma - 1.09) ^ Ii(i) * (Sigma + 0.0000366) ^ Ji(i);
    end
    h4_s = eta * 1700;
elseif (s > 4.41202148223476 & s <= 5.85 )==1
    %Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
    %Page 19, Eq 5
    %hV2c3b_s(s)
    Ii = [0, 0, 0, 1, 1, 5, 6, 7, 8, 8, 12, 16, 22, 22, 24, 36];
    Ji = [0, 3, 4, 0, 12, 36, 12, 16, 2, 20, 32, 36, 2, 32, 7, 20];
    ni = [1.04351280732769, -2.27807912708513, 1.80535256723202, 0.420440834792042, -105721.24483466, 4.36911607493884E+24, -328032702839.753, -6.7868676080427E+15, 7439.57464645363, -3.56896445355761E+19, 1.67590585186801E+31, -3.55028625419105E+37, 396611982166.538, -4.14716268484468E+40, 3.59080103867382E+18, -1.16994334851995E+40];
    Sigma = s / 5.9;
    eta = 0;
    for i = 1 : 16
        eta = eta + ni(i) * (Sigma - 1.02) ^ Ii(i) * (Sigma - 0.726) ^ Ji(i);
    end
    h4_s = eta ^ 4 * 2800;
elseif (s > 5.85 & s < 9.155759395)==1
    %Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
    %Page 20, Eq 6
    Ii = [1, 1, 2, 2, 4, 4, 7, 8, 8, 10, 12, 12, 18, 20, 24, 28, 28, 28, 28, 28, 32, 32, 32, 32, 32, 36, 36, 36, 36, 36];
    Ji = [8, 24, 4, 32, 1, 2, 7, 5, 12, 1, 0, 7, 10, 12, 32, 8, 12, 20, 22, 24, 2, 7, 12, 14, 24, 10, 12, 20, 22, 28];
    ni = [-524.581170928788, -9269472.18142218, -237.385107491666, 21077015581.2776, -23.9494562010986, 221.802480294197, -5104725.33393438, 1249813.96109147, 2000084369.96201, -815.158509791035, -157.612685637523, -11420042233.2791, 6.62364680776872E+15, -2.27622818296144E+18, -1.71048081348406E+31, 6.60788766938091E+15, 1.66320055886021E+22, -2.18003784381501E+29, -7.87276140295618E+29, 1.51062329700346E+31, 7957321.70300541, 1.31957647355347E+15, -3.2509706829914E+23, -4.18600611419248E+25, 2.97478906557467E+34, -9.53588761745473E+19, 1.66957699620939E+24, -1.75407764869978E+32, 3.47581490626396E+34, -7.10971318427851E+38];
    Sigma1 = s / 5.21;
    Sigma2 = s / 9.2;
    eta = 0;
    for i = 1 : 30
        eta = eta + ni(i) * (1 / Sigma1 - 0.513) ^ Ii(i) * (Sigma2 - 0.524) ^ Ji(i);
    end
    h4_s = exp(eta) * 2800;
else
    h4_s = -99999;
end
 
function p4_s = p4_s(s)
%Uses h4_s and p_hs for the diffrent regions to determine p4_s
h_sat = h4_s(s);
if (s > -0.0001545495919 & s <= 3.77828134)==1
    p4_s = p1_hs(h_sat, s);
elseif (s > 3.77828134 & s <= 5.210887663)==1
    p4_s = p3_hs(h_sat, s);
elseif (s > 5.210887663 & s < 9.155759395)==1
    p4_s = p2_hs(h_sat, s);
else
    p4_s = -99999;
end
 
function h4L_p = h4L_p(p)
if (p > 0.000611657 & p < 22.06395)==1
    Ts = T4_p(p);
    if p < 16.529 
        h4L_p = h1_pT(p, Ts);
    else
        %Iterate to find the the backward solution of p3sat_h
        Low_Bound = 1670.858218;
        High_Bound = 2087.23500164864;
        ps=-1000;
        while abs(p - ps) > 0.00001
            hs = (Low_Bound + High_Bound) / 2;
            ps = p3sat_h(hs);
            if ps > p 
                High_Bound = hs;
            else
                Low_Bound = hs;
            end
        end
 
        h4L_p = hs;
    end
else
    h4L_p = -99999;
end
 
function h4V_p = h4V_p(p)
if (p > 0.000611657 & p < 22.06395)==1
    Ts = T4_p(p);
    if p < 16.529 
        h4V_p = h2_pT(p, Ts);
    else
        %Iterate to find the the backward solution of p3sat_h
        Low_Bound = 2087.23500164864;
        High_Bound = 2563.592004+5;
        ps=-1000;
        while abs(p - ps) > 0.000001
            hs = (Low_Bound + High_Bound) / 2;
            ps = p3sat_h(hs);
            if ps < p 
                High_Bound = hs;
            else
                Low_Bound = hs;
            end
        end
        h4V_p = hs;
    end
else
    h4V_p = -99999;
end
 
function x4_ph = x4_ph(p, h)
%Calculate vapour fraction from hL and hV for given p
h4v = h4V_p(p);
h4L = h4L_p(p);
if h > h4v 
    x4_ph = 1;
elseif h < h4L 
    x4_ph = 0;
else
    x4_ph = (h - h4L) / (h4v - h4L);
end
 
function x4_ps = x4_ps(p, s)
if p < 16.529 
    ssv = s2_pT(p, T4_p(p));
    ssL = s1_pT(p, T4_p(p));
else
    ssv = s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p));
    ssL = s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p));
end
if s < ssL 
    x4_ps = 0;
elseif s > ssv 
    x4_ps = 1;
else
    x4_ps = (s - ssL) / (ssv - ssL);
end
 
function T4_hs = T4_hs(h, s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Chapter 5.3 page 30.
%The if 97 function is only valid for part of region4. Use iteration outsida.
Ii = [0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 8, 10, 10, 12, 14, 14, 16, 16, 18, 18, 18, 20, 28];
Ji = [0, 3, 12, 0, 1, 2, 5, 0, 5, 8, 0, 2, 3, 4, 0, 1, 1, 2, 4, 16, 6, 8, 22, 1, 20, 36, 24, 1, 28, 12, 32, 14, 22, 36, 24, 36];
ni = [0.179882673606601, -0.267507455199603, 1.162767226126, 0.147545428713616, -0.512871635973248, 0.421333567697984, 0.56374952218987, 0.429274443819153, -3.3570455214214, 10.8890916499278, -0.248483390456012, 0.30415322190639, -0.494819763939905, 1.07551674933261, 7.33888415457688E-02, 1.40170545411085E-02, -0.106110975998808, 1.68324361811875E-02, 1.25028363714877, 1013.16840309509, -1.51791558000712, 52.4277865990866, 23049.5545563912, 2.49459806365456E-02, 2107964.67412137, 366836848.613065, -144814105.365163, -1.7927637300359E-03, 4899556021.00459, 471.262212070518, -82929439019.8652, -1715.45662263191, 3557776.82973575, 586062760258.436, -12988763.5078195, 31724744937.1057];
if (s > 5.210887825 & s < 9.15546555571324)==1
    Sigma = s / 9.2;
    eta = h / 2800;
    teta = 0;
    for i = 1 : 36
        teta = teta + ni(i) * (eta - 0.119) ^ Ii(i) * (Sigma - 1.07) ^ Ji(i);
    end
    T4_hs = teta * 550;
else
    %function psat_h
    if (s > -0.0001545495919 && s <= 3.77828134)==1
        Low_Bound = 0.000611;
        High_Bound = 165.291642526045;
        hL=-1000;
        while abs(hL - h) > 0.00001 && abs(High_Bound - Low_Bound) > 0.0001
            PL = (Low_Bound + High_Bound) / 2;
            Ts = T4_p(PL);
            hL = h1_pT(PL, Ts);
            if hL > h 
                High_Bound = PL;
            else
                Low_Bound = PL;
            end
        end
    elseif s > 3.77828134 && s <= 4.41202148223476 
        PL = p3sat_h(h);
    elseif s > 4.41202148223476 && s <= 5.210887663 
        PL = p3sat_h(h);
    end
    Low_Bound = 0.000611;
    High_Bound = PL;
    sss=-1000;
    while (abs(s - sss) > 0.000001 & abs(High_Bound - Low_Bound) > 0.0000001)==1
        p = (Low_Bound + High_Bound) / 2;
        %Calculate s4_ph
        Ts = T4_p(p);
        xs = x4_ph(p, h);
        if p < 16.529 
            s4v = s2_pT(p, Ts);
            s4L = s1_pT(p, Ts);
        else
            v4v = v3_ph(p, h4V_p(p));
            s4v = s3_rhoT(1 / v4v, Ts);
            v4L = v3_ph(p, h4L_p(p));
            s4L = s3_rhoT(1 / v4L, Ts);
        end
        sss = (xs * s4v + (1 - xs) * s4L);
 
        if sss < s 
            High_Bound = p;
        else
            Low_Bound = p;
        end
    end
    T4_hs = T4_p(p);
end
%***********************************************************************************************************
%*2.5 functions for region 5
function h5_pT = h5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tau = 0;
for i = 1 : 6
    gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1);
end
gammar_tau = 0;
for i = 1 : 5
    gammar_tau = gammar_tau + nir(i) * Pi ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1);
end
h5_pT = R * T * tau * (gamma0_tau + gammar_tau);
 
 
 
function v5_pT = v5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_pi = 1 / Pi;
gammar_pi = 0;
for i = 1 : 5
    gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
end
v5_pT = R * T / p * Pi * (gamma0_pi + gammar_pi) / 1000;
 
 
function u5_pT = u5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_pi = 1 / Pi;
gamma0_tau = 0;
for i = 1 : 6
    gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1);
end
gammar_pi = 0;
gammar_tau = 0;
for i = 1 : 5
    gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
    gammar_tau = gammar_tau + nir(i) * Pi ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1);
end
u5_pT = R * T * (tau * (gamma0_tau + gammar_tau) - Pi * (gamma0_pi + gammar_pi));
 
function Cp5_pT = Cp5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tautau = 0;
for i = 1 : 6
    gamma0_tautau = gamma0_tautau + ni0(i) * Ji0(i) * (Ji0(i) - 1) * tau ^ (Ji0(i) - 2);
end
gammar_tautau = 0;
for i = 1 : 5
    gammar_tautau = gammar_tautau + nir(i) * Pi ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2);
end
Cp5_pT = -R * tau ^ 2 * (gamma0_tautau + gammar_tautau);
 
 
function s5_pT = s5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0 = log(Pi);
gamma0_tau = 0;
for i = 1 : 6
    gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1);
    gamma0 = gamma0 + ni0(i) * tau ^ Ji0(i);
end
gammar = 0;
gammar_tau = 0;
for i = 1 : 5
    gammar = gammar + nir(i) * Pi ^ Iir(i) * tau ^ Jir(i);
    gammar_tau = gammar_tau + nir(i) * Pi ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1);
end
s5_pT = R * (tau * (gamma0_tau + gammar_tau) - (gamma0 + gammar));
 
function Cv5_pT = Cv5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tautau = 0;
for i = 1 : 6
    gamma0_tautau = gamma0_tautau + ni0(i) * (Ji0(i) - 1) * Ji0(i) * tau ^ (Ji0(i) - 2);
end
gammar_pi = 0;
gammar_pitau = 0;
gammar_pipi = 0;
gammar_tautau = 0;
for i = 1 : 5
    gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
    gammar_pitau = gammar_pitau + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * Jir(i) * tau ^ (Jir(i) - 1);
    gammar_pipi = gammar_pipi + nir(i) * Iir(i) * (Iir(i) - 1) * Pi ^ (Iir(i) - 2) * tau ^ Jir(i);
    gammar_tautau = gammar_tautau + nir(i) * Pi ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2);
end
Cv5_pT = R * (-(tau ^ 2 *(gamma0_tautau + gammar_tautau)) - (1 + Pi * gammar_pi - tau * Pi * gammar_pitau)^2 / (1 - Pi ^ 2 * gammar_pipi));
 
function w5_pT = w5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tautau = 0;
for i = 1 : 6
    gamma0_tautau = gamma0_tautau + ni0(i) * (Ji0(i) - 1) * Ji0(i) * tau ^ (Ji0(i) - 2);
end
gammar_pi = 0;
gammar_pitau = 0;
gammar_pipi = 0;
gammar_tautau = 0;
for i = 1 : 5
    gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
    gammar_pitau = gammar_pitau + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * Jir(i) * tau ^ (Jir(i) - 1);
    gammar_pipi = gammar_pipi + nir(i) * Iir(i) * (Iir(i) - 1) * Pi ^ (Iir(i) - 2) * tau ^ Jir(i);
    gammar_tautau = gammar_tautau + nir(i) * Pi ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2);
end
w5_pT = (1000 * R * T * (1 + 2 * Pi * gammar_pi + Pi ^ 2 * gammar_pi ^ 2) / ((1 - Pi ^ 2 * gammar_pipi) + (1 + Pi * gammar_pi - tau * Pi * gammar_pitau) ^ 2 / (tau ^ 2 * (gamma0_tautau + gammar_tautau)))) ^ 0.5;
 
 
function T5_ph = T5_ph(p, h)
%Solve with half interval method
Low_Bound = 1073.15;
High_Bound = 2273.15;
hs=h-1;
while abs(h - hs) > 0.00001
    Ts = (Low_Bound + High_Bound) / 2;
    hs = h5_pT(p, Ts);
    if hs > h 
        High_Bound = Ts;
    else
        Low_Bound = Ts;
    end
end
T5_ph = Ts;
 
 
function T5_ps = T5_ps(p, s)
%Solve with half interval method
Low_Bound = 1073.15;
High_Bound = 2273.15;
ss=s-1;
while abs(s - ss) > 0.00001
    Ts = (Low_Bound + High_Bound) / 2;
    ss = s5_pT(p, Ts);
    if ss > s 
        High_Bound = Ts;
    else
        Low_Bound = Ts;
    end
end
T5_ps = Ts;
 
function T5_prho=T5_prho(p,rho)
  %Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  %Solve with half interval method
 
    Low_Bound = 1073.15;
    High_Bound = 2073.15;
    rhos=-1000;
  while abs(rho - rhos) > 0.000001
    Ts = (Low_Bound + High_Bound) / 2;
    rhos = 1 / v2_pT(p, Ts);
    if rhos < rho
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
    T5_prho = Ts;
 
%***********************************************************************************************************
%*3 Region Selection
%***********************************************************************************************************
%*3.1 Regions as a function of pT
function region_pT = region_pT(p, T)
if T > 1073.15 && p < 10 && T < 2273.15 && p > 0.000611 
    region_pT = 5;
elseif T <= 1073.15 && T > 273.15 && p <= 100 && p > 0.000611 
    if T > 623.15 
        if p > B23p_T(T) 
            region_pT = 3;
            if T < 647.096 
                ps = p4_T(T);
                if abs(p - ps) < 0.00001
                    region_pT = 4;
                end
            end
        else
            region_pT = 2;
        end
    else
        ps = p4_T(T);
        if abs(p - ps) < 0.00001 
            region_pT = 4;
        elseif p > ps
            region_pT = 1;
        else
            region_pT = 2;
        end
    end
else
    region_pT = 0; %**Error, Outside valid area
end
%***********************************************************************************************************
%*3.2 Regions as a function of ph
function region_ph = region_ph(  p,   h)
%Check if outside pressure limits
if p < 0.000611657 || p > 100 
    region_ph = 0;
    return
end
 
%Check if outside low h.
if h < 0.963 * p + 2.2  %Linear adaption to h1_pt()+2 to speed up calcualations.
    if h < h1_pT(p, 273.15) 
        region_ph = 0;
        return
    end
end
 
if p < 16.5292  %Bellow region 3,Check  region 1,4,2,5
    %Check Region 1
    Ts = T4_p(p);
    hL = 109.6635 * log(p) + 40.3481 * p + 734.58; %Approximate function for hL_p
    if abs(h - hL) < 100  %if approximate is not god enough use real function
        hL = h1_pT(p, Ts);
    end
    if h <= hL 
        region_ph = 1;
        return
    end
    %Check Region 4
    hV = 45.1768 * log(p) - 20.158 * p + 2804.4; %Approximate function for hV_p
    if abs(h - hV) < 50  %if approximate is not god enough use real function
        hV = h2_pT(p, Ts);
    end
    if h < hV 
        region_ph = 4;
        return
    end
    %Check upper limit of region 2 Quick Test
    if h < 4000 
        region_ph = 2;
        return
    end
    %Check region 2 (Real value)
    h_45 = h2_pT(p, 1073.15);
    if h <= h_45 
        region_ph = 2;
        return
    end
    %Check region 5
    if p > 10 
        region_ph = 0;
        return
    end
    h_5u = h5_pT(p, 2273.15);
    if h < h_5u 
        region_ph = 5;
        return
    end
    region_ph = 0;
    return
else %for p>16.5292
    %Check if in region1
    if h < h1_pT(p, 623.15) 
        region_ph = 1;
        return
    end
    %Check if in region 3 or 4 (Bellow Reg 2)
    if h < h2_pT(p, B23T_p(p)) 
        %Region 3 or 4
        if p > p3sat_h(h) 
            region_ph = 3;
            return
        else
            region_ph = 4;
            return
        end
    end
    %Check if region 2
    if h < h2_pT(p, 1073.15) 
        region_ph = 2;
        return
    end
end
region_ph = 0;
 
 
%***********************************************************************************************************
%*3.3 Regions as a function of ps
function region_ps = region_ps(  p,   s)
if p < 0.000611657 || p > 100 || s < 0 || s > s5_pT(p, 2273.15) 
    region_ps = 0;
    return
end
 
%Check region 5
if s > s2_pT(p, 1073.15) 
    if p <= 10 
        region_ps = 5;
        return
    else
        region_ps = 0;
        return
    end
end
 
%Check region 2
if p > 16.529 
    ss = s2_pT(p, B23T_p(p)); %Between 5.047  & 5.261. Use to speed up!
else
    ss = s2_pT(p, T4_p(p));
end
if s > ss 
    region_ps = 2;
    return
end
 
%Check region 3
ss = s1_pT(p, 623.15);
if p > 16.529  && s > ss 
    if p > p3sat_s(s) 
        region_ps = 3;
        return
    else
        region_ps = 4;
        return
    end
end
 
%Check region 4 (Not inside region 3)
if p < 16.529  && s > s1_pT(p, T4_p(p)) 
    region_ps = 4;
    return
end
 
%Check region 1
if p > 0.000611657  && s > s1_pT(p, 273.15) 
    region_ps = 1;
    return
end
region_ps = 1;
 
%***********************************************************************************************************
%*3.4 Regions as a function of hs
function region_hs = region_hs(  h,   s)
if s < -0.0001545495919 
    region_hs = 0;
    return
end
%Check linear adaption to p=0.000611. if bellow region 4.
hMin = (((-0.0415878 - 2500.89262) / (-0.00015455 - 9.155759)) * s);
if s < 9.155759395  && h < hMin 
    region_hs = 0;
    return
end
 
%******Kolla 1 eller 4. (+liten bit över B13)
if s >= -0.0001545495919  && s <= 3.77828134 
    if h < h4_s(s) 
        region_hs = 4;
        return
    elseif s < 3.397782955  %100MPa line is limiting
        TMax = T1_ps(100, s);
        hMax = h1_pT(100, TMax);
        if h < hMax 
            region_hs = 1;
            return
        else
            region_hs = 0;
            return
        end
    else %The point is either in region 4,1,3. Check B23
        hB = hB13_s(s);
        if h < hB 
            region_hs = 1;
            return
        end
        TMax = T3_ps(100, s);
        vmax = v3_ps(100, s);
        hMax = h3_rhoT(1 / vmax, TMax);
        if h < hMax 
            region_hs = 3;
            return
        else
            region_hs = 0;
            return
        end
    end
end
 
%******Kolla region 2 eller 4. (Övre delen av område b23-> max)
if s >= 5.260578707  && s <= 11.9212156897728 
    if s > 9.155759395  %Above region 4
        Tmin = T2_ps(0.000611, s);
        hMin = h2_pT(0.000611, Tmin);
        %function adapted to h(1073.15,s)
        hMax = -0.07554022 * s ^ 4 + 3.341571 * s ^ 3 - 55.42151 * s ^ 2 + 408.515 * s + 3031.338;
        if h > hMin  && h < hMax 
            region_hs = 2;
            return
        else
            region_hs = 0;
            return
        end
    end
 
    hV = h4_s(s);
 
    if h < hV   %Region 4. Under region 3.
        region_hs = 4;
        return
    end
    if s < 6.04048367171238 
        TMax = T2_ps(100, s);
        hMax = h2_pT(100, TMax);
    else
        %function adapted to h(1073.15,s)
        hMax = -2.988734 * s ^ 4 + 121.4015 * s ^ 3 - 1805.15 * s ^ 2 + 11720.16 * s - 23998.33;
    end
    if h < hMax   %Region 2. Över region 4.
        region_hs = 2;
        return
    else
        region_hs = 0;
        return
    end
end
 
%Kolla region 3 eller 4. Under kritiska punkten.
if s >= 3.77828134  && s <= 4.41202148223476 
    hL = h4_s(s);
    if h < hL 
        region_hs = 4;
        return
    end
    TMax = T3_ps(100, s);
    vmax = v3_ps(100, s);
    hMax = h3_rhoT(1 / vmax, TMax);
    if h < hMax 
        region_hs = 3;
        return
    else
        region_hs = 0;
        return
    end
end
 
%Kolla region 3 eller 4 från kritiska punkten till övre delen av b23
if s >= 4.41202148223476  && s <= 5.260578707 
    hV = h4_s(s);
    if h < hV 
        region_hs = 4;
        return
    end
    %Kolla om vi är under b23 giltighetsområde.
    if s <= 5.048096828 
        TMax = T3_ps(100, s);
        vmax = v3_ps(100, s);
        hMax = h3_rhoT(1 / vmax, TMax);
        if h < hMax 
            region_hs = 3;
            return
        else
            region_hs = 0;
            return
        end
    else %Inom området för B23 i s led.
        if (h > 2812.942061)  %Ovanför B23 i h_led
            if s > 5.09796573397125 
                TMax = T2_ps(100, s);
                hMax = h2_pT(100, TMax);
                if h < hMax 
                    region_hs = 2;
                    return
                else
                    region_hs = 0;
                    return
                end
            else
                region_hs = 0;
                return
            end
        end
        if (h < 2563.592004)    %Nedanför B23 i h_led men vi har redan kollat ovanför hV2c3b
            region_hs = 3;
            return
        end
        %Vi är inom b23 området i både s och h led.
        Tact = TB23_hs(h, s);
        pact = p2_hs(h, s);
        pBound = B23p_T(Tact);
        if pact > pBound 
            region_hs = 3;
            return
        else
            region_hs = 2;
            return
        end
    end
end
region_hs = 0;
 
%***********************************************************************************************************
%*3.5 Regions as a function of p and rho
function Region_prho = Region_prho(p,rho)
v = 1 / rho;
if p < 0.000611657 || p > 100
    Region_prho = 0;
    return
end
if p < 16.5292  %Bellow region 3, Check region 1,4,2
    if v < v1_pT(p, 273.15) %Observe that this is not actually min of v. Not valid Water of 4°C is ligther.
        Region_prho = 0;
        return
    end
    if v <= v1_pT(p, T4_p(p))
        Region_prho = 1;
        return
    end
    if v < v2_pT(p, T4_p(p)) 
        Region_prho = 4;
        return
    end
    if v <= v2_pT(p, 1073.15)
        Region_prho = 2;
        return
    end 
    if p > 10 %Above region 5
        Region_prho = 0;
        return
    end 
    if v <= v5_pT(p, 2073.15) 
        Region_prho = 5;
        return
    end
else %Check region 1,3,4,3,2 (Above the lowest point of region 3.)
    if v < v1_pT(p, 273.15) %Observe that this is not actually min of v. Not valid Water of 4°C is ligther.
        Region_prho = 0;
        return
    end
    if v < v1_pT(p, 623.15)
        Region_prho = 1;
        return
    end
    %Check if in region 3 or 4 (Bellow Reg 2)
    if v < v2_pT(p, B23T_p(p))
        %Region 3 or 4
        if p > 22.064 %Above region 4
            Region_prho = 3;
            return
        end
        if v < v3_ph(p, h4L_p(p)) || v > v3_ph(p, h4V_p(p)) %Uses iteration!!
            Region_prho = 3;
            return
        else
            Region_prho = 4;
            return
        end 
    end 
    %Check if region 2
    if v < v2_pT(p, 1073.15)
        Region_prho = 2;
        return
    end
end 
Region_prho = 0;
 
 
 
%***********************************************************************************************************
%*4 Region Borders
%***********************************************************************************************************
%***********************************************************************************************************
%*4.1 Boundary between region 2 and 3.
function B23p_T = B23p_T(T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%1997
%Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
%Eq 5, Page 5
B23p_T = 348.05185628969 - 1.1671859879975 * T + 1.0192970039326E-03 * T ^ 2;
 
function B23T_p = B23T_p(p)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%1997
%Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
%Eq 6, Page 6
B23T_p = 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ^ 0.5;
 
 
%***********************************************************************************************************
%*4.2 Region 3. pSat_h  & pSat_s
function p3sat_h = p3sat_h(h)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h)  & T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water  & Steam
%2004
%Section 4 Boundary Equations psat(h)  & psat(s) for the Saturation Lines of Region 3
%Se pictures Page 17, Eq 10, Table 17, Page 18
Ii = [0, 1, 1, 1, 1, 5, 7, 8, 14, 20, 22, 24, 28, 36];
Ji = [0, 1, 3, 4, 36, 3, 0, 24, 16, 16, 3, 18, 8, 24];
ni = [0.600073641753024, -9.36203654849857, 24.6590798594147, -107.014222858224, -91582131580576.8, -8623.32011700662, -23.5837344740032, 2.52304969384128E+17, -3.89718771997719E+18, -3.33775713645296E+22, 35649946963.6328, -1.48547544720641E+26, 3.30611514838798E+18, 8.13641294467829E+37];
hs = h / 2600;
ps = 0;
for i = 1:14
    ps = ps + ni(i) * (hs - 1.02) ^ Ii(i) * (hs - 0.608) ^ Ji(i);
end
p3sat_h = ps * 22;
 
 
function p3sat_s = p3sat_s(s)
Ii = [0, 1, 1, 4, 12, 12, 16, 24, 28, 32];
Ji = [0, 1, 32, 7, 4, 14, 36, 10, 0, 18];
ni = [0.639767553612785, -12.9727445396014, -2.24595125848403E+15, 1774667.41801846, 7170793495.71538, -3.78829107169011E+17, -9.55586736431328E+34, 1.87269814676188E+23, 119254746466.473, 1.10649277244882E+36];
Sigma = s / 5.2;
Pi = 0;
for i = 1:10
    Pi = Pi + ni(i) * (Sigma - 1.03) ^ Ii(i) * (Sigma - 0.699) ^ Ji(i);
end
p3sat_s = Pi * 22;
 
 
%***********************************************************************************************************
%4.3 Region boundary 1to3  & 3to2 as a functions of s
function hB13_s = hB13_s(s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Chapter 4.5 page 23.
Ii = [0, 1, 1, 3, 5, 6];
Ji = [0, -2, 2, -12, -4, -3];
ni = [0.913965547600543, -4.30944856041991E-05, 60.3235694765419, 1.17518273082168E-18, 0.220000904781292, -69.0815545851641];
Sigma = s / 3.8;
eta = 0;
for i = 1 : 6
    eta = eta + ni(i) * (Sigma - 0.884) ^ Ii(i) * (Sigma - 0.864) ^ Ji(i);
end
hB13_s = eta * 1700;
 
 
function TB23_hs = TB23_hs(h, s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Chapter 4.6 page 25.
Ii = [-12, -10, -8, -4, -3, -2, -2, -2, -2, 0, 1, 1, 1, 3, 3, 5, 6, 6, 8, 8, 8, 12, 12, 14, 14];
Ji = [10, 8, 3, 4, 3, -6, 2, 3, 4, 0, -3, -2, 10, -2, -1, -5, -6, -3, -8, -2, -1, -12, -1, -12, 1];
ni = [6.2909626082981E-04, -8.23453502583165E-04, 5.15446951519474E-08, -1.17565945784945, 3.48519684726192, -5.07837382408313E-12, -2.84637670005479, -2.36092263939673, 6.01492324973779, 1.48039650824546, 3.60075182221907E-04, -1.26700045009952E-02, -1221843.32521413, 0.149276502463272, 0.698733471798484, -2.52207040114321E-02, 1.47151930985213E-02, -1.08618917681849, -9.36875039816322E-04, 81.9877897570217, -182.041861521835, 2.61907376402688E-06, -29162.6417025961, 1.40660774926165E-05, 7832370.62349385];
Sigma = s / 5.3;
eta = h / 3000;
teta = 0;
for i = 1 : 25
    teta = teta + ni(i) * (eta - 0.727) ^ Ii(i) * (Sigma - 0.864) ^ Ji(i);
end
TB23_hs = teta * 900;
 
 
%***********************************************************************************************************
%*5 Transport properties
%***********************************************************************************************************
%*5.1 Viscosity (IAPWS formulation 1985, Revised 2003)
%***********************************************************************************************************
function my_AllRegions_pT = my_AllRegions_pT(  p,   T)
h0 = [0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447];
h1 = [0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0];
h2 = [-0.2818107, -1.070786, -1.263184, 0, 0, 0];
h3 = [0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0];
h4 = [-0.0417661, 0, 0, 0.1600435, 0, 0];
h5 = [0, -0.01578386, 0, 0, 0, 0];
h6 = [0, 0, 0, -0.003629481, 0, 0];
 
%Calcualte density.
switch region_pT(p, T)
case 1
    rho = 1 / v1_pT(p, T);
case 2
    rho = 1 / v2_pT(p, T);
case 3
    hs = h3_pT(p, T);
    rho = 1 / v3_ph(p, hs);
case 4
    rho = NaN;
case 5
    rho = 1 / v5_pT(p, T);
otherwise
    my_AllRegions_pT = NaN;
    return
end
 
rhos = rho / 317.763;
Ts = T / 647.226;
ps = p / 22.115;
 
%Check valid area
if T > 900 + 273.15 || (T > 600 + 273.15  && p > 300) || (T > 150 + 273.15  && p > 350) || p > 500 
    my_AllRegions_pT = NaN;
    return
end
my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3));
Sum = 0;
for i = 0 : 5
    Sum = Sum + h0(i+1) * (1 / Ts - 1) ^ i + h1(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6;
end
my1 = exp(rhos * Sum);
mys = my0 * my1;
my_AllRegions_pT = mys * 0.000055071;
 
 
function my_AllRegions_ph = my_AllRegions_ph(  p,   h)
h0 = [0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447];
h1 = [0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0];
h2 = [-0.2818107, -1.070786, -1.263184, 0, 0, 0];
h3 = [0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0];
h4 = [-0.0417661, 0, 0, 0.1600435, 0, 0];
h5 = [0, -0.01578386, 0, 0, 0, 0];
h6 = [0, 0, 0, -0.003629481, 0, 0];
 
%Calcualte density.
switch region_ph(p, h)
case 1
    Ts = T1_ph(p, h);
    T = Ts;
    rho = 1 / v1_pT(p, Ts);
case 2
    Ts = T2_ph(p, h);
    T = Ts;
    rho = 1 / v2_pT(p, Ts);
case 3
    rho = 1 / v3_ph(p, h);
    T = T3_ph(p, h);
case 4
    xs = x4_ph(p, h);
    if p < 16.529 
        v4v = v2_pT(p, T4_p(p));
        v4L = v1_pT(p, T4_p(p));
    else
        v4v = v3_ph(p, h4V_p(p));
        v4L = v3_ph(p, h4L_p(p));
    end
    rho = 1 / (xs * v4v + (1 - xs) * v4L);
    T = T4_p(p);
case 5
    Ts = T5_ph(p, h);
    T = Ts;
    rho = 1 / v5_pT(p, Ts);
otherwise
    my_AllRegions_ph = NaN;
    return
end
rhos = rho / 317.763;
Ts = T / 647.226;
ps = p / 22.115;
%Check valid area
if T > 900 + 273.15 || (T > 600 + 273.15  && p > 300) || (T > 150 + 273.15  && p > 350) || p > 500 
    my_AllRegions_ph = NaN;
    return
end
my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3));
 
Sum = 0;
for i = 0 : 5
    Sum = Sum + h0(i+1) * (1 / Ts - 1) ^ i + h1(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6;
end
my1 = exp(rhos * Sum);
mys = my0 * my1;
my_AllRegions_ph = mys * 0.000055071;
 
%***********************************************************************************************************
%*5.2 Thermal Conductivity (IAPWS formulation 1985)
function tc_ptrho = tc_ptrho(  p,   T,   rho)
%Revised release on the IAPWS formulation 1985 for the Thermal Conductivity of ordinary water
%IAPWS September 1998
%Page 8
%ver2.6 Start corrected bug
 if T < 273.15
   tc_ptrho = NaN; %Out of range of validity (para. B4)
   return
 elseif T < 500 + 273.15
   if p > 100
     tc_ptrho = NaN; %Out of range of validity (para. B4)
     return
   end
 elseif T <= 650 + 273.15
   if p > 70
     tc_ptrho = NaN; %Out of range of validity (para. B4)
     return
   end
 else T <= 800 + 273.15
   if p > 40
      tc_ptrho = NaN; %Out of range of validity (para. B4)
      return
   end
 end
%ver2.6 End corrected bug
T = T / 647.26;
rho = rho / 317.7;
tc0 = T ^ 0.5 * (0.0102811 + 0.0299621 * T + 0.0156146 * T ^ 2 - 0.00422464 * T ^ 3);
tc1 = -0.39707 + 0.400302 * rho + 1.06 * exp(-0.171587 * (rho + 2.39219) ^ 2);
dT = abs(T - 1) + 0.00308976;
Q = 2 + 0.0822994 / dT ^ (3 / 5);
if T >= 1 
    s = 1 / dT;
else
    s = 10.0932 / dT ^ (3 / 5);
end
tc2 = (0.0701309 / T ^ 10 + 0.011852) * rho ^ (9 / 5) * exp(0.642857 * (1 - rho ^ (14 / 5))) + 0.00169937 * s * rho ^ Q * exp((Q / (1 + Q)) * (1 - rho ^ (1 + Q))) - 1.02 * exp(-4.11717 * T ^ (3 / 2) - 6.17937 / rho ^ 5);
tc_ptrho = tc0 + tc1 + tc2;
 
%***********************************************************************************************************
%5.3 Surface Tension
function Surface_Tension_T = Surface_Tension_T(  T)
%IAPWS Release on Surface Tension of Ordinary Water Substance,
%September 1994
tc = 647.096; %K
B = 0.2358;    %N/m
bb = -0.625;
my = 1.256;
if T < 0.01 || T > tc 
    Surface_Tension_T = NaN; %"Out of valid region"
    return
end
tau = 1 - T / tc;
Surface_Tension_T = B * tau ^ my * (1 + bb * tau);
 
 
%***********************************************************************************************************
%*6 Units                                                                                      *
%***********************************************************************************************************
function toSIunit_p = toSIunit_p( Ins )
%Translate bar to MPa
toSIunit_p = Ins / 10;
function fromSIunit_p = fromSIunit_p( Ins )
%Translate bar to MPa
fromSIunit_p = Ins * 10;
function toSIunit_T = toSIunit_T( Ins )
%Translate degC to Kelvon
toSIunit_T = Ins + 273.15;
function fromSIunit_T = fromSIunit_T( Ins )
%Translate Kelvin to degC
fromSIunit_T = Ins - 273.15;
function toSIunit_h = toSIunit_h( Ins )
toSIunit_h = Ins;
function  fromSIunit_h = fromSIunit_h( Ins )
fromSIunit_h = Ins;
function toSIunit_v = toSIunit_v( Ins )
toSIunit_v = Ins;
function fromSIunit_v = fromSIunit_v( Ins )
fromSIunit_v = Ins;
function toSIunit_s = toSIunit_s( Ins )
toSIunit_s = Ins;
function fromSIunit_s = fromSIunit_s( Ins )
fromSIunit_s = Ins;
function toSIunit_u = toSIunit_u( Ins )
toSIunit_u = Ins;
function fromSIunit_u = fromSIunit_u( Ins )
fromSIunit_u = Ins;
function toSIunit_Cp = toSIunit_Cp( Ins )
toSIunit_Cp = Ins;
function fromSIunit_Cp = fromSIunit_Cp( Ins )
fromSIunit_Cp = Ins;
function toSIunit_Cv = toSIunit_Cv( Ins )
toSIunit_Cv = Ins;
function fromSIunit_Cv = fromSIunit_Cv( Ins )
fromSIunit_Cv = Ins;
function toSIunit_w = toSIunit_w( Ins )
toSIunit_w = Ins;
function fromSIunit_w = fromSIunit_w( Ins )
fromSIunit_w = Ins;
function toSIunit_tc = toSIunit_tc( Ins )
toSIunit_tc = Ins;
function fromSIunit_tc = fromSIunit_tc( Ins )
fromSIunit_tc = Ins;
function toSIunit_st = toSIunit_st( Ins )
toSIunit_st = Ins;
function fromSIunit_st = fromSIunit_st( Ins )
fromSIunit_st = Ins;
function toSIunit_x = toSIunit_x( Ins )
toSIunit_x = Ins;
function fromSIunit_x = fromSIunit_x( Ins )
fromSIunit_x = Ins;
function toSIunit_vx = toSIunit_vx( Ins )
toSIunit_vx = Ins;
function fromSIunit_vx = fromSIunit_vx( Ins )
fromSIunit_vx = Ins;
function toSIunit_my = toSIunit_my( Ins )
toSIunit_my = Ins;
function fromSIunit_my = fromSIunit_my( Ins )
fromSIunit_my = Ins;
 
%***********************************************************************************************************
%*7 Verification                                                                                      *
%***********************************************************************************************************
 
function err = check()
err=0;
%*********************************************************************************************************
%* 7.1 Verifiy region 1
%IF-97 Table 5, Page 9
p=[30/10,800/10,30/10];
T=[300,300,500];
Fun={'v1_pT','h1_pT','u1_pT','s1_pT','Cp1_pT','w1_pT'};
IF97=[0.00100215168,0.000971180894,0.001202418;...
        115.331273,184.142828,975.542239;...
        112.324818,106.448356,971.934985;...
        0.392294792,0.368563852,2.58041912;...
        4.17301218,4.01008987,4.65580682;...
        1507.73921,1634.69054,1240.71337];
R1=zeros(6,3);
for i=1:3
    for j=1:6
        R1(j,i)=eval([char(Fun(j)),'(',num2str(p(i)),',',num2str(T(i)),');']);
    end
end
Region1_error=abs((R1-IF97)./IF97)
err=err+sum(sum(Region1_error>1E-8))
 
%IF-97 Table 7, Page 11
p=[30/10,800/10,800/10];
h=[500,500,1500];
IF97=[391.798509,378.108626,611.041229];
R1=zeros(1,3);
for i=1:3
    R1(i)=T1_ph(p(i),h(i));
end
T1_ph_error=abs((R1-IF97)./IF97)
err=err+sum(sum(T1_ph_error>1E-8))
 
%IF-97 Table 9, Page 12
p=[30/10,800/10,800/10];
s=[0.5,0.5,3];
IF97=[307.842258,309.979785,565.899909];
R1=zeros(1,3);
for i=1:3
    R1(i)=T1_ps(p(i),s(i));
end
T1_ps_error=abs((R1-IF97)./IF97)
err=err+sum(sum(T1_ps_error>1E-8))
 
%Supplementary Release on Backward Equations 
%for Pressure as a Function of Enthalpy and Entropy p(h,s) 
%Table 3, Page 6
h=[0.001,90,1500];
s=[0,0,3.4];
IF97=[0.0009800980612,91.929547272,58.68294423];
R1=zeros(1,3);
for i=1:3
    R1(i)=p1_hs(h(i),s(i));
end
p1_hs_error=abs((R1-IF97)./IF97)
err=err+sum(sum(p1_hs_error>1E-8))
 
%*********************************************************************************************************
%* 7.2 Verifiy region 2
% IF-97 Table 15, Page 17
p=[0.035/10,0.035/10,300/10];
T=[300,700,700];
Fun={'v2_pT','h2_pT','u2_pT','s2_pT','Cp2_pT','w2_pT'};
IF97=[39.4913866,92.3015898,0.00542946619;...
        2549.91145,3335.68375,2631.49474;...
        2411.6916,3012.62819,2468.61076;...
        8.52238967,10.1749996,5.17540298;...
        1.91300162,2.08141274,10.3505092;...
        427.920172,644.289068,480.386523];
R2=zeros(6,3);
for i=1:3
    for j=1:6
        R2(j,i)=eval([char(Fun(j)),'(',num2str(p(i)),',',num2str(T(i)),');']);
    end
end
Region2_error=abs((R2-IF97)./IF97)
err=err+sum(sum(Region2_error>1E-8))
 
%IF-97 Table 24, Page 25
p=[0.01/10,30/10,30/10,50/10,50/10,250/10,400/10,600/10,600/10];
h=[3000,3000,4000,3500,4000,3500,2700,2700,3200];
IF97=[534.433241,575.37337,1010.77577,801.299102,1015.31583,875.279054,743.056411,791.137067,882.75686];
R2=zeros(1,9);
for i=1:9
    R2(i)=T2_ph(p(i),h(i));
end
T2_ph_error=abs((R2-IF97)./IF97)
err=err+sum(sum(T2_ph_error>1E-8))
 
%IF-97 Table 29, Page 29
p=[1/10,1/10,25/10,80/10,80/10,900/10,200/10,800/10,800/10];
s=[7.5,8,8,6,7.5,6,5.75,5.25,5.75];
IF97=[399.517097,514.127081,1039.84917,600.48404,1064.95556,1038.01126,697.992849,854.011484,949.017998];
R2=zeros(1,9);
for i=1:9
    R2(i)=T2_ps(p(i),s(i));
end
T2_ps_error=abs((R2-IF97)./IF97)
err=err+sum(sum(T2_ps_error>1E-8))
 
%Supplementary Release on Backward Equations for Pressure as a Function of Enthalpy and Entropy p(h,s) 
%Table 3, Page 6
h=[2800,2800,4100,2800,3600,3600,2800,2800,3400];
s=[6.5,9.5,9.5,6,6,7,5.1,5.8,5.8];
IF97=[1.371012767,0.001879743844,0.1024788997,4.793911442,83.95519209,7.527161441,94.3920206,8.414574124,83.76903879];
R2=zeros(1,9);
for i=1:9
    R2(i)=p2_hs(h(i),s(i));
end
p2_hs_error=abs((R2-IF97)./IF97)
err=err+sum(sum(p2_hs_error>1E-8))
 
%*********************************************************************************************************
%* 7.3 Verifiy region 3
% IF-97 Table 33, Page 32
T=[650,650,750];
rho=[500,200,500];
Fun={'p3_rhoT','h3_rhoT','u3_rhoT','s3_rhoT','Cp3_rhoT','w3_rhoT'};
IF97=[25.5837018,22.2930643,78.3095639;...
        1863.43019,2375.12401,2258.68845;...
        1812.26279,2263.65868,2102.06932;...
        4.05427273,4.85438792,4.46971906;...
        13.8935717,44.6579342,6.34165359;...
        502.005554,383.444594,760.696041];
R3=zeros(6,3);
for i=1:3
    for j=1:6
        R3(j,i)=eval([char(Fun(j)),'(',num2str(rho(i)),',',num2str(T(i)),');']);
    end
end
Region3_error=abs((R3-IF97)./IF97)
err=err+sum(sum(Region3_error>1E-8))
 
%T3_ph
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
h=[1700,2000,2100,2500,2400,2700];
IF97=[629.3083892,690.5718338,733.6163014,641.8418053,735.1848618,842.0460876];
R3=zeros(1,6);
for i=1:6
    R3(i)=T3_ph(p(i),h(i));
end
T3_ph_error=abs((R3-IF97)./IF97)
err=err+sum(sum(T3_ph_error>1E-8))
 
%v3_ph
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
h=[1700,2000,2100,2500,2400,2700];
IF97=[0.001749903962,0.001908139035,0.001676229776,0.006670547043,0.0028012445,0.002404234998];
R3=zeros(1,6);
for i=1:6
    R3(i)=v3_ph(p(i),h(i));
end
v3_ph_error=abs((R3-IF97)./IF97)
err=err+sum(sum(v3_ph_error>1E-7))
 
%T3_ps
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
s=[3.7,3.5,4,5,4.5,5];
IF97=[620.8841563,618.1549029,705.6880237,640.1176443,716.3687517,847.4332825];
R3=zeros(1,6);
for i=1:6
    R3(i)=T3_ps(p(i),s(i));
end
T3_ps_error=abs((R3-IF97)./IF97)
err=err+sum(sum(T3_ps_error>1E-8))
 
%v3_ps
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
s=[3.7,3.5,4,5,4.5,5];
IF97=[0.001639890984,0.001423030205,0.001555893131,0.006262101987,0.002332634294,0.002449610757];
R3=zeros(1,6);
for i=1:6
    R3(i)=v3_ps(p(i),s(i));
end
v3_ps_error=abs((R3-IF97)./IF97)
err=err+sum(sum(v3_ps_error>1E-8))
 
%p3_hs
h=[1700,2000,2100,2500,2400,2700];
s=[3.8,4.2,4.3,5.1,4.7,5];
IF97=[25.55703246,45.40873468,60.7812334,17.20612413,63.63924887,88.39043281];
R3=zeros(1,6);
for i=1:6
    R3(i)=p3_hs(h(i),s(i));
end
p3_hs_error=abs((R3-IF97)./IF97)
err=err+sum(sum(p3_hs_error>1E-8))
 
%h3_pT (Iteration)
p=[255.83702,222.93064,783.09564]./10;
T=[650,650,750];
IF97=[1863.271389,2375.696155,2258.626582];
R3=zeros(1,3);
for i=1:3
    R3(i)=h3_pT(p(i),T(i));
end
h3_pT_error=abs((R3-IF97)./IF97)
err=err+sum(sum(h3_pT_error>1E-6)) %Decimals in IF97
 
%*********************************************************************************************************
%* 7.4 Verifiy region 4
%Saturation pressure, If97, Table 35, Page 34
T=[300,500,600];
IF97=[0.0353658941, 26.3889776, 123.443146]/10;
R3=zeros(1,3);
for i=1:3
    R4(i)=p4_T(T(i));
end
p4_t_error=abs((R4-IF97)./IF97)
err=err+sum(sum( p4_t_error>1E-7))
 
T=[1,10,100]/10;
IF97=[372.755919,453.035632,584.149488];
R3=zeros(1,3);
for i=1:3
    R4(i)=T4_p(T(i));
end
T4_p_error=abs((R4-IF97)./IF97)
err=err+sum(sum( T4_p_error>1E-7))
 
s=[1,2,3,3.8,4,4.2,7,8,9,5.5,5,4.5];
IF97=[308.5509647,700.6304472,1198.359754,1685.025565,1816.891476,1949.352563,2723.729985,2599.04721,2511.861477,2687.69385,2451.623609,2144.360448];
R3=zeros(1,12);
for i=1:12
    R4(i)=h4_s(s(i));
end
h4_s_error=abs((R4-IF97)./IF97)
err=err+sum(sum( h4_s_error>1E-7)) 
 
%*********************************************************************************************************
%* 7.5 Verifiy region 5
% IF-97 Table 42, Page 39
T=[1500,1500,2000];
p=[5,80,80]/10;
Fun={'v5_pT','h5_pT','u5_pT','s5_pT','Cp5_pT','w5_pT'};
IF97=[1.38455354,0.0865156616,0.115743146;...
        5219.76332,5206.09634,6583.80291;...
        4527.48654,4513.97105,5657.85774;...
        9.65408431,8.36546724,9.15671044;...
        2.61610228,2.64453866,2.8530675;...
        917.071933,919.708859,1054.35806];
R5=zeros(6,3);
for i=1:3
    for j=1:6
        R5(j,i)=eval([char(Fun(j)),'(',num2str(p(i)),',',num2str(T(i)),');']);
    end
end
Region5_error=abs((R5-IF97)./IF97)
err=err+sum(sum(Region5_error>1E-8))
 
%T5_ph (Iteration)
p=[0.5,8,8];
h=[5219.76331549428,5206.09634477373,6583.80290533381];
IF97=[1500,1500,2000];
R5=zeros(1,3);
for i=1:3
    R5(i)=T5_ph(p(i),h(i));
end
T5_ph_error=abs((R5-IF97)./IF97)
err=err+sum(sum(T5_ph_error>1E-7)) %Decimals in IF97
 
%T5_ps (Iteration)
p=[0.5,8,8];
s=[9.65408430982588,8.36546724495503,9.15671044273249];
IF97=[1500,1500,2000];
R5=zeros(1,3);
for i=1:3
    R5(i)=T5_ps(p(i),s(i));
end
T5_ps_error=abs((R5-IF97)./IF97)
err=err+sum(sum(T5_ps_error>1E-4)) %Decimals in IF97
 
%*********************************************************************************************************
%* 7.6 Verifiy calling funtions
 
fun={'Tsat_p','T_ph','T_ps','T_hs','psat_T','p_hs','hV_p','hL_p','hV_T','hL_T','h_pT','h_ps','h_px','h_prho','h_Tx','vV_p','vL_p','vV_T','vL_T','v_pT','v_ph','v_ps','rhoV_p','rhoL_p','rhoV_T','rhoL_T','rho_pT','rho_ph','rho_ps','sV_p','sL_p','sV_T','sL_T','s_pT','s_ph','uV_p','uL_p','uV_T','uL_T','u_pT','u_ph','u_ps','CpV_p','CpL_p','CpV_T','CpL_T','Cp_pT','Cp_ph','Cp_ps','CvV_p','CvL_p','CvV_T','CvL_T','Cv_pT','Cv_ph','Cv_ps','wV_p','wL_p','wV_T','wL_T','w_pT','w_ph','w_ps','my_pT','my_ph','my_ps','tcL_p','tcV_p','tcL_T','tcV_T','tc_pT','tc_ph','tc_hs','st_T','st_p','x_ph','x_ps','vx_ph','vx_ps'};
In1={'1','1','1','100','100','84','1','1','100','100','1','1','1','1','100','1','1','100','100','1','1','1','1','1','100','100','1','1','1','0.006117','0.0061171','0.0001','100','1','1','1','1','100','100','1','1','1','1','1','100','100','1','1','1','1','1','100','100','1','1','1','1','1','100','100','1','1','1','1','1','1','1','1','25','25','1','1','100','100','1','1','1','1','1'};
In2={'0','100','1','0.2','0','0.296','0','0','0','0','20','1','0.5','2','0.5','0','0','0','0','100','1000','5','0','0','0','0','100','1000','1','0','0','0','0','20','84.01181117','0','0','0','0','100','1000','1','0','0','0','0','100','200','1','0','0','0','0','100','200','1','0','0','0','0','100','200','1','100','100','1','0','0','0','0','25','100','0.34','0','0','1000','4','418','4'};
Control={'99.60591861','23.84481908','73.70859421','13.84933511','1.014179779','2.295498269','2674.949641','417.4364858','2675.572029','419.099155','84.01181117','308.6107171','1546.193063','1082.773391','1547.33559210927','1.694022523','0.001043148','1.671860601','0.001043455','1.695959407','0.437925658','1.03463539','0.590310924','958.6368897','0.598135993','958.3542773','0.589636754','2.283492601','975.6236788','9.155465556','1.8359E-05','9.155756716','1.307014328','0.296482921','0.296813845','2505.547389','417.332171','2506.015308','418.9933299','2506.171426','956.2074342','308.5082185','2.075938025','4.216149431','2.077491868','4.216645119','2.074108555','4.17913573168802','4.190607038','1.552696979','3.769699683','1.553698696','3.76770022','1.551397249','4.035176364','3.902919468','472.0541571','1545.451948','472.2559492','1545.092249','472.3375235','1542.682475','1557.8585','1.22704E-05','0.000914003770302108','0.000384222','0.677593822','0.024753668','0.607458162','0.018326723','0.607509806','0.605710062','0.606283124','0.0589118685876641','0.058987784','0.258055424','0.445397961','0.288493093','0.999233827'};
 
for i=1:length(fun)
 
    T=['XSteam(''',char(fun(i)),''',',char(In1(i)),',',char(In2(i)),')'];
    Res=eval(T);
    Error=(Res-(str2num(char(Control(i)))))/str2num(char(Control(i)));
    Check=[T,'=',num2str(Res),' - (Control)',char(Control(i)),'=',num2str(Error)]
    if Error>1E-5
        err=err+1;
        error('To large error')
    end   
end

MATLAB code that uses Xsteam to generate steam tables and graphs

This version of the code uses an Excel spreadsheet that must be placed in the same folder as the executable code. The spreadsheet is called parameters, and the first row contains the pressures. The columns contain the corresponding temperatures. Do not include the saturation temperature -- that is automatically included. The and the MATLAB code is called TableGraphMaker.m.^[2]

Copy of parameters

click to view or hide

25	50	100	150	200	250	300	350	400	450	500
200	200	200	200	200	200	200	200	200	200	200
250	250	250	250	250	250	250	250	250	300	300
300	300	300	300	300	300	300	300	300	350	350
350	350	350	350	350	350	350	350	350	400	400
400	400	400	375	375	370	375	375	380	425	425
450	450	450	400	390	380	400	400	400	450	450
500	500	500	425	400	390	415	415	425	500	500
			450	425	400	425	425	450
			500	450	425	450	450	500
				500	450	500	500
					500

Copy of TableGraphMaker.m

click to view or hide

clear all;close all;clc;
 
Tmin =200; Tmax=500;
Tphase = linspace(Tmin,373.9,500);
 
linewidth=1.2;%linewidth
markersize = 5;
% This makes a T_versus_s plot
 
% There are two user selected paramters that often seemed necessary.  Make
% them as small as possible.  The problem is that calling XSteam at a
% saturation value sometimes yields odd results.
 
%These cell arrays are used to call the various state variables
Call ={'v_pT', 'u_pT','h_pT','s_pT' };
%Label ={'v (m^3/kg)', 'u (kJ/kg)','h (kJ/kg)','s(kJ/kg)','T (&deg;C)', 'p (bars)' };
Label ={'v', 'u','h','s','&deg;C','bars' };
plotcall ={'v_pT', 'u_pT','h_pT','s_pT' }; 
phasecall ={'vL_T','uL_T','hL_T','sL_T'}; %For liquid-gas phase
phasecallG ={'vV_T','uV_T','hV_T','sV_T'}; %For liquid-gas phase
tempfix = .09; %ideal is zero.  This is added to T so that T critical calculates parameters
presfix=.000001; %added so that pressure is above critical
 
sigfigs = '%10.2E'; %Use  '%10.3E' for  big table and '%10.3E' for compact table.
%Here are some paramters we don't seem to need:  %Rwater=461.526/100000;  %Tcritical = 374; %Pcritical = 220.6395 %bar
%% Read parameters and replace NAN by -1
parameters= xlsread('parameters.xlsx'); %Omit saturation temperatures or they will be repeated
[NTfirst,Np]=size(parameters);
for nT=1:NTfirst %Fills non numbers with -1
    for np=1:Np
        if isnan(  parameters(nT,np) )==1 % i.e., we are above saturation pressure  '
            parameters(nT,np)=-1;
        end
    end
end
% Define p as pressure array
p=parameters(1,:); %Defines p;
 
%%  fill the first row of parameters with saturation temperatures, when they exist
% Also, define T: T(1) is large to ensure evaluation by XSteam, so we must
% subtract tempfix everytime we print the saturation temperature.
for np = 1:Np
    if isnan(   XSteam('Tsat_p',p(np))  )==0 %  if Tsat exists
        parameters(1,np)=XSteam('Tsat_p',p(np))+tempfix;%added tempfix to ensure return%
    else parameters(1,np)=-1;
    end
end
 
%% Make wikitable header
fout=fopen('wikitable.txt', 'w+');% creates wikitext
for np=1:Np
    T=parameters(:,np); %need to define the temperatures after
    stringp=num2str(p(np),'%10.2E'); % is the pressure
    if T(1)>0
        stringT = num2str(T(1)-tempfix,'%10.1f'); % is the saturation temperature
        commentInHeader{np} = [' bars  &nbsp; &nbsp; (T<sub>sat</sub> = ',stringT,' &deg;C)===='];
        belowCrit=true;
    else %above critical pressure
        belowCrit=false;
        commentInHeader{np} = ' bars &nbsp; &nbsp; (above critical)====';
    end
    fprintf(fout,'%s\n',['====P = ',stringp,commentInHeader{np}]);
    %fprintf(fout,'%s\n','{| class="wikitable" style="text-align:center;  width:600px;" ');
    fprintf(fout,'%s\n','{| class="wikitable" style="text-align:center; " ');
    fprintf(fout,'%s\n','|-');
    % fprintf(fout,'%s\n',['|',' ','||', ' ', '||', ' ', '||' ,' ', '||', stringp ]);fprintf(fout,'%s\n','|-');
    string = ['|', Label{5}, '||' , Label{1}, '||', Label{2}, '||', Label{3}, '||', Label{4} ];
    fprintf(fout,'%s\n',string);
    %% print entries
    for nT = 1:NTfirst
        printline=0;
 
        if nT==1 && T(nT)>0, printline=1; end;
        if nT>1 && T(nT)>T(1), printline=1; end;
 
        if printline==1
            stringT = num2str(T(nT),'%10.0f');
            if nT==1 stringT='Sat'; end;
            stringv = num2str(XSteam('v_pT',p(np), T(nT)),sigfigs);
            stringu = num2str(XSteam('u_pT', p(np), T(nT)),sigfigs);
            stringh = num2str(XSteam('h_pT', p(np), T(nT)),sigfigs);
            strings = num2str(XSteam('s_pT', p(np), T(nT)),sigfigs);
            fprintf(fout,'%s\n','|-');
            fprintf(fout,'%s\n',['|',stringT,'||', stringv, '||', stringu, '||' ,stringh, '||', strings ]);
            if nT==1
                T(nT)=T(nT)+tempfix; %restore, just in case
            end
        end
    end
    fprintf(fout,'%s\n','|}');
end
fclose(fout);
%%
 
for count=1:4 %try this breakup
    what2plot = plotcall{count};
 
    for np=1:Np %create Tplot and Xplot
        clear Tplot Xplot; %and build up the temperatures for a contour
        T=parameters(:,np);
        %plotthis=0;
        firstplot=1; %default values
        for nT = 1:NTfirst   %begin search of a point worth plotting
            if T(nT)>=T(1),
                % plothis=1;
                if firstplot==1 %first element in both arrays
 
                    Xplot = XSteam(plotcall{count}, p(np), T(nT) );
                    Tplot =T(nT);
                    firstplot = 0; %on next iteration do the following else
                else
                    Xnext = XSteam(plotcall{count}, p(np), T(nT) );
                    Tnext =T(nT);
                    Xplot = [Xplot,Xnext];
                    Tplot = [Tplot,Tnext];
                end %if firstplot==1 xxx else xxx
            end %if plotthis
        end; %for nT
        %%Make one  countour of the subplot
        subplot(2,2,count);
        if np == 1
          plot(Xplot,Tplot,':','LineWidth',linewidth) 
        else
        plot(Xplot,Tplot,'-o','LineWidth',linewidth,'MarkerSize',markersize)   
        end
        hold on; %Go to the next countour
 
    end %for np % we are done with all the countours
 
    % Add Liquid and Gas phase transition
    for n = 1:size(Tphase,2)
        Xphase(n) = XSteam(phasecall{count},Tphase(n) );
    end
    for n = 1:size(Tphase,2)
        XphaseG(n) = XSteam(phasecallG{count},Tphase(n));
    end
    plot(Xphase,Tphase,'r', 'Linewidth',linewidth)
    plot(XphaseG,Tphase,'m','Linewidth',linewidth)
 
    xlabel(Label{count});
    ylabel(['T=Temperature ( ^{\circ}C)']);
    if count==1 xlim([0 .1]); else xlim auto; end;
    ylim([Tmin Tmax+1]);
    hold off;
 
end %for count

Footnotes and references

↑ XSteam available at http://www.mathworks.com/matlabcentral/fileexchange/9817-x-steam--thermodynamic-properties-of-water-and-steam
↑ aka "student"

This article is issued from Wikiversity - version of the Tuesday, March 24, 2015. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.

25	50	100	150	200	250	300	350	400	450	500
200	200	200	200	200	200	200	200	200	200	200
250	250	250	250	250	250	250	250	250	300	300
300	300	300	300	300	300	300	300	300	350	350
350	350	350	350	350	350	350	350	350	400	400
400	400	400	375	375	370	375	375	380	425	425
450	450	450	400	390	380	400	400	400	450	450
500	500	500	425	400	390	415	415	425	500	500
			450	425	400	425	425	450
			500	450	425	450	450	500
				500	450	500	500
					500

25	50	100	150	200	250	300	350	400	450	500
200	200	200	200	200	200	200	200	200	200	200
250	250	250	250	250	250	250	250	250	300	300
300	300	300	300	300	300	300	300	300	350	350
350	350	350	350	350	350	350	350	350	400	400
400	400	400	375	375	370	375	375	380	425	425
450	450	450	400	390	380	400	400	400	450	450
500	500	500	425	400	390	415	415	425	500	500
			450	425	400	425	425	450
			500	450	425	450	450	500
				500	450	500	500
					500

25	50	100	150	200	250	300	350	400	450	500
200	200	200	200	200	200	200	200	200	200	200
250	250	250	250	250	250	250	250	250	300	300
300	300	300	300	300	300	300	300	300	350	350
350	350	350	350	350	350	350	350	350	400	400
400	400	400	375	375	370	375	375	380	425	425
450	450	450	400	390	380	400	400	400	450	450
500	500	500	425	400	390	415	415	425	500	500
			450	425	400	425	425	450
			500	450	425	450	450	500
				500	450	500	500
					500