MATLAB/Engineering thermodynamics

This page is devoted to XSteam.m, which is a MATLAB code available available on the internet.



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 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:

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

Returns the saturation temperature of water at 1 bar.

Backup copy of XSteam.m
%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.

Copy of TableGraphMaker.m
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     (Tsat = ',stringT,' &deg;C)====']; belowCrit=true; else %above critical pressure belowCrit=false; commentInHeader{np} = ' bars    (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