//***************************************************************************************************************************
//* 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   
//================================================================================
//* Modified for SciLab by Bondarenko Oleksiy bond@nmri.go.jp                          
//****************************************************************************************************************************
//
// 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.
//                XSteam('test') Verifies calling funtions
//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.

clear
clearglobal
funcprot(0);

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=convstr(fun,'l');

select fun
//***********************************************************************************************************
    //*1.2 Temperature
case 'tsat_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657  & p < 22.06395) then 
        Out = fromSIunit_T(T4_p(p));
    else
        Out = %nan;
    end
case 'tsat_s' then
    s = toSIunit_s(In1);
    if (s > -0.0001545495919  & s < 9.155759395) then 
        ps = p4_s(s);
        Out = fromSIunit_T(T4_p(ps));
    else
        Out = %nan;
    end
case 't_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        Out = fromSIunit_T(T1_ph(p, h));
    case 2 then
        Out = fromSIunit_T(T2_ph(p, h));
    case 3 then
        Out = fromSIunit_T(T3_ph(p, h));
    case 4 then
        Out = fromSIunit_T(T4_p(p));
    case 5 then
        Out = fromSIunit_T(T5_ph(p, h));
    else
        Out = %nan;
    end
case 't_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Out = fromSIunit_T(T1_ps(p, s));
    case 2 then
        Out = fromSIunit_T(T2_ps(p, s));
    case 3 then
        Out = fromSIunit_T(T3_ps(p, s));
    case 4 then
        Out = fromSIunit_T(T4_p(p));
    case 5 then
        Out = fromSIunit_T(T5_ps(p, s));
    else
        Out = %nan;
    end 
    
case 't_hs' then
    h = toSIunit_h(In1);
    s = toSIunit_s(In2);
    Region = region_hs(h, s);
    select Region
    case 1 then
        p1 = p1_hs(h, s);
        Out = fromSIunit_T(T1_ph(p1, h));
    case 2 then
        p2 = p2_hs(h, s);
        Out = fromSIunit_T(T2_ph(p2, h));
    case 3 then
        p3 = p3_hs(h, s);
        Out = fromSIunit_T(T3_ph(p3, h));
    case 4 then
        Out = fromSIunit_T(T4_hs(h, s));
    case 5 then
        error('functions of hs is not avlaible in region 5');
    else
        Out = %nan;
    end
    //***********************************************************************************************************
    //*1.3 Pressure (p)
case 'psat_t' then
    T = toSIunit_T(In1);
    if (T < 647.096  & T > 273.15) then 
        Out = fromSIunit_p(p4_T(T));
    else
        Out = %nan;
    end
    
case 'psat_s' then
    s = toSIunit_s(In1);
    if (s > -0.0001545495919 & s < 9.155759395) then 
        Out = fromSIunit_p(p4_s(s));
    else
        Out = %nan;
    end
    
case 'p_hs' then
    h = toSIunit_h(In1);
    s = toSIunit_s(In2);
    Region = region_hs(h, s);
    select Region
    case 1 then
        Out = fromSIunit_p(p1_hs(h, s));
    case 2 then
        Out = fromSIunit_p(p2_hs(h, s));
    case 3 then
        Out = fromSIunit_p(p3_hs(h, s));
    case 4 then
        tSat = T4_hs(h, s);
        Out = fromSIunit_p(p4_T(tSat));
    case 5 then
        error('functions of hs is not avlaible in region 5');
    else
        Out = %nan;
    end
    
    case 'p_hrho' then
    h=In1;
    rho=In2;
//Not valid for water or supercritical 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 then
      High_Bound = ps;
    else
      Low_Bound = ps;
    end
    ps = (Low_Bound + High_Bound) / 2;
  end
    Out = ps;

    //***********************************************************************************************************
    //*1.4 Enthalpy (h)
case 'hv_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        Out = fromSIunit_h(h4V_p(p));
    else
        Out = %nan;
    end
    
case 'hl_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        Out = fromSIunit_h(h4L_p(p));
    else
        Out = %nan;
    end
    
case 'hv_t' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        p = p4_T(T);
        Out = fromSIunit_h(h4V_p(p));
    else
        Out = %nan;
    end
    
case 'hl_t' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        p = p4_T(T);
        Out = fromSIunit_h(h4L_p(p));
    else
        Out = %nan;
    end
    
    
case 'h_pt' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_h(h1_pT(p, T));
    case 2 then
        Out = fromSIunit_h(h2_pT(p, T));
    case 3 then
        Out = fromSIunit_h(h3_pT(p, T));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_h(h5_pT(p, T));
    else
        Out = %nan;
    end
    
case 'h_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Out = fromSIunit_h(h1_pT(p, T1_ps(p, s)));
    case 2 then
        Out = fromSIunit_h(h2_pT(p, T2_ps(p, s)));
    case 3 then
        Out = fromSIunit_h(h3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)));
    case 4 then
        xs = x4_ps(p, s);
        Out = fromSIunit_h(xs * h4V_p(p) + (1 - xs) * h4L_p(p));
    case 5 then
        Out = fromSIunit_h(h5_pT(p, T5_ps(p, s)));
    else
        Out = %nan;
    end
    
case 'h_px' then
    p = toSIunit_p(In1);
    x = toSIunit_x(In2);
    if (x > 1 | x < 0 | p >= 22.064) then 
        Out = %nan;
        return
    end
    hL = h4L_p(p);
    hV = h4V_p(p);
    Out = hL + x * (hV - hL);
    
case 'h_prho' then
    p = toSIunit_p(In1);
    rho = 1 / toSIunit_v(1 / In2);
    Region = Region_prho(p, rho);
    select Region
        case 1 then
            Out = fromSIunit_h(h1_pT(p, T1_prho(p, rho)));
        case 2 then
            Out = fromSIunit_h(h2_pT(p, T2_prho(p, rho)));
        case 3 then
            Out = fromSIunit_h(h3_rhoT(rho, T3_prho(p, rho)));
        case 4 then
            if p < 16.529 then
                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 then
            Out = fromSIunit_h(h5_pT(p, T5_prho(p, rho)));
        else
            Out = %nan;
    end

case 'h_tx' then
    T = toSIunit_T(In1);
    x = toSIunit_x(In2);
    if (x > 1 | x < 0 | T >= 647.096) then
        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' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            Out = fromSIunit_v(v2_pT(p, T4_p(p)));
        else
            Out = fromSIunit_v(v3_ph(p, h4V_p(p)));
        end
    else
        Out = %nan;
    end
     
case 'vl_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then 
            Out = fromSIunit_v(v1_pT(p, T4_p(p)));
        else
            Out = fromSIunit_v(v3_ph(p, h4L_p(p)));
        end
    else
        Out = %nan;
    end
        
case 'vv_t' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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
    
case 'vl_t' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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
      
case 'v_pt' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_v(v1_pT(p, T));
    case 2 then
        Out = fromSIunit_v(v2_pT(p, T));
    case 3 then
        Out = fromSIunit_v(v3_ph(p, h3_pT(p, T)));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_v(v5_pT(p, T));
    else
        Out = %nan;
    end
      
case 'v_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region 
    case 1 then
        Out = fromSIunit_v(v1_pT(p, T1_ph(p, h)));
    case 2 then
        Out = fromSIunit_v(v2_pT(p, T2_ph(p, h)));
    case 3 then
        Out = fromSIunit_v(v3_ph(p, h));
    case 4 then
        xs = x4_ph(p, h);
        if p < 16.529 then
            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 then
        Ts = T5_ph(p, h);
        Out = fromSIunit_v(v5_pT(p, Ts));
    else
        Out = %nan;
    end
 
case 'v_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Out = fromSIunit_v(v1_pT(p, T1_ps(p, s)));
    case 2 then
        Out = fromSIunit_v(v2_pT(p, T2_ps(p, s)));
    case 3 then
        Out = fromSIunit_v(v3_ps(p, s));
    case 4 then
        xs = x4_ps(p, s);
        if p < 16.529 then
            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 then
        Ts = T5_ps(p, s);
        Out = fromSIunit_v(v5_pT(p, Ts));
    else
        Out = %nan;
    end
     
    //***********************************************************************************************************
    //*1.6 Density (rho)
    // Density is calculated as 1/v. 
case 'rhov_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            Out = 1/fromSIunit_v(v2_pT(p, T4_p(p)));
        else
            Out = 1/fromSIunit_v(v3_ph(p, h4V_p(p)));
        end
    else
        Out = %nan;
    end

case 'rhol_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then 
            Out = 1/fromSIunit_v(v1_pT(p, T4_p(p)));
        else
            Out = 1/fromSIunit_v(v3_ph(p, h4L_p(p)));
        end
    else
        Out = %nan;
    end
        
case 'rhov_t' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            Out = 1/fromSIunit_v(v2_pT(p4_T(T), T));
        else
            Out = 1/fromSIunit_v(v3_ph(p4_T(T), h4V_p(p4_T(T))));
        end
    else
        Out = %nan;
    end
    
case 'rhol_t' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            Out = 1/fromSIunit_v(v1_pT(p4_T(T), T));
        else
            Out = 1/fromSIunit_v(v3_ph(p4_T(T), h4L_p(p4_T(T))));
        end
    else
        Out = %nan;
    end
      
case 'rho_pt' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = 1/fromSIunit_v(v1_pT(p, T));
    case 2 then
        Out = 1/fromSIunit_v(v2_pT(p, T));
    case 3 then
        Out = 1/fromSIunit_v(v3_ph(p, h3_pT(p, T)));
    case 4 then
        Out = %nan;
    case 5 then
        Out = 1/fromSIunit_v(v5_pT(p, T));
    else
        Out = %nan;
    end
      
case 'rho_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region 
    case 1 then
        Out = 1/fromSIunit_v(v1_pT(p, T1_ph(p, h)));
    case 2 then
        Out = 1/fromSIunit_v(v2_pT(p, T2_ph(p, h)));
    case 3 then
        Out = 1/fromSIunit_v(v3_ph(p, h));
    case 4 then
        xs = x4_ph(p, h);
        if p < 16.529 then
            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 = 1/fromSIunit_v((xs * v4v + (1 - xs) * v4L));
    case 5 then
        Ts = T5_ph(p, h);
        Out = 1/fromSIunit_v(v5_pT(p, Ts));
    else
        Out = %nan;
    end
 
case 'rho_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Out = 1/fromSIunit_v(v1_pT(p, T1_ps(p, s)));
    case 2 then
        Out = 1/fromSIunit_v(v2_pT(p, T2_ps(p, s)));
    case 3 then
        Out = 1/fromSIunit_v(v3_ps(p, s));
    case 4 then
        xs = x4_ps(p, s);
        if p < 16.529 then
            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 = 1/fromSIunit_v((xs * v4v + (1 - xs) * v4L));
    case 5 then
        Ts = T5_ps(p, s);
        Out = 1/fromSIunit_v(v5_pT(p, Ts));
    else
        Out = %nan;
    end	
	
    //***********************************************************************************************************
    //*1.7 Specific entropy (s)
case 'sv_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_s(s1_pT(p, T));
    case 2 then
        Out = fromSIunit_s(s2_pT(p, T));
    case 3 then
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_s(s3_rhoT(rhos, T));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_s(s5_pT(p, T));
    else
        Out = %nan;
    end
    
case 's_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        T = T1_ph(p, h);
        Out = fromSIunit_s(s1_pT(p, T));
    case 2 then
        T = T2_ph(p, h);
        Out = fromSIunit_s(s2_pT(p, T));
    case 3 then
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_s(s3_rhoT(rhos, Ts));
    case 4 then
        Ts = T4_p(p);
        xs = x4_ph(p, h);
        if p < 16.529 then 
            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 then
        T = T5_ph(p, h);
        Out = fromSIunit_s(s5_pT(p, T));
    else
        Out = %nan;
    end
    
    
    //***********************************************************************************************************
    //*1.8 Specific internal energy (u)
case 'uv_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then
        if p < 16.529 then
            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' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then
        if p < 16.529 then 
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then
        if T <= 623.15 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_u(u1_pT(p, T));
    case 2 then
        Out = fromSIunit_u(u2_pT(p, T));
    case 3 then
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_u(u3_rhoT(rhos, T));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_u(u5_pT(p, T));
    else
        Out = %nan;
    end
    
case 'u_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        Ts = T1_ph(p, h);
        Out = fromSIunit_u(u1_pT(p, Ts));
    case 2 then
        Ts = T2_ph(p, h);
        Out = fromSIunit_u(u2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_u(u3_rhoT(rhos, Ts));
    case 4 then
        Ts = T4_p(p);
        xs = x4_ph(p, h);
        if p < 16.529 then 
            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 then
        Ts = T5_ph(p, h);
        Out = fromSIunit_u(u5_pT(p, Ts));
    else
        Out = %nan;
    end
    
case 'u_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Ts = T1_ps(p, s);
        Out = fromSIunit_u(u1_pT(p, Ts));
    case 2 then
        Ts = T2_ps(p, s);
        Out = fromSIunit_u(u2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_u(u3_rhoT(rhos, Ts));
    case 4 then
        if p < 16.529 then
            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 then
        Ts = T5_ps(p, s);
        Out = fromSIunit_u(u5_pT(p, Ts));
    else
        Out = %nan;
    end
    //***********************************************************************************************************
    //*1.9 Specific isobaric heat capacity (Cp)
case 'cpv_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_Cp(Cp1_pT(p, T));
    case 2 then
        Out = fromSIunit_Cp(Cp2_pT(p, T));
    case 3 then
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_Cp(Cp3_rhoT(rhos, T));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_Cp(Cp5_pT(p, T));
    else
        Out = %nan;
    end
    
case 'cp_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        Ts = T1_ph(p, h);
        Out = fromSIunit_Cp(Cp1_pT(p, Ts));
    case 2 then
        Ts = T2_ph(p, h);
        Out = fromSIunit_Cp(Cp2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_Cp(Cp3_rhoT(rhos, Ts));
    case 4 then
        Out = %nan;
    case 5 then
        Ts = T5_ph(p, h);
        Out = fromSIunit_Cp(Cp5_pT(p, Ts));
    else
        Out = %nan;
    end
    
case 'cp_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Ts = T1_ps(p, s);
        Out = fromSIunit_Cp(Cp1_pT(p, Ts));
    case 2 then
        Ts = T2_ps(p, s);
        Out = fromSIunit_Cp(Cp2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_Cp(Cp3_rhoT(rhos, Ts));
    case 4 then
        Out = %nan;
    case 5 then
        Ts = T5_ps(p, s);
        Out = fromSIunit_Cp(Cp5_pT(p, Ts));
    else
        Out = %nan;
    end
    
    //***********************************************************************************************************
    //*1.10 Specific isochoric heat capacity (Cv)
case 'cvv_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then
        if T <= 623.15 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_Cv(Cv1_pT(p, T));
    case 2 then
        Out = fromSIunit_Cv(Cv2_pT(p, T));
    case 3 then
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_Cv(Cv3_rhoT(rhos, T));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_Cv(Cv5_pT(p, T));
    else
        Out = %nan;
    end
    
case 'cv_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        Ts = T1_ph(p, h);
        Out = fromSIunit_Cv(Cv1_pT(p, Ts));
    case 2 then
        Ts = T2_ph(p, h);
        Out = fromSIunit_Cv(Cv2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_Cv(Cv3_rhoT(rhos, Ts));
    case 4 then
        Out = %nan;
    case 5 then
        Ts = T5_ph(p, h);
        Out = fromSIunit_Cv(Cv5_pT(p, Ts));
    else
        Out = %nan;
    end
    
    
case 'cv_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Ts = T1_ps(p, s);
        Out = fromSIunit_Cv(Cv1_pT(p, Ts));
    case 2 then
        Ts = T2_ps(p, s);
        Out = fromSIunit_Cv(Cv2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_Cv(Cv3_rhoT(rhos, Ts));
    case 4 then
        Out = %nan;  //(xs * CvVp + (1 - xs) * CvLp) / Cv_scale - Cv_offset
    case 5 then
        Ts = T5_ps(p, s);
        Out = fromSIunit_Cv(Cv5_pT(p, Ts));
    else
        Out = CVErr(xlErrValue);
    end
    
    
    //***********************************************************************************************************
    //*1.11 Speed of sound
case 'wv_p' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    p = toSIunit_p(In1);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    T = toSIunit_T(In1);
    if (T > 273.15 & T < 647.096) then 
        if T <= 623.15 then
            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' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
    case 1 then
        Out = fromSIunit_w(w1_pT(p, T));
    case 2 then
        Out = fromSIunit_w(w2_pT(p, T));
    case 3 then
        hs = h3_pT(p, T);
        rhos = 1 / v3_ph(p, hs);
        Out = fromSIunit_w(w3_rhoT(rhos, T));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_w(w5_pT(p, T));
    else
        Out = %nan;
    end
    
    
case 'w_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        Ts = T1_ph(p, h);
        Out = fromSIunit_w(w1_pT(p, Ts));
    case 2 then
        Ts = T2_ph(p, h);
        Out = fromSIunit_w(w2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ph(p, h);
        Ts = T3_ph(p, h);
        Out = fromSIunit_w(w3_rhoT(rhos, Ts));
    case 4 then
        Out = %nan;
    case 5 then
        Ts = T5_ph(p, h);
        Out = fromSIunit_w(w5_pT(p, Ts));
    else
        Out = %nan;
    end
    
case 'w_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    Region = region_ps(p, s);
    select Region
    case 1 then
        Ts = T1_ps(p, s);
        Out = fromSIunit_w(w1_pT(p, Ts));
    case 2 then
        Ts = T2_ps(p, s);
        Out = fromSIunit_w(w2_pT(p, Ts));
    case 3 then
        rhos = 1 / v3_ps(p, s);
        Ts = T3_ps(p, s);
        Out = fromSIunit_w(w3_rhoT(rhos, Ts));
    case 4 then
        Out = %nan; %(xs * wVp + (1 - xs) * wLp) / w_scale - w_offset
    case 5 then
        Ts = T5_ps(p, s);
        Out = fromSIunit_w(w5_pT(p, Ts));
    else
        Out = %nan;
    end
    //***********************************************************************************************************
    //*1.12 Viscosity
case 'my_pt' then
    p = toSIunit_p(In1);
    T = toSIunit_T(In2);
    Region = region_pT(p, T);
    select Region
	case 1 then
	    Out = fromSIunit_my(my_AllRegions_pT(p, T));
	case 2 then
	    Out = fromSIunit_my(my_AllRegions_pT(p, T));
	case 3 then
	    Out = fromSIunit_my(my_AllRegions_pT(p, T));	
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_my(my_AllRegions_pT(p, T));
    else
        Out = %nan;
    end
    
case 'my_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    Region = region_ph(p, h);
    select Region
    case 1 then
        Out = fromSIunit_my(my_AllRegions_ph(p, h));
    case 2 then
        Out = fromSIunit_my(my_AllRegions_ph(p, h));
    case 3 then
        Out = fromSIunit_my(my_AllRegions_ph(p, h));
    case 4 then
        Out = %nan;
    case 5 then
        Out = fromSIunit_my(my_AllRegions_ph(p, h));		
    else
        Out = %nan;
    end
    
case 'my_ps' then
    h = XSteam('h_ps',In1, In2);
    Out = XSteam('my_ph',In1, h);
    
    //***********************************************************************************************************
    //*1.13 Prandtl
case 'pr_pt' then
  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' then
  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' then
    T = toSIunit_T(In1);
    Out = fromSIunit_st(Surface_Tension_T(T));
    
case 'st_p' then
    T = XSteam('Tsat_p',In1);
    T = toSIunit_T(T);
    Out = fromSIunit_st(Surface_Tension_T(T));
    
    //***********************************************************************************************************
    //*1.16 Thermal conductivity
case 'tcl_p' then
    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' then
    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' then
    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' then
    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' then
    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' then
    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' then
    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' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    if (p > 0.000611657 & p < 22.06395) then 
        Out = fromSIunit_x(x4_ph(p, h));
    else
        Out = %nan;
    end
    
case 'x_ps' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    if (p > 0.000611657 & p < 22.06395) then 
        Out = fromSIunit_x(x4_ps(p, s));
    else
        Out = %nan;
    end
    
    //***********************************************************************************************************
    //*1.18 Vapour Volume Fraction
case 'vx_ph' then
    p = toSIunit_p(In1);
    h = toSIunit_h(In2);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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' then
    p = toSIunit_p(In1);
    s = toSIunit_s(In2);
    if (p > 0.000611657 & p < 22.06395) then 
        if p < 16.529 then
            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 'test' then
   Out=check();
    
else
    error('Unknown calling function to XSteam Fun, See help XSteam for valid calling functions');
end

endfunction

//***********************************************************************************************************
//*2 IAPWS IF 97 Calling functions                                                                          *
//***********************************************************************************************************
//
//***********************************************************************************************************
//*2.1 Functions for region 1

//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
function fun = v1_pT(p, T)
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

gamma_der_pi = - sum(n1 .* I1 .* (ones(n1)*7.1 - Pi) .^ (I1 - 1) .* (ones(n1)*tau - 1.222) .^ J1);

fun = R * T / p * Pi * gamma_der_pi / 1000;
endfunction

//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
function fun = h1_pT(p, T)
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
gamma_der_tau = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* J1 .* (ones(n1)*tau - 1.222) .^ (J1 - 1));

fun = R * T * tau * gamma_der_tau;
endfunction

//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
function fun = u1_pT(p, T)
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

gamma_der_tau = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* J1 .* (ones(n1)*tau - 1.222) .^ (J1 - 1));
gamma_der_pi = - sum(n1 .* I1 .* (ones(n1)*7.1 - Pi) .^ (I1 - 1) .* (ones(n1)*tau - 1.222) .^ J1);


fun = R * T * (tau * gamma_der_tau - Pi * gamma_der_pi);
endfunction

//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
function fun = s1_pT(p, T)
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;

//gam_ma = 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));
//    gam_ma = gam_ma + n1(i) * (7.1 - Pi) ^ I1(i) * (tau - 1.222) ^ J1(i);
//end
gam_ma = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* (ones(n1)*tau - 1.222) .^ J1);
gamma_der_tau = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* J1 .* (ones(n1)*tau - 1.222) .^ (J1 - 1));

fun = R * tau * gamma_der_tau - R * gam_ma;
endfunction

//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
function fun = Cp1_pT(p, T)
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

gamma_der_tautau = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* J1 .* (J1 - 1) .* (ones(n1)*tau - 1.222) .^ (J1 - 2));

fun = -R * tau ^ 2 * gamma_der_tautau;
endfunction

//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
function fun = Cv1_pT(p, T)
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

gamma_der_pi = -sum(n1 .* I1 .* (ones(n1)*7.1 - Pi) .^ (I1 - 1) .* (ones(n1)*tau - 1.222) .^ J1);
gamma_der_pipi = sum(n1 .* I1 .* (I1 - 1) .* (ones(n1)*7.1 - Pi) .^ (I1 - 2) .* (ones(n1)*tau - 1.222) .^ J1);
gamma_der_pitau = -sum(n1 .* I1 .* (ones(n1)*7.1 - Pi) .^ (I1 - 1) .* J1 .* (ones(n1)*tau - 1.222) .^ (J1 - 1));
gamma_der_tautau = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* J1 .* (J1 - 1) .* (ones(n1)*tau - 1.222) .^ (J1 - 2));

fun = R * (-tau ^ 2 * gamma_der_tautau + (gamma_der_pi - tau * gamma_der_pitau) ^ 2 / gamma_der_pipi);
endfunction

//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
function fun = w1_pT(p, T)
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

gamma_der_pi = -sum(n1 .* I1 .* (ones(n1)*7.1 - Pi) .^(I1 - 1) .*(ones(n1)*tau - 1.222) .^J1);
gamma_der_pipi = sum(n1 .* I1 .* (I1 - 1) .* (ones(n1)*7.1 - Pi) .^ (I1 - 2) .* (ones(n1)*tau - 1.222) .^ J1);
gamma_der_pitau = -sum(n1 .* I1 .* (ones(n1)*7.1 - Pi) .^ (I1 - 1) .* J1 .* (ones(n1)*tau - 1.222) .^ (J1 - 1));
gamma_der_tautau = sum(n1 .* (ones(n1)*7.1 - Pi) .^ I1 .* J1 .* (J1 - 1) .* (ones(n1)*tau - 1.222) .^ (J1 - 2));

fun = (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;
endfunction

//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
function fun = T1_ph(p, h)
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

fun = sum(n1 .*(Pi*ones(n1)).^I1 .*(eta*ones(n1)+1).^J1);
endfunction

//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
function fun = T1_ps(p, s)
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

fun = sum(n1 .* (ones(n1)*Pi) .^ I1 .* (ones(n1)*Sigma + 2) .^ J1);
endfunction

//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
function fun = p1_hs(h, s)
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

fun = sum(n1 .* (ones(n1)*eta + 0.05) .^ I1 .* (ones(n1)*Sigma + 0.05) .^ J1) * 100;
endfunction

  //Solve by iteration. Observe that for low temperatures this equation has 2 solutions.
  //Solve with half interval method
function fun = T1_prho(p ,rho)
  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 then
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
  fun = Ts;
endfunction

//***********************************************************************************************************
//*2.2 functions for region 2

//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
function fun = v2_pT(p, T)
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

fun = R * T / p * Pi * (g0_pi + sum( nr .* Ir .* (ones(nr)*Pi) .^ (Ir - 1) .* (ones(nr)*tau - 0.5) .^ Jr)) / 1000;
endfunction

//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
function fun = h2_pT(p, T)
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

fun = R * T * tau * (sum(n0 .* J0 .* (ones(n0)*tau) .^ (J0 - 1)) + sum(nr .* (ones(nr)*Pi) .^ Ir .* Jr .* (ones(nr)*tau - 0.5) .^ (Jr - 1)));
endfunction

//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
function fun = u2_pT(p, T)
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

g0_tau = sum(n0 .* J0 .* (ones(n0)*tau) .^ (J0 - 1));
gr_pi = sum(nr .* Ir .* (ones(nr)*Pi) .^ (Ir - 1) .* (ones(nr)*tau - 0.5) .^ Jr);
gr_tau = sum(nr .* (ones(nr)*Pi) .^ Ir .* Jr .* (ones(nr)*tau - 0.5) .^ (Jr - 1));

fun = R * T * (tau * (g0_tau + gr_tau) - Pi * (g0_pi + gr_pi));
endfunction

//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
function fun = s2_pT(p, T)
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 = 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
g0 = log(Pi) + sum(n0 .* (ones(n0)*tau) .^ J0);
g0_tau = sum(n0 .* J0 .* (ones(n0)*tau) .^ (J0 - 1));
gr = sum( nr .* (ones(nr)*Pi) .^ Ir .* (ones(nr)*tau - 0.5) .^ Jr);
gr_tau = sum(nr .* (ones(nr)*Pi) .^ Ir .* Jr .* (ones(nr)*tau - 0.5) .^ (Jr - 1));

fun = R * (tau * (g0_tau + gr_tau) - (g0 + gr));
endfunction

//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
function fun = Cp2_pT(p, T)
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

g0_tautau = sum(n0 .* J0 .* (J0 - 1) .* (ones(n0)*tau) .^ (J0 - 2));
gr_tautau = sum(nr .* (ones(nr)*Pi) .^ Ir .* Jr .* (Jr - 1) .* (ones(nr)*tau - 0.5) .^ (Jr - 2));

fun = -R * tau ^ 2 * (g0_tautau + gr_tautau);
endfunction

//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
function fun = Cv2_pT(p, T)
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

g0_tautau = sum(n0 .* J0 .* (J0 - 1) .* (ones(n0)*tau) .^ (J0 - 2));
gr_pi = sum(nr .* Ir .* (ones(nr)*Pi) .^ (Ir - 1) .* (ones(nr)*tau - 0.5) .^ Jr);
gr_pipi = sum( nr .* Ir .* (Ir - 1) .* (ones(nr)*Pi) .^ (Ir - 2) .* (ones(nr)*tau - 0.5) .^ Jr);
gr_pitau = sum(nr .* Ir .* (ones(nr)*Pi) .^ (Ir - 1) .* Jr .* (ones(nr)*tau - 0.5) .^ (Jr - 1));
gr_tautau = sum(nr .* (ones(nr)*Pi) .^ Ir .* Jr .* (Jr - 1) .* (ones(nr)*tau - 0.5) .^ (Jr - 2));

fun = R * (-tau ^ 2 * (g0_tautau + gr_tautau) - (1 + Pi * gr_pi - tau * Pi * gr_pitau) ^ 2 / (1 - Pi ^ 2 * gr_pipi));
endfunction

//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
function fun = w2_pT(p, T)
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

g0_tautau = sum(n0 .* J0 .* (J0 - 1) .* (ones(n0)*tau) .^ (J0 - 2));
gr_pi = sum(nr .* Ir .* (ones(nr)*Pi) .^ (Ir - 1) .* (ones(nr)*tau - 0.5) .^ Jr);
gr_pipi = sum(nr .* Ir .* (Ir - 1) .* (ones(nr)*Pi) .^ (Ir - 2) .* (ones(nr)*tau - 0.5) .^ Jr);
gr_pitau = sum( nr .* Ir .* (ones(nr)*Pi) .^ (Ir - 1) .* Jr .* (ones(nr)*tau - 0.5) .^ (Jr - 1));
gr_tautau = sum(nr .* (ones(nr)*Pi) .^ Ir .* Jr .* (Jr - 1) .* (ones(nr)*tau - 0.5) .^ (Jr - 2));

fun = (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;
endfunction

//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
function fun = T2_ph(p, h)

if p < 4 then
    sub_reg = 1;
else
    if p < (905.84278514723 - 0.67955786399241 * h + 1.2809002730136E-04 * h ^ 2) then
        sub_reg = 2;
    else
        sub_reg = 3;
    end
end

select sub_reg
case 1 then
    //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
    fun = sum(ni .* (ones(ni)*p) .^ Ii .* (ones(ni)*hs - 2.1) .^ Ji);
case 2 then
    //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
    fun = sum(ni .* (ones(ni)*p - 2) .^ Ii .* (ones(ni)*hs - 2.6) .^ Ji);
else
    //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
    fun = sum(ni .* (ones(ni)*p + 25) .^ Ii .* (ones(ni)*hs - 1.8) .^ Ji);
end
endfunction

//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
function fun = T2_ps(p, s)

if p < 4 then
    sub_reg = 1;
else
    if s < 5.85 then
        sub_reg = 3;
    else
        sub_reg = 2;
    end
end
select sub_reg
case 1 then
    //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
    fun = sum(ni .* (ones(ni)*Pi) .^ Ii .* (ones(ni)*Sigma - 2) .^ Ji);
case 2 then
    //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
    fun = sum(ni .* (ones(ni)*Pi) .^ Ii .* (10 *ones(ni) - Sigma) .^ Ji);
else
    //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
    fun = sum(ni .* (ones(ni)*Pi) .^ Ii .* (2*ones(ni) - Sigma) .^ Ji);
end
endfunction

//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
function fun = p2_hs(h, s)

if h < (-3498.98083432139 + 2575.60716905876 * s - 421.073558227969 * s ^ 2 + 27.6349063799944 * s ^ 3) then
    sub_reg = 1;
else
    if s < 5.85 then
        sub_reg = 3;
    else
        sub_reg = 2;
    end
end
select sub_reg
case 1 then
    //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
    Pi = sum(ni .* (ones(ni)*eta - 0.5) .^ Ii .* (ones(ni)*Sigma - 1.2) .^ Ji);
    fun = Pi ^ 4 * 4;
case 2 then
    //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
     Pi = sum(ni .* (ones(ni)*eta - 0.6) .^ Ii .* (ones(ni)*Sigma - 1.01) .^ Ji);
    fun = Pi ^ 4 * 100;
else
    //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
    Pi = sum(ni .* (ones(ni)*eta - 0.7) .^ Ii .* (ones(ni)*Sigma - 1.1) .^ Ji);
    fun = Pi ^ 4 * 100;
end
endfunction

  //Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  //Solve with half interval method
function fun=T2_prho(p,rho) 

  if p < 16.5292 then
    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 then
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
    fun = Ts;
endfunction

//***********************************************************************************************************
//*2.3 functions for region 3

//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
function fun = p3_rhoT(rho, T)
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;
fidelta = ni(1)/delta + sum(ni(2:$) .* Ii(2:$) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 1) .* (ones(ni(2:$))*tau) .^ Ji(2:$));
fun = rho * R * T * delta * fidelta / 1000;
endfunction

//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
function fun = u3_rhoT(rho, T)
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
fitau = sum(ni(2:$) .* (ones(ni(2:$))*delta) .^ Ii(2:$) .* Ji(2:$) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 1));
fun = R * T * (tau * fitau);
endfunction

//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
function fun = h3_rhoT(rho, T)
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 = ni(1) / delta + sum(ni(2:$) .* Ii(2:$) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 1) .* (ones(ni(2:$))*tau) .^ Ji(2:40));
fitau = sum(ni(2:$) .* (ones(ni(2:$))*delta) .^ Ii(2:$) .* Ji(2:$) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 1));
fun = R * T * (tau * fitau + delta * fidelta);
endfunction


//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
function fun = s3_rhoT(rho, T)
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 = ni(1) * log(delta) + sum(ni(2:$) .* (ones(ni(2:$))*delta) .^ Ii(2:$) .* (ones(ni(2:$))*tau) .^ Ji(2:$));
fitau = sum(ni(2:$) .* (ones(ni(2:$))*delta) .^ Ii(2:$) .* Ji(2:$) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 1));
fun = R * (tau * fitau - fi);
endfunction

//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
function fun = Cp3_rhoT(rho, T)
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 = ni(1) / delta + sum(ni(2:$) .* Ii(2:$) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 1) .* (ones(ni(2:$))*tau) .^ Ji(2:$));
fideltadelta = sum( ni(2:$) .* Ii(2:$) .* (Ii(2:$) - 1) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 2) .* (ones(ni(2:$))*tau) .^ Ji(2:$)) - ni(1) / (delta ^ 2);
fitautau = sum( ni(2:$) .* (ones(ni(2:$))*delta) .^ Ii(2:$) .* Ji(2:$) .* (Ji(2:$) - 1) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 2));
fideltatau = sum(ni(2:$) .* Ii(2:$) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 1) .* Ji(2:$) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 1));

fun = R * (-tau ^ 2 * fitautau + (delta * fidelta - delta * tau * fideltatau) ^ 2 / (2 * delta * fidelta + delta ^ 2 * fideltadelta));
endfunction

//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
function fun = Cv3_rhoT(rho, T)
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

fitautau = sum(ni .* (ones(ni)*delta) .^ Ii .* Ji .* (Ji - 1) .* (ones(ni)*tau) .^ (Ji - 2));
fun = R * -(tau * tau * fitautau);
endfunction

//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
function fun = w3_rhoT(rho, T)
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 = ni(1) / delta + sum(ni(2:$) .* Ii(2:$) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 1) .* (ones(ni(2:$))*tau) .^ Ji(2:$));
fideltadelta = sum(ni(2:$) .* Ii(2:$) .* (Ii(2:$) - 1) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 2) .* (ones(ni(2:$))*tau) .^ Ji(2:$)) - ni(1) / (delta ^ 2);
fitautau = sum(ni(2:$) .* (ones(ni(2:$))*delta) .^ Ii(2:$) .* Ji(2:$) .* (Ji(2:$) - 1) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 2));
fideltatau = sum(ni(2:$) .* Ii(2:$) .* (ones(ni(2:$))*delta) .^ (Ii(2:$) - 1) .* Ji(2:$) .* (ones(ni(2:$))*tau) .^ (Ji(2:$) - 1));

fun = (1000 * R * T * (2 * delta * fidelta + delta ^ 2 * fideltadelta - (delta * fidelta - delta * tau * fideltatau) ^ 2 / (tau ^ 2 * fitautau))) ^ 0.5;
endfunction

//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
function fun = T3_ph(p, h)

h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3;

if h < h3ab then
    //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
    fun = sum(ni .* (ones(ni)*ps + 0.24) .^ Ii .* (ones(ni)*hs - 0.615) .^ Ji) * 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

    fun = sum(ni .* (ones(ni)*ps + 0.298) .^ Ii .* (ones(ni)*hs - 0.72) .^ Ji) * 860;
end
endfunction

//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
function fun = v3_ph(p, h)

h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3;

if h < h3ab then
    //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

    fun = sum(ni .* (ones(ni)*ps + 0.128) .^ Ii .* (ones(ni)*hs - 0.727) .^ Ji) * 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

    fun = sum(ni .* (ones(ni)*ps + 0.0661) .^ Ii .* (ones(ni)*hs - 0.72) .^ Ji) * 0.0088;
end
endfunction

//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
function fun = T3_ps(p, s)

if s <= 4.41202148223476 then
    //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

    fun = sum(ni .* (ones(ni)*Pi + 0.24) .^ Ii .* (ones(ni)*Sigma - 0.703) .^ Ji) * 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

    fun = sum(ni .* (ones(ni)*Pi + 0.76) .^ Ii .* (ones(ni)*Sigma - 0.818) .^ Ji) * 860;
end
endfunction

//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
function fun = v3_ps(p, s)

if s <= 4.41202148223476 then
    //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

    fun = sum(ni .* (ones(ni)*Pi + 0.187) .^ Ii .* (ones(ni)*Sigma - 0.755) .^ Ji) * 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

    fun = sum(ni .* (ones(ni)*Pi + 0.298) .^ Ii .* (ones(ni)*Sigma - 0.816) .^ Ji) * 0.0088;
end
endfunction

//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
function fun = p3_hs(h, s)

if s < 4.41202148223476 then
    //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

    fun = sum(ni .* (ones(ni)*eta - 1.01) .^ Ii .* (ones(ni)*Sigma - 0.75) .^ Ji) * 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

    fun = 16.6 / sum(ni .* (ones(ni)*eta - 0.681) .^ Ii .* (ones(ni)*Sigma - 0.792) .^ Ji);
end
endfunction

//Not avalible with if 97
//Solve function T3_ph-T=0 with half interval method.
function fun = h3_pT(p, T)

//ver2.6 Start corrected bug
if p < 22.06395 then   //Bellow tripple point
  Ts = T4_p(p);    //Saturation temperature
  if T <= Ts then     //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 then
        High_Bound = hs;
    else
        Low_Bound = hs;
    end
end
fun = hs;
endfunction

  //Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  //Solve with half interval method
function fun = T3_prho(p, rho)

  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 then
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
  fun = Ts;
endfunction  


//***********************************************************************************************************
//*2.4 functions for region 4

//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
function fun = p4_T(T)
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;
fun = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4;
endfunction

//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
function fun = T4_p(p)
bet_a = p ^ 0.25;
E = bet_a ^ 2 - 17.073846940092 * bet_a + 14.91510861353;
f = 1167.0521452767 * bet_a ^ 2 + 12020.82470247 * bet_a - 4823.2657361591;
G = -724213.16703206 * bet_a ^ 2 - 3232555.0322333 * bet_a + 405113.40542057;
D = 2 * G / (-f - (f ^ 2 - 4 * E * G) ^ 0.5);
fun = (650.17534844798 + D - ((650.17534844798 + D) ^ 2 - 4 * (-0.23855557567849 + 650.17534844798 * D)) ^ 0.5) / 2;
endfunction

//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
//See picture page 14
function fun = h4_s(s)

if (s > -0.0001545495919 & s <= 3.77828134) then
    //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

    fun = sum(ni .* (ones(ni)*Sigma - 1.09) .^ Ii .* (ones(ni)*Sigma + 0.0000366) .^ Ji) * 1700;

elseif (s > 3.77828134 & s <= 4.41202148223476 ) then
    //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

    fun = sum(ni .* (ones(ni)*Sigma - 1.09) .^ Ii .* (ones(ni)*Sigma + 0.0000366) .^ Ji) * 1700;

elseif (s > 4.41202148223476 & s <= 5.85 ) then
    //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

    fun = sum(ni .* (ones(ni)*Sigma - 1.02) .^ Ii .* (ones(ni)*Sigma - 0.726) .^ Ji) ^ 4 * 2800;

elseif (s > 5.85 & s < 9.155759395) then
    //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

    fun = exp(sum(ni .* (ones(ni)*1/Sigma1 - 0.513) .^ Ii .* (ones(ni)*Sigma2 - 0.524) .^ Ji)) * 2800;

else
    fun = -99999;
end
endfunction

//Uses h4_s and p_hs for the diffrent regions to determine p4_s
function fun = p4_s(s)
h_sat = h4_s(s);
if (s > -0.0001545495919 & s <= 3.77828134) then
    fun = p1_hs(h_sat, s);
elseif (s > 3.77828134 & s <= 5.210887663) then
    fun = p3_hs(h_sat, s);
elseif (s > 5.210887663 & s < 9.155759395) then
    fun = p2_hs(h_sat, s);
else
    fun = -99999;
end
endfunction

function fun = h4L_p(p)
if (p > 0.000611657 & p < 22.06395) then
    Ts = T4_p(p);
    if p < 16.529 then
        fun = 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 then
                High_Bound = hs;
            else
                Low_Bound = hs;
            end
        end
        
        fun = hs;
    end
else
    fun = -99999;
end
endfunction


function fun = h4V_p(p)
if (p > 0.000611657 & p < 22.06395) then
    Ts = T4_p(p);
    if p < 16.529 then
        fun = 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 then
                High_Bound = hs;
            else
                Low_Bound = hs;
            end
        end
        fun = hs;
    end
else
    fun = -99999;
end
endfunction


function fun = 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 then
    fun = 1;
elseif h < h4L then
    fun = 0;
else
    fun = (h - h4L) / (h4v - h4L);
end
endfunction


function fun = x4_ps(p, s)
if p < 16.529 then
    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 then
    fun = 0;
elseif s > ssv then
    fun = 1;
else
    fun = (s - ssL) / (ssv - ssL);
end
endfunction

//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.
function fun = T4_hs(h, s)
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) then
    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

    fun = sum( ni .* (ones(ni)*eta - 0.119) .^ Ii * (ones(ni)*Sigma - 1.07) .^ Ji) * 550;
else
    //function psat_h
    if (s > -0.0001545495919 & s <= 3.77828134) then
        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 then
                High_Bound = PL;
            else
                Low_Bound = PL;
            end
        end
    elseif (s > 3.77828134 & s <= 4.41202148223476) then 
        PL = p3sat_h(h);
    elseif (s > 4.41202148223476 & s <= 5.210887663) then 
        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)
        p = (Low_Bound + High_Bound) / 2;
        //Calculate s4_ph
        Ts = T4_p(p);
        xs = x4_ph(p, h);
        if p < 16.529 then
            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 then
            High_Bound = p;
        else
            Low_Bound = p;
        end
    end
    fun = T4_p(p);
end
endfunction

//***********************************************************************************************************
//*2.5 functions for region 5

//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
function fun = h5_pT(p, T)
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
gamma0_tau = sum(ni0 .* Ji0 .* (ones(ni0)*tau) .^ (Ji0 - 1));
gammar_tau = sum(nir .* (ones(nir)*Pi) .^ Iir .* Jir .* (ones(nir)*tau) .^ (Jir - 1));

fun = R * T * tau * (gamma0_tau + gammar_tau);
endfunction

//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
function fun = v5_pT(p, T)
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
gammar_pi = sum(nir .* Iir .* (ones(nir)*Pi) .^ (Iir - 1) .* (ones(nir)*tau) .^ Jir);

fun = R * T / p * Pi * (gamma0_pi + gammar_pi) / 1000;
endfunction

//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
function fun = u5_pT(p, T)
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

gamma0_tau = sum(ni0 .* Ji0 .* (ones(ni0)*tau) .^ (Ji0 - 1));
gammar_pi = sum(nir .* Iir .* (ones(nir)*Pi) .^ (Iir - 1) .* (ones(nir)*tau) .^ Jir);
gammar_tau = sum(nir .* (ones(nir)*Pi) .^ Iir .* Jir .* (ones(nir)*tau) .^ (Jir - 1));

fun = R * T * (tau * (gamma0_tau + gammar_tau) - Pi * (gamma0_pi + gammar_pi));

endfunction

//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
function fun = Cp5_pT(p, T)
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
gamma0_tautau = sum(ni0 .* Ji0 .* (Ji0 - 1) .* (ones(ni0)*tau) .^ (Ji0 - 2));
gammar_tautau = sum(nir .* (ones(nir)*Pi) .^ Iir .* Jir .* (Jir - 1) .* (ones(nir)*tau) .^ (Jir - 2));

fun = -R * tau ^ 2 * (gamma0_tautau + gammar_tautau);
endfunction

//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
function fun = s5_pT(p, T)
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

gamma0 = log(Pi)+sum(ni0 .* (ones(ni0)*tau) .^ Ji0);
gamma0_tau = sum(ni0 .* Ji0 .* (ones(ni0)*tau) .^ (Ji0 - 1));
gammar = sum(nir .* (ones(nir)*Pi) .^ Iir .* (ones(nir)*tau) .^ Jir);
gammar_tau = sum(nir .* (ones(nir)*Pi) .^ Iir .* Jir .* (ones(nir)*tau) .^ (Jir - 1));

fun = R * (tau * (gamma0_tau + gammar_tau) - (gamma0 + gammar));
endfunction

//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
function fun = Cv5_pT(p, T)
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

gamma0_tautau = sum(ni0 .* (Ji0 - 1) .* Ji0 .* (ones(ni0)*tau) .^ (Ji0 - 2));
gammar_pi = sum(nir .* Iir .* (ones(nir)*Pi) .^ (Iir - 1) .* (ones(nir)*tau) .^ Jir);
gammar_pitau = sum(nir .* Iir .* (ones(nir)*Pi) .^ (Iir - 1) .* Jir .* (ones(nir)*tau) .^ (Jir - 1));
gammar_pipi = sum(nir .* Iir .* (Iir - 1) .* (ones(nir)*Pi) .^ (Iir - 2) .* (ones(nir)*tau) .^ Jir);
gammar_tautau = sum(nir .* (ones(nir)*Pi) .^ Iir .* Jir .* (Jir - 1) .* (ones(nir)*tau) .^ (Jir - 2));

fun = R * (-(tau ^ 2 *(gamma0_tautau + gammar_tautau)) - (1 + Pi * gammar_pi - tau * Pi * gammar_pitau)^2 / (1 - Pi ^ 2 * gammar_pipi));
endfunction


//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
function fun = w5_pT(p, T)
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

gamma0_tautau = sum(ni0 .* (Ji0 - 1) .* Ji0 .* (ones(ni0)*tau) .^ (Ji0 - 2));
gammar_pi = sum(nir .* Iir .* (ones(nir)*Pi) .^ (Iir - 1) .* (ones(nir)*tau) .^ Jir);
gammar_pitau = sum(nir .* Iir .* (ones(nir)*Pi) .^ (Iir - 1) .* Jir .* (ones(nir)*tau) .^ (Jir - 1));
gammar_pipi = sum(nir .* Iir .* (Iir - 1) .* (ones(nir)*Pi) .^ (Iir - 2) .* (ones(nir)*tau) .^ Jir);
gammar_tautau = sum(nir .* (ones(nir)*Pi) .^ Iir .* Jir .* (Jir - 1) .* (ones(nir)*tau) .^ (Jir - 2));

fun = (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;
endfunction

function fun = 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 then
        High_Bound = Ts;
    else
        Low_Bound = Ts;
    end
end
fun = Ts;
endfunction

function fun = 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 then
        High_Bound = Ts;
    else
        Low_Bound = Ts;
    end
end
fun = Ts;
endfunction


function fun=T5_prho(p,rho)
  //Solve by iteration. Observe that for 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 then
      High_Bound = Ts;
    else
      Low_Bound = Ts;
    end
  end
    fun = Ts;
endfunction

//***********************************************************************************************************
//*3 Region Selection
//***********************************************************************************************************
//*3.1 Regions as a function of pT
function fun = region_pT(p, T)
if (T > 1073.15 & p < 10 & T < 2273.15 & p > 0.000611) then 
    fun = 5;
elseif (T <= 1073.15 & T > 273.15 & p <= 100 & p > 0.000611) then 
    if T > 623.15 then
        if p > B23p_T(T) then 
            fun = 3;
            if T < 647.096 then
                ps = p4_T(T);
                if abs(p - ps) < 0.00001 then
                    fun = 4;
                end
            end
        else
            fun = 2;
        end
    else
        ps = p4_T(T);
        if abs(p - ps) < 0.00001 then
            fun = 4;
        elseif p > ps then
            fun = 1;
        else
            fun = 2;
        end
    end
else
    fun = 0; //**Error, Outside valid area
end
endfunction

//***********************************************************************************************************
//*3.2 Regions as a function of ph
function fun = region_ph(  p,   h)
//Check if outside pressure limits
if (p < 0.000611657 | p > 100) then 
    fun = 0;
    return
end

//Check if outside low h.
if h < 0.963 * p + 2.2 then //Linear adaption to h1_pt()+2 to speed up calcualations.
    if h < h1_pT(p, 273.15) then 
        fun = 0;
        return
    end
end

if p < 16.5292 then //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 then //if approximate is not god enough use real function
        hL = h1_pT(p, Ts);
    end
    if h <= hL then
        fun = 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 then //if approximate is not god enough use real function
        hV = h2_pT(p, Ts);
    end
    if h < hV then
        fun = 4;
        return
    end
    //Check upper limit of region 2 Quick Test
    if h < 4000 then
        fun = 2;
        return
    end
    //Check region 2 (Real value)
    h_45 = h2_pT(p, 1073.15);
    if h <= h_45 then
        fun = 2;
        return
    end
    //Check region 5
    if p > 10 then
        fun = 0;
        return
    end
    h_5u = h5_pT(p, 2273.15);
    if h < h_5u then
        fun = 5;
        return
    end
    fun = 0;
    return
else //for p>16.5292
    //Check if in region1
    if h < h1_pT(p, 623.15) then
        fun = 1;
        return
    end
    //Check if in region 3 or 4 (Bellow Reg 2)
    if h < h2_pT(p, B23T_p(p)) then
        //Region 3 or 4
        if p > p3sat_h(h) then
            fun = 3;
            return
        else
            fun = 4;
            return
        end
    end
    //Check if region 2
    if h < h2_pT(p, 1073.15) then
        fun = 2;
        return
    end
end
fun = 0;
endfunction

//***********************************************************************************************************
//*3.3 Regions as a function of ps
function fun = region_ps(  p,   s)
if (p < 0.000611657 | p > 100 | s < 0 | s > s5_pT(p, 2273.15)) then 
    fun = 0;
    return
end

//Check region 5
if s > s2_pT(p, 1073.15) then
    if p <= 10 then
        fun = 5;
        return
    else
        fun = 0;
        return
    end
end

//Check region 2
if p > 16.529 then
    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 then
    fun = 2;
    return
end

//Check region 3
ss = s1_pT(p, 623.15);
if (p > 16.529 & s > ss) then 
    if p > p3sat_s(s) then
        fun = 3;
        return
    else
        fun = 4;
        return
    end
end

//Check region 4 (Not inside region 3)
if (p < 16.529 & s > s1_pT(p, T4_p(p))) then 
    fun = 4;
    return
end

//Check region 1
if (p > 0.000611657 & s > s1_pT(p, 273.15)) then 
    fun = 1;
    return
end
fun = 1;
endfunction

//***********************************************************************************************************
//*3.4 Regions as a function of hs
function fun = region_hs(  h,   s)
if s < -0.0001545495919 then
    fun = 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 ) then
    fun = 0;
    return
end

//******Kolla 1 eller 4. (+liten bit över B13)
if (s >= -0.0001545495919 & s <= 3.77828134 ) then
    if h < h4_s(s) then
        fun = 4;
        return
    elseif s < 3.397782955 then //100MPa line is limiting
        TMax = T1_ps(100, s);
        hMax = h1_pT(100, TMax);
        if h < hMax then
            fun = 1;
            return
        else
            fun = 0;
            return
        end
    else //The point is either in region 4,1,3. Check B23
        hB = hB13_s(s);
        if h < hB then
            fun = 1;
            return
        end
        TMax = T3_ps(100, s);
        vmax = v3_ps(100, s);
        hMax = h3_rhoT(1 / vmax, TMax);
        if h < hMax then
            fun = 3;
            return
        else
            fun = 0;
            return
        end
    end
end

//******Kolla region 2 eller 4. (Övre delen av område b23-> max)
if (s >= 5.260578707 & s <= 11.9212156897728) then 
    if s > 9.155759395 then //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 ) then
            fun = 2;
            return
        else
            fun = 0;
            return
        end
    end
    
    hV = h4_s(s);
    
    if h < hV  then //Region 4. Under region 3.
        fun = 4;
        return
    end
    if s < 6.04048367171238 then
        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 then //Region 2. Över region 4.
        fun = 2;
        return
    else
        fun = 0;
        return
    end
end

//Kolla region 3 eller 4. Under kritiska punkten.
if (s >= 3.77828134 & s <= 4.41202148223476) then 
    hL = h4_s(s);
    if h < hL then
        fun = 4;
        return
    end
    TMax = T3_ps(100, s);
    vmax = v3_ps(100, s);
    hMax = h3_rhoT(1 / vmax, TMax);
    if h < hMax then
        fun = 3;
        return
    else
        fun = 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) then 
    hV = h4_s(s);
    if h < hV then
        fun = 4;
        return
    end
    //Kolla om vi är under b23 giltighetsområde.
    if s <= 5.048096828 then
        TMax = T3_ps(100, s);
        vmax = v3_ps(100, s);
        hMax = h3_rhoT(1 / vmax, TMax);
        if h < hMax then
            fun = 3;
            return
        else
            fun = 0;
            return
        end
    else //Inom området för B23 i s led.
        if (h > 2812.942061) then //Ovanför B23 i h_led
            if s > 5.09796573397125 then
                TMax = T2_ps(100, s);
                hMax = h2_pT(100, TMax);
                if h < hMax then
                    fun = 2;
                    return
                else
                    fun = 0;
                    return
                end
            else
                fun = 0;
                return
            end
        end
        if h < 2563.592004 then    //Nedanför B23 i h_led men vi har redan kollat ovanför hV2c3b
            fun = 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 then
            fun = 3;
            return
        else
            fun = 2;
            return
        end
    end
end
fun = 0;
endfunction

//***********************************************************************************************************
//*3.5 Regions as a function of p and rho
function fun = Region_prho(p,rho)
v = 1 / rho;
if (p < 0.000611657 | p > 100) then
    fun = 0;
    return
end
if p < 16.5292 then //Bellow region 3, Check region 1,4,2
    if v < v1_pT(p, 273.15) then//Observe that this is not actually min of v. Not valid Water of 4°C is ligther.
        fun = 0;
        return
    end
    if v <= v1_pT(p, T4_p(p)) then
        fun = 1;
        return
    end
    if v < v2_pT(p, T4_p(p)) then
        fun = 4;
        return
    end
    if v <= v2_pT(p, 1073.15) then
        fun = 2;
        return
    end 
    if p > 10 then //Above region 5
        fun = 0;
        return
    end 
    if v <= v5_pT(p, 2073.15) then
        fun = 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) then //Observe that this is not actually min of v. Not valid Water of 4°C is ligther.
        fun = 0;
        return
    end
    if v < v1_pT(p, 623.15) then
        fun = 1;
        return
    end
    //Check if in region 3 or 4 (Bellow Reg 2)
    if v < v2_pT(p, B23T_p(p)) then
        //Region 3 or 4
        if p > 22.064 then //Above region 4
            fun = 3;
            return
        end
        if (v < v3_ph(p, h4L_p(p)) | v > v3_ph(p, h4V_p(p))) then //Uses iteration!!
            fun = 3;
            return
        else
            fun = 4;
            return
        end 
    end 
    //Check if region 2
    if v < v2_pT(p, 1073.15) then
        fun = 2;
        return
    end
end 
fun = 0;
endfunction


//***********************************************************************************************************
//*4 Region Borders
//***********************************************************************************************************
//***********************************************************************************************************
//*4.1 Boundary between region 2 and 3.
//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

function fun = B23p_T(T)

fun = 348.05185628969 - 1.1671859879975 * T + 1.0192970039326E-03 * T ^ 2;

endfunction

//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
function fun = B23T_p(p)

fun = 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ^ 0.5;

endfunction

//***********************************************************************************************************
//*4.2 Region 3. pSat_h  & pSat_s

//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
function fun = p3sat_h(h)
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

fun = sum(ni .* (ones(ni)*hs - 1.02) .^ Ii .* (ones(ni)*hs - 0.608) .^ Ji) * 22;
endfunction

function fun = 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

fun = sum(ni .* (ones(ni)*Sigma - 1.03) .^ Ii .* (ones(ni)*Sigma - 0.699) .^ Ji) * 22;
endfunction

//***********************************************************************************************************
//4.3 Region boundary 1to3  & 3to2 as a functions of s

//Supplementary Release on Backward Equations ( ) , p h s for Region 3,
//Chapter 4.5 page 23.
function fun = hB13_s(s)
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

fun = sum(ni .* (ones(ni)*Sigma - 0.884) .^ Ii .* (ones(ni)*Sigma - 0.864) .^ Ji) * 1700;
endfunction

//Supplementary Release on Backward Equations ( ) , p h s for Region 3,
//Chapter 4.6 page 25.
function fun = TB23_hs(h, s)
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

fun = sum(ni .* (ones(ni)*eta - 0.727) .^ Ii .* (ones(ni)*Sigma - 0.864) .^ Ji) * 900;
endfunction

//***********************************************************************************************************
//*5 Transport properties
//***********************************************************************************************************
//*5.1 Viscosity (IAPWS formulation 1985, Revised 2003)
//***********************************************************************************************************
function fun = 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.
select region_pT(p, T)
case 1 then
    rho = 1 / v1_pT(p, T);
case 2 then
    rho = 1 / v2_pT(p, T);
case 3 then
    hs = h3_pT(p, T);
    rho = 1 / v3_ph(p, hs);
case 4 then
    rho = %nan;
case 5 then
    rho = 1 / v5_pT(p, T);
else
    fun = %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) then 
    fun = %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

Sum = sum(h0 .* (ones(h0)*1 / Ts - 1) .^ [0:5] + h1 .* (ones(h1)*1 / Ts - 1) .^ [0:5] .* (ones(h1)*rhos - 1) .^ 1 ..
                                                                               + h2 .* (ones(h2)*1 / Ts - 1) .^ [0:5] .* (ones(h2)*rhos - 1) .^ 2 ..
                                                                               + h3 .* (ones(h3)*1 / Ts - 1) .^ [0:5] .* (ones(h3)*rhos - 1) .^ 3 ..
                                                                               + h4 .* (ones(h4)*1 / Ts - 1) .^ [0:5] .* (ones(h4)*rhos - 1) .^ 4 ..
                                                                               + h5 .* (ones(h5)*1 / Ts - 1) .^ [0:5] .* (ones(h5)*rhos - 1) .^ 5 ..
                                                                               + h6 .* (ones(h6)*1 / Ts - 1) .^ [0:5] .* (ones(h6)*rhos - 1) .^ 6);
my1 = exp(rhos * Sum);
mys = my0 * my1;
fun = mys * 0.000055071;
endfunction

function fun = 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.
select region_ph(p, h)
case 1 then
    Ts = T1_ph(p, h);
    T = Ts;
    rho = 1 / v1_pT(p, Ts);
case 2 then
    Ts = T2_ph(p, h);
    T = Ts;
    rho = 1 / v2_pT(p, Ts);
case 3 then
    rho = 1 / v3_ph(p, h);
    T = T3_ph(p, h);
case 4 then
    xs = x4_ph(p, h);
    if p < 16.529 then
        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 then
    Ts = T5_ph(p, h);
    T = Ts;
    rho = 1 / v5_pT(p, Ts);
else
    fun = %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 ) then
    fun = %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

Sum = sum(h0 .* (ones(h0)*1 / Ts - 1) .^ [0:5] + h1 .* (ones(h1)*1 / Ts - 1) .^ [0:5] .* (ones(h1)*rhos - 1) .^ 1 ..
                                                                               + h2 .* (ones(h2)*1 / Ts - 1) .^ [0:5] .* (ones(h2)*rhos - 1) .^ 2 ..
                                                                               + h3 .* (ones(h3)*1 / Ts - 1) .^ [0:5] .* (ones(h3)*rhos - 1) .^ 3 ..
                                                                               + h4 .* (ones(h4)*1 / Ts - 1) .^ [0:5] .* (ones(h4)*rhos - 1) .^ 4 ..
                                                                               + h5 .* (ones(h5)*1 / Ts - 1) .^ [0:5] .* (ones(h5)*rhos - 1) .^ 5 ..
                                                                               + h6 .* (ones(h6)*1 / Ts - 1) .^ [0:5] .* (ones(h6)*rhos - 1) .^ 6);
                                                                               
my1 = exp(rhos * Sum);
mys = my0 * my1;
fun = mys * 0.000055071;
endfunction

//***********************************************************************************************************
//*5.2 Thermal Conductivity (IAPWS formulation 1985)

//Revised release on the IAPWS formulation 1985 for the Thermal Conductivity of ordinary water
//IAPWS September 1998
//Page 8
function fun = tc_ptrho(  p,   T,   rho)

//ver2.6 Start corrected bug
 if T < 273.15 then
   fun = %nan; //Out of range of validity (para. B4)
   return
 elseif T < 500 + 273.15 then
   if p > 100 then
     fun = %nan; //Out of range of validity (para. B4)
     return
   end
 elseif T <= 650 + 273.15 then
   if p > 70 then
     fun = %nan; //Out of range of validity (para. B4)
     return
   end
 else T <= 800 + 273.15 then
   if p > 40 then
      fun = %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 then
    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);
fun = tc0 + tc1 + tc2;
endfunction

//***********************************************************************************************************
//5.3 Surface Tension

//IAPWS Release on Surface Tension of Ordinary Water Substance,
//September 1994
function fun = Surface_Tension_T(  T)
tc = 647.096; //K
B = 0.2358;    //N/m
bb = -0.625;
my = 1.256;
if (T < 0.01 | T > tc) then 
    fun = %nan; //"Out of valid region"
    return
end
tau = 1 - T / tc;
fun = B * tau ^ my * (1 + bb * tau);
endfunction

//***********************************************************************************************************
//*6 Units                                                                                      *
//***********************************************************************************************************
function fun = toSIunit_p( Ins )
//Translate bar to MPa
fun = Ins / 10;
endfunction

function fun = fromSIunit_p( Ins )
//Translate bar to MPa
fun = Ins * 10;
endfunction

function fun = toSIunit_T( Ins )
//Translate degC to Kelvon
fun = Ins + 273.15;
endfunction

function fun = fromSIunit_T( Ins )
//Translate Kelvin to degC
fun = Ins - 273.15;
endfunction

function fun = toSIunit_h( Ins )
fun = Ins;
endfunction

function  fun = fromSIunit_h( Ins )
fun = Ins;
endfunction

function fun = toSIunit_v( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_v( Ins )
fun = Ins;
endfunction

function fun = toSIunit_s( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_s( Ins )
fun = Ins;
endfunction

function fun = toSIunit_u( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_u( Ins )
fun = Ins;
endfunction

function fun = toSIunit_Cp( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_Cp( Ins )
fun = Ins;
endfunction

function fun = toSIunit_Cv( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_Cv( Ins )
fun = Ins;
endfunction

function fun = toSIunit_w( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_w( Ins )
fun = Ins;
endfunction

function fun = toSIunit_tc( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_tc( Ins )
fun = Ins;
endfunction

function fun = toSIunit_st( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_st( Ins )
fun = Ins;
endfunction

function fun = toSIunit_x( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_x( Ins )
fun = Ins;
endfunction

function fun = toSIunit_vx( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_vx( Ins )
fun = Ins;
endfunction

function fun = toSIunit_my( Ins )
fun = Ins;
endfunction

function fun = fromSIunit_my( Ins )
fun = Ins;
endfunction

//***********************************************************************************************************
//*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
        execstr('R1(j,i)='+Fun(j)+'('+string(p(i))+','+string(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
        execstr('R2(j,i)='+Fun(j)+'('+string(p(i))+','+string(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
        execstr('R3(j,i)='+Fun(j)+'('+string(rho(i))+','+string(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;
R4=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];
R4=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];
R4=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
        execstr('R5(j,i)='+Fun(j)+'('+string(p(i))+','+string(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:size(fun,'c')
    
    T=['XSteam(''',fun(i),''',',In1(i),',',In2(i),')'];
    execstr('Res='+'XSteam'+'('+fun(i)+','+In1(i)+','+In2(i)+')');
    Error=(Res-evstr(Control(i)))/evstr(Control(i));
    Check(i,:)=[T,'=',string(Res),' - (Control)',Control(i),'=',string(Error)];
    if Error>1E-5
        err=err+1;
        error('To large error')
    end   
end

err = Check;

endfunction
