Understand the concept of Rankine cycle by referring to standard thermal engineering text book.

• Create a Rankine Cycle Simulator using OCTAVE/MATLAB.
• Your code should calculate the state points of the Rankine Cycle (any type of your choice) based on user inputs. 
• Then, plot the corresponding T-s and h-s plots for the given set of inputs. Title and axes labels are a must, while legends could be used if necessary.
• The presentation of your code will be taken into consideration while grading. Labelling sections of your code using comments is appreciated.
• Attach a report explaining the working of a rankine cycle simulator and also the research done by you.

Required Inputs:

Turbine Inlet temperature = 400° C

Turbine Inlet Pressure = 30 bars

Turbine Outlet Pressure = 0.05 bars = Condenser Inlet Pressure

Calculate the followings,

• Net work output
• Back work ratio

Rankine cycle:

     The Rankine cycle is an idealized thermodynamic cycle describing the process by which certain heat engines, such as steam turbines or reciprocating steam engines, allow mechanical work to be extracted from a fluid as it moves between a heat source and heat sink. The Rankine cycle is named after William John Macquorn Rankine, a Scottish polymath professor at Glasgow University.

Heat energy is supplied to the system via a boiler where the working fluid (typically water) is converted to a high pressure gaseous state (steam) in order to turn a turbine. After passing over the turbine the fluid is allowed to condense back into a liquid state as waste heat energy is rejected before being returned to boiler, completing the cycle. Friction losses throughout the system are often neglected for the purpose of simplifying calculations as such losses are usually much less significant than thermodynamic losses, especially in larger systems.

There are four processes in the Rankine cycle. The states are identified by numbers (in brown) in the T–s diagram.

• Process 1–2: The working fluid is pumped from low to high pressure. As the fluid is a liquid at this stage, the pump requires little input energy. Process 1-2 is isentropic compression.
• Process 2–3: The high-pressure liquid enters a boiler, where it is heated at constant pressure by an external heat source to become a dry saturated vapour. The input energy required can be easily calculated graphically, using an enthalpy–entropy chart (h–s chart, or Mollier diagram), or numerically, using steam tables or software.  Process 2-3 is constant pressure heat addition in boiler.
• Process 3–4: The dry saturated vapour expands through a turbine, generating power. This decreases the temperature and pressure of the vapour, and some condensation may occur. The output in this process can be easily calculated using the chart or tables noted above. Process 3-4 is isentropic expansion.
• Process 4–1: The wet vapour then enters a condenser, where it is condensed at a constant pressure to become a saturated liquid. Process 4-1 is constant pressure heat rejection in condenser.
  In an ideal Rankine cycle the pump and turbine would be isentropic, i.e., the pump and turbine would generate no entropy and hence maximize the net work output. Processes 1–2 and 3–4 would be represented by vertical lines on the T–s diagram and more closely resemble that of the Carnot cycle. The Rankine cycle shown here prevents the state of the working fluid from ending up in the superheated vapor region after the expansion in the turbine,  which reduces the energy removed by the condensers.

The actual vapor power cycle differs from the ideal Rankine cycle because of irreversibilities in the inherent components caused by fluid friction and heat loss to the surroundings; fluid friction causes pressure drops in the boiler, the condenser, and the piping between the components, and as a result the steam leaves the boiler at a lower pressure; heat loss reduces the net work output, thus heat addition to the steam in the boiler is required to maintain the same level of net work output.

Xsteam Program(given):

%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

Code explanation:

• The basic data is printed about the Rankine cycle.
• Then the inputs provided are taken is as inputs.
• With the XSteam libarary file other inputs are taken from that library.
• The syntax for inputing the file data is alos provided.
• With the basis formula provided and formula taken from wikipedia the cycle value are calculated.
• Later the T-s and h-s plots for the given set of inputs are ploted.
• Then the required output calculated using the formulas.

Program to run Rankine Cycle:

clear all
close all
clc

fprintf('Ideal Rankine Cycle Simulator')
fprintf('\n\nProcess 1-2 is isentropic compression in the pump')
fprintf('\nProcess 2-3 is constant pressure heat addition in the boiler ')
fprintf('\nProcess 3-4 is isentropic expansion in the steam turbine')
fprintf('\nProcess 4-1 is constant pressure heat rejection in conderser')


fprintf('\nInput required')
p3 = input('\n1.  Enter the pressre at turbine inlet(in bar) = ')
T3 = input('\n2.  Enter the temperature at turbine inlet(in Celcius) = ')
p4 = input('\n3.  Enter the Pressure at Turbine Outlet / Condenser inlet(in bar) = ')

p1=p4;
T1 = XSteam('Tsat_p', p1);
h1 = XSteam('hL_p',p1);
s1 = XSteam('sL_p',p1);
v1 = XSteam('vL_p',p1);
sf1 = XSteam('sL_p', p1);

p2=p3;
s2=s1;
T2 = XSteam('T_ps',p2,s2);
v2 = XSteam('v_pT', p2,s2);
h2 = XSteam('h_ps',p2,s2);

Tsat_2 =XSteam('Tsat_p',p2);
sf2 = XSteam('sL_T',Tsat_2);
hf2 = XSteam('hL_T',Tsat_2);
sg2 = XSteam('sV_T', Tsat_2);
hg2 = XSteam('hV_T',Tsat_2);
Cp = XSteam('Cp_ps',p2,s2);
Cv = XSteam('Cv_ps',p2,s2);
C = Cp/Cv;

w_p = h2-h1;

h3 = XSteam('h_pt',p3,T3);
s3 =XSteam('s_pT',p3,T3);
v3 = XSteam('v_pT',p3,T3);
x3 = XSteam('vx_ph',p3,h3);

s4=s3;
x4 = XSteam('x_ps',p4,s4);
h4 = XSteam('h_ps',p4,s4);
T4 = XSteam('T_ph',p4,h4);
sg4 = XSteam('sV_p',p4);

q_in = h3-h2;

q_out = h4-h1;

w_t = h3-h4;

w_net = w_t-w_p;

Eff_th = (w_net/q_in)*100;

BWR = w_p/w_t;

t_s = linspace(5,400,500);
for i= 1:length(t_s)
    s(i) = XSteam('sL_T',t_s(i));
    S(i) = XSteam('sV_T', t_s(i));
    h(i) = XSteam('hL_T',t_s(i));
    H(i) = XSteam('hV_T',t_s(i));
end

figure(1)
hold on
line([s1 s2],[T1 T2],'color','r','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([s2 sf2],[T2 Tsat_2],'color','b','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([s3 s4],[T3 T4],'color','g','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([s4 s1],[T4 T1],'color','m','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([sf2 sg2],[Tsat_2 Tsat_2],'color','b','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([sg2 s3],[Tsat_2 T3],'color','g','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
plot(s,t_s,'color','black','linewidth',2);
plot(S,t_s,'color','black','linewidth',2);
xlabel('Entropy [KJ/K]');
ylabel('Temperature (deg-celcius)')
title('T-s diagram of ideal Rankine Cycle')
legend('Isentropic Compression in pump','Isobaric heating in boiler','Isentropic expansion in turbine','Isobaric heat rejection in condenser')

figure(2)
hold on
line([s1 s2],[h1 h2],'color','r','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([s2 sf2],[h2 hf2],'color','b','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([s3 s4],[h3 h4],'color','g','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([s4 s1],[h4 h1],'color','m','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([sf2 sg2],[hf2 hg2],'color','b','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
line([sg2 s3],[hg2 h3],'color','g','Marker','x','linewidth',2,'Markersize',10,'MarkerEdgeColor','k')
plot(s,h,'color','black','linewidth',2);
plot(S,H,'color','black','linewidth',2);
xlabel('Entropy [KJ/K]');
ylabel('Enthalpy [KJ/Kg')
title('h-s diagram of ideal Rankine Cycle')
legend('Isentropic Compression in pump','Isobaric heating in biler','Isentropic expansion in turbine','Isobaric heat rejection in condenser')

fprintf('\nOutput:\n')
fprintf('\nA.  State1(pump inlet) properties:')
fprintf('\n     1.Pressure=%d bar',p1)
fprintf('\n     2.Temperature=%d K',T1)
fprintf('\n     3.Enthalpy=%d KJ/Kg',h1)
fprintf('\n     4.Entropy=%d KJ/k',s1)

fprintf('\n\nB.  State2(Boiler inlet) properties:')
fprintf('\n     1.Pressure=%d bar',p2)
fprintf('\n     2.Temperature=%d K',T2)
fprintf('\n     3.Enthalpy=%d KJ/Kg',h2)
fprintf('\n     4.Entropy=%d KJ/k',s2)

fprintf('\n\nC.  Properties at Boiler pressure Saturation:')
fprintf('\n     1.Staturatoin Temperature=%d K',Tsat_2)
fprintf('\n     2.  Saturated Liquid enthalpy=%d KJ/Kg',hf2)
fprintf('\n     3.  Saturated Liquid entropy=%d KJ/K',sf2)
fprintf('\n     4.  Saturated vapour enthalpy=%d KJ/kg',hg2)
fprintf('\n     5.  Saturated vapour entropy=%d KJ/k',sg2)

fprintf('\n\nD.  State3(Turbine inlet) properties:')
fprintf('\n     1.  Pressure=%d bar',p3)
fprintf('\n     2.  Temperature=%d K',T3)
fprintf('\n     3.  Enthalpy=%d KJ/Kg',h3)
fprintf('\n     4.  Entropy=%d KJ/k',s3)
fprintf('\n     5.  Vapour volume Fraction = %d', x3)

fprintf('\n\nE.  State4(Condenser inlet) properties:')
fprintf('\n     1.  Pressure=%d bar',p4)
fprintf('\n     2.  Temperature=%d K',T4)
fprintf('\n     3.  Enthalpy=%d KJ/Kg',h4)
fprintf('\n     4.  Entropy=%d KJ/k',s4)
fprintf('\n     5.  Vapour volume Fraction = %d', x4)

fprintf('\n\nF.  Cycle Results:')
fprintf('\n     1.  Pump Work =%d KJ/Kg',w_p)
fprintf('\n     2.  Hest Boiler =%d KJ/Kg',q_in)
fprintf('\n     3.  Turbine work =%d KJ/Kg',w_t)
fprintf('\n     4.  Net work output =%d KJ/kg',w_net)
fprintf('\n     5.  Thermal Efficiency = %d ', Eff_th)
fprintf('\n     6.  Back Work ratio= %d ', BWR)

Program output:

Graph:

T-s plot:

h-s plot:

Error Faced during programing:

      There were couple of error made while programming but they were typical typo errors.  Hence, they were easily overcomed.

Screenshots showing errors thrown:

1.

2.

Inference:

     Thus a Rankine Cycle Simulator is created using Matlab and attached with deatiled report.