SHARE
TWEET

Untitled

a guest Apr 10th, 2019 139 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. Option Explicit
  2. 'Niezbędne informacje:
  3. 'A/ Jednostki:
  4. ' - entalpia, h kJ/kg
  5. ' - entropia, s kJ/kgK
  6. ' - temperatura, t C
  7. ' - ciśnienie, p bar
  8. '
  9. 'B/ Funkcje:
  10. ' - entalpia w funkcji ciśnienia i temperatury: h_pT(p,T)
  11. ' - entropia w funkcji ciśnienia i temperatury: s_pT(p,T)
  12. ' - stopień suchości pary w funkcji ciśnienia i entropii: x_ps(p,s)
  13. ' - entalpia w funkcji ciśnienia i entropii: h_ps(p,s)
  14. ' - entalpia pary nasyconej suchej w funkcji ciśnienia: hV_p(p)
  15. ' - entalpia wody nasyconej w funkcji temperatury: hL_t(t)
  16. ' - entalpia wody nasyconej w funkcji ciśnienia: hL_p(p)
  17.  
  18.  
  19.  
  20. '***********************************************************************************************************
  21. '* Water and steam properties according to IAPWS IF-97                                                     *
  22. '* By Magnus Holmgren, www.x-eng.com                                                                       *
  23. '* The steam tables are free and provided as is.                                                           *
  24. '* We take no responsibilities for any errors in the code or damage thereby.                               *
  25. '* You are free to use, modify and distribute the code as long as authorship is properly acknowledged.     *
  26. '* Please notify me at magnus@x-eng.com if the code is used in commercial applications                     *
  27. '***********************************************************************************************************'
  28. '
  29. ' The code is also avalibale for matlab at www.x-eng.com
  30. '
  31. '*Contents.
  32. '*1 Calling functions
  33. '*1.1
  34. '*1.2 Temperature (T)
  35. '*1.3 Pressure (p)
  36. '*1.4 Enthalpy (h)
  37. '*1.5 Specific Volume (v)
  38. '*1.6 Density (rho)
  39. '*1.7 Specific entropy (s)
  40. '*1.8 Specific internal energy (u)
  41. '*1.9 Specific isobaric heat capacity (Cp)
  42. '*1.10 Specific isochoric heat capacity (Cv)
  43. '*1.11 Speed of sound
  44. '*1.12 Viscosity
  45. '*1.13 Prandtl
  46. '*1.14 Kappa
  47. '*1.15 Surface tension
  48. '*1.16 Heat conductivity
  49. '*1.17 Vapour fraction
  50. '*1.18 Vapour Volume Fraction
  51. '
  52. '*2 IAPWS IF 97 Calling functions
  53. '*2.1 Functions for region 1
  54. '*2.2 Functions for region 2
  55. '*2.3 Functions for region 3
  56. '*2.4 Functions for region 4
  57. '*2.5 Functions for region 5
  58. '
  59. '*3 Region Selection
  60. '*3.1 Regions as a function of pT
  61. '*3.2 Regions as a function of ph
  62. '*3.3 Regions as a function of ps
  63. '*3.4 Regions as a function of hs
  64. '*3.5 Regions as a function of p and rho
  65. '
  66. '4 Region Borders
  67. '4.1 Boundary between region 1 and 3.
  68. '4.2 Region 3. pSat_h and pSat_s
  69. '4.3 Region boundary 1to3 and 3to2 as a functions of s
  70. '
  71. '5 Transport properties
  72. '5.1 Viscosity (IAPWS formulation 1985)
  73. '5.2 Thermal Conductivity (IAPWS formulation 1985)
  74. '5.3 Surface Tension
  75. '
  76. '6 Units
  77.  
  78.  
  79. '***********************************************************************************************************
  80. '*1 Calling functions                                                                                      *
  81. '***********************************************************************************************************
  82.  
  83. '***********************************************************************************************************
  84. '*1.1
  85.  
  86.  
  87. '***********************************************************************************************************
  88. '*1.2 Temperature
  89. Function Tsat_p(ByVal p As Double) As Double
  90.  p = toSIunit_p(p)
  91.  If p >= 0.000611657 And p <= 22.06395 + 0.001 Then '0.001 Added to enable the tripple point.
  92.   Tsat_p = fromSIunit_T(T4_p(p))
  93.  Else
  94.    Tsat_p = CVErr(xlErrValue)
  95.  End If
  96. End Function
  97. Function Tsat_s(ByVal s As Double) As Double
  98.  s = toSIunit_s(s)
  99.  If s > -0.0001545495919 And s < 9.155759395 Then
  100.    Tsat_s = fromSIunit_T(T4_p(p4_s(s)))
  101.  Else
  102.    Tsat_s = CVErr(xlErrValue)
  103.  End If
  104. End Function
  105.  
  106. Function T_ph(ByVal p As Double, ByVal h As Double) As Double
  107.  p = toSIunit_p(p)
  108.  h = toSIunit_h(h)
  109.  Select Case region_ph(p, h)
  110.  Case 1
  111.    T_ph = fromSIunit_T(T1_ph(p, h))
  112.  Case 2
  113.    T_ph = fromSIunit_T(T2_ph(p, h))
  114.  Case 3
  115.    T_ph = fromSIunit_T(T3_ph(p, h))
  116.  Case 4
  117.    T_ph = fromSIunit_T(T4_p(p))
  118.  Case 5
  119.    T_ph = fromSIunit_T(T5_ph(p, h))
  120.  Case Else
  121.   T_ph = CVErr(xlErrValue)
  122.  End Select
  123. End Function
  124. Function T_ps(ByVal p As Double, ByVal s As Double) As Double
  125.  p = toSIunit_p(p)
  126.  s = toSIunit_s(s)
  127.  Select Case region_ps(p, s)
  128.  Case 1
  129.    T_ps = fromSIunit_T(T1_ps(p, s))
  130.  Case 2
  131.    T_ps = fromSIunit_T(T2_ps(p, s))
  132.  Case 3
  133.    T_ps = fromSIunit_T(T3_ps(p, s))
  134.  Case 4
  135.    T_ps = fromSIunit_T(T4_p(p))
  136.  Case 5
  137.    T_ps = fromSIunit_T(T5_ps(p, s))
  138.  Case Else
  139.   T_ps = CVErr(xlErrValue)
  140.  End Select
  141. End Function
  142. Function T_hs(ByVal h As Double, ByVal s As Double) As Double
  143.  h = toSIunit_h(h)
  144.  s = toSIunit_s(s)
  145.  Select Case Region_hs(h, s)
  146.  Case 1
  147.    T_hs = fromSIunit_T(T1_ph(p1_hs(h, s), h))
  148.  Case 2
  149.    T_hs = fromSIunit_T(T2_ph(p2_hs(h, s), h))
  150.  Case 3
  151.    T_hs = fromSIunit_T(T3_ph(p3_hs(h, s), h))
  152.  Case 4
  153.    T_hs = fromSIunit_T(T4_hs(h, s))
  154.  Case 5
  155.    T_hs = CVErr(xlErrValue) 'Functions of hs is not implemented in region 5
  156. Case Else
  157.    T_hs = CVErr(xlErrValue)
  158.  End Select
  159. End Function
  160. '***********************************************************************************************************
  161. '*1.3 Pressure (p)
  162. Function psat_T(ByVal t As Double) As Double
  163.  t = toSIunit_T(t)
  164.  If t <= 647.096 And t > 273.15 Then
  165.    psat_T = fromSIunit_p(p4_T(t))
  166.  Else
  167.    psat_T = CVErr(xlErrValue)
  168.  End If
  169. End Function
  170. Function psat_s(ByVal s As Double) As Double
  171.  s = toSIunit_s(s)
  172.  If s > -0.0001545495919 And s < 9.155759395 Then
  173.    psat_s = fromSIunit_p(p4_s(s))
  174.  Else
  175.    psat_s = CVErr(xlErrValue)
  176.  End If
  177. End Function
  178. Function p_hs(ByVal h As Double, ByVal s As Double) As Double
  179.  h = toSIunit_h(h)
  180.  s = toSIunit_s(s)
  181.  Select Case Region_hs(h, s)
  182.  Case 1
  183.    p_hs = fromSIunit_p(p1_hs(h, s))
  184.  Case 2
  185.    p_hs = fromSIunit_p(p2_hs(h, s))
  186.  Case 3
  187.    p_hs = fromSIunit_p(p3_hs(h, s))
  188.  Case 4
  189.    p_hs = fromSIunit_p(p4_T(T4_hs(h, s)))
  190.  Case 5
  191.    p_hs = CVErr(xlErrValue) 'Functions of hs is not implemented in region 5
  192. Case Else
  193.   p_hs = CVErr(xlErrValue)
  194.  End Select
  195. End Function
  196. Function p_hrho(ByVal h As Double, ByVal rho As Double) As Double
  197. '*Not valid for water or sumpercritical since water rho does not change very much with p.
  198. '*Uses iteration to find p.
  199.  Dim High_Bound As Double
  200.   Dim Low_Bound As Double
  201.   Dim p As Double
  202.   Dim rhos As Double
  203.   High_Bound = fromSIunit_p(100)
  204.   Low_Bound = fromSIunit_p(0.000611657)
  205.   p = fromSIunit_p(10)
  206.   rhos = 1 / v_ph(p, h)
  207.   Do While Abs(rho - rhos) > 0.0000001
  208.     rhos = 1 / v_ph(p, h)
  209.     If rhos >= rho Then
  210.       High_Bound = p
  211.     Else
  212.       Low_Bound = p
  213.     End If
  214.     p = (Low_Bound + High_Bound) / 2
  215.     Loop
  216.     p_hrho = p
  217. End Function
  218.  
  219.  
  220. '***********************************************************************************************************
  221. '*1.4 Enthalpy (h)
  222. Function hV_p(ByVal p As Double) As Double
  223.  p = toSIunit_p(p)
  224.  If p > 0.000611657 And p < 22.06395 Then
  225.    hV_p = fromSIunit_h(h4V_p(p))
  226.  Else
  227.    hV_p = CVErr(xlErrValue)
  228.  End If
  229. End Function
  230. Function hL_p(ByVal p As Double) As Double
  231.  p = toSIunit_p(p)
  232.  If p > 0.000611657 And p < 22.06395 Then
  233.    hL_p = fromSIunit_h(h4L_p(p))
  234.  Else
  235.    hL_p = CVErr(xlErrValue)
  236.  End If
  237. End Function
  238. Function hV_T(ByVal t As Double) As Double
  239.  t = toSIunit_T(t)
  240.  If t > 273.15 And t < 647.096 Then
  241.   hV_T = fromSIunit_h(h4V_p(p4_T(t)))
  242.  Else
  243.   hV_T = CVErr(xlErrValue)
  244.  End If
  245. End Function
  246. Function hL_T(ByVal t As Double) As Double
  247.  t = toSIunit_T(t)
  248.  If t > 273.15 And t < 647.096 Then
  249.   hL_T = fromSIunit_h(h4L_p(p4_T(t)))
  250. Else
  251.   hL_T = CVErr(xlErrValue)
  252.  End If
  253. End Function
  254.  
  255. Function h_pT(ByVal p As Double, ByVal t As Double) As Double
  256.  p = toSIunit_p(p)
  257.  t = toSIunit_T(t)
  258.  Select Case region_pT(p, t)
  259.  Case 1
  260.    h_pT = fromSIunit_h(h1_pT(p, t))
  261.  Case 2
  262.    h_pT = fromSIunit_h(h2_pT(p, t))
  263.  Case 3
  264.    h_pT = fromSIunit_h(h3_pT(p, t))
  265.  Case 4
  266.    h_pT = CVErr(xlErrValue)
  267.  Case 5
  268.    h_pT = fromSIunit_h(h5_pT(p, t))
  269.  Case Else
  270.   h_pT = CVErr(xlErrValue)
  271.  End Select
  272. End Function
  273. Function h_ps(ByVal p As Double, ByVal s As Double) As Double
  274.  Dim xs As Double
  275.  p = toSIunit_p(p)
  276.  s = toSIunit_s(s)
  277.  Select Case region_ps(p, s)
  278.  Case 1
  279.    h_ps = fromSIunit_h(h1_pT(p, T1_ps(p, s)))
  280.  Case 2
  281.    h_ps = fromSIunit_h(h2_pT(p, T2_ps(p, s)))
  282.  Case 3
  283.    h_ps = fromSIunit_h(h3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)))
  284.  Case 4
  285.    xs = x4_ps(p, s)
  286.    h_ps = fromSIunit_h(xs * h4V_p(p) + (1 - xs) * h4L_p(p))
  287.  Case 5
  288.    h_ps = fromSIunit_h(h5_pT(p, T5_ps(p, s)))
  289.  Case Else
  290.   h_ps = CVErr(xlErrValue)
  291.  End Select
  292. End Function
  293. Function h_px(ByVal p As Double, ByVal x As Double) As Double
  294.  Dim hL As Double
  295.  Dim hV As Double
  296.  p = toSIunit_p(p)
  297.  x = toSIunit_x(x)
  298.  If x > 1 Or x < 0 Or p >= 22.064 Then
  299.    h_px = CVErr(xlErrValue)
  300.    Exit Function
  301.  End If
  302.  hL = h4L_p(p)
  303.  hV = h4V_p(p)
  304.  h_px = fromSIunit_h(hL + x * (hV - hL))
  305. End Function
  306. Function h_Tx(ByVal t As Double, ByVal x As Double) As Double
  307.  Dim hL As Double
  308.  Dim hV As Double
  309.  Dim p As Double
  310.  t = toSIunit_T(t)
  311.  x = toSIunit_x(x)
  312.  If x > 1 Or x < 0 Or t >= 647.096 Then
  313.    h_Tx = CVErr(xlErrValue)
  314.    Exit Function
  315.  End If
  316.  p = p4_T(t)
  317.  hL = h4L_p(p)
  318.  hV = h4V_p(p)
  319.  h_Tx = fromSIunit_h(hL + x * (hV - hL))
  320. End Function
  321. Function h_prho(ByVal p As Double, ByVal rho As Double) As Double
  322.   Dim hL, hV, vL, vV, x As Double
  323.   p = toSIunit_p(p)
  324.   rho = 1 / toSIunit_v(1 / rho)
  325.   Select Case Region_prho(p, rho)
  326.   Case 1
  327.     h_prho = fromSIunit_h(h1_pT(p, T1_prho(p, rho)))
  328.   Case 2
  329.     h_prho = fromSIunit_h(h2_pT(p, T2_prho(p, rho)))
  330.   Case 3
  331.     h_prho = fromSIunit_h(h3_rhoT(rho, T3_prho(p, rho)))
  332.   Case 4
  333.    If p < 16.529 Then
  334.      vV = v2_pT(p, T4_p(p))
  335.      vL = v1_pT(p, T4_p(p))
  336.    Else
  337.      vV = v3_ph(p, h4V_p(p))
  338.      vL = v3_ph(p, h4L_p(p))
  339.     End If
  340.     hV = h4V_p(p)
  341.     hL = h4L_p(p)
  342.   x = (1 / rho - vL) / (vV - vL)
  343.   h_prho = fromSIunit_h((1 - x) * hL + x * hV)
  344.  Case 5
  345.    h_prho = fromSIunit_h(h5_pT(p, T5_prho(p, rho)))
  346.  Case Else
  347.    h_prho = CVErr(xlErrValue)
  348.  End Select
  349. End Function
  350.  
  351.  
  352. '***********************************************************************************************************
  353. '*1.5 Specific Volume (v)
  354. Function vV_p(ByVal p As Double) As Double
  355.  p = toSIunit_p(p)
  356.  If p > 0.000611657 And p < 22.06395 Then
  357.   If p < 16.529 Then
  358.    vV_p = fromSIunit_v(v2_pT(p, T4_p(p)))
  359.   Else
  360.    vV_p = fromSIunit_v(v3_ph(p, h4V_p(p)))
  361.   End If
  362.  Else
  363.    vV_p = CVErr(xlErrValue)
  364.  End If
  365. End Function
  366. Function vL_p(ByVal p As Double) As Double
  367.  p = toSIunit_p(p)
  368.  If p > 0.000611657 And p < 22.06395 Then
  369.   If p < 16.529 Then
  370.    vL_p = fromSIunit_v(v1_pT(p, T4_p(p)))
  371.   Else
  372.    vL_p = fromSIunit_v(v3_ph(p, h4L_p(p)))
  373.   End If
  374.  Else
  375.    vL_p = CVErr(xlErrValue)
  376.  End If
  377. End Function
  378. Function vV_T(ByVal t As Double) As Double
  379.  t = toSIunit_T(t)
  380.  If t > 273.15 And t < 647.096 Then
  381.   If t <= 623.15 Then
  382.    vV_T = fromSIunit_v(v2_pT(p4_T(t), t))
  383.   Else
  384.    vV_T = fromSIunit_v(v3_ph(p4_T(t), h4V_p(p4_T(t))))
  385.   End If
  386.  Else
  387.    vV_T = CVErr(xlErrValue)
  388.  End If
  389. End Function
  390. Function vL_T(ByVal t As Double) As Double
  391.  t = toSIunit_T(t)
  392.  If t > 273.15 And t < 647.096 Then
  393.   If t <= 623.15 Then
  394.    vL_T = fromSIunit_v(v1_pT(p4_T(t), t))
  395.   Else
  396.    vL_T = fromSIunit_v(v3_ph(p4_T(t), h4L_p(p4_T(t))))
  397.   End If
  398.  Else
  399.    vL_T = CVErr(xlErrValue)
  400.  End If
  401. End Function
  402. Function v_pT(ByVal p As Double, ByVal t As Double) As Double
  403.  p = toSIunit_p(p)
  404.  t = toSIunit_T(t)
  405.  Select Case region_pT(p, t)
  406.  Case 1
  407.    v_pT = fromSIunit_v(v1_pT(p, t))
  408.  Case 2
  409.    v_pT = fromSIunit_v(v2_pT(p, t))
  410.  Case 3
  411.    v_pT = fromSIunit_v(v3_ph(p, h3_pT(p, t)))
  412.  Case 4
  413.    v_pT = CVErr(xlErrValue)
  414.  Case 5
  415.    v_pT = fromSIunit_v(v5_pT(p, t))
  416.  Case Else
  417.   v_pT = CVErr(xlErrValue)
  418.  End Select
  419. End Function
  420.  
  421. Function v_ph(ByVal p As Double, ByVal h As Double) As Double
  422.  Dim xs As Double
  423.  Dim v4V As Double
  424.  Dim v4L As Double
  425.  p = toSIunit_p(p)
  426.  h = toSIunit_h(h)
  427.  Select Case region_ph(p, h)
  428.  Case 1
  429.    v_ph = fromSIunit_v(v1_pT(p, T1_ph(p, h)))
  430.  Case 2
  431.    v_ph = fromSIunit_v(v2_pT(p, T2_ph(p, h)))
  432.  Case 3
  433.    v_ph = fromSIunit_v(v3_ph(p, h))
  434.  Case 4
  435.    xs = x4_ph(p, h)
  436.    If p < 16.529 Then
  437.      v4V = v2_pT(p, T4_p(p))
  438.      v4L = v1_pT(p, T4_p(p))
  439.    Else
  440.      v4V = v3_ph(p, h4V_p(p))
  441.      v4L = v3_ph(p, h4L_p(p))
  442.     End If
  443.     v_ph = fromSIunit_v((xs * v4V + (1 - xs) * v4L))
  444.  Case 5
  445.    v_ph = fromSIunit_v(v5_pT(p, T5_ph(p, h)))
  446.  Case Else
  447.   v_ph = CVErr(xlErrValue)
  448.  End Select
  449. End Function
  450. Function v_ps(ByVal p As Double, ByVal s As Double) As Double
  451.  Dim xs As Double
  452.  Dim v4V As Double
  453.  Dim v4L As Double
  454.  p = toSIunit_p(p)
  455.  s = toSIunit_s(s)
  456.  Select Case region_ps(p, s)
  457.  Case 1
  458.    v_ps = fromSIunit_v(v1_pT(p, T1_ps(p, s)))
  459.  Case 2
  460.    v_ps = fromSIunit_v(v2_pT(p, T2_ps(p, s)))
  461.  Case 3
  462.    v_ps = fromSIunit_v(v3_ps(p, s))
  463.  Case 4
  464.    xs = x4_ps(p, s)
  465.    If p < 16.529 Then
  466.      v4V = v2_pT(p, T4_p(p))
  467.      v4L = v1_pT(p, T4_p(p))
  468.    Else
  469.      v4V = v3_ph(p, h4V_p(p))
  470.      v4L = v3_ph(p, h4L_p(p))
  471.     End If
  472.     v_ps = fromSIunit_v((xs * v4V + (1 - xs) * v4L))
  473.  Case 5
  474.    v_ps = fromSIunit_v(v5_pT(p, T5_ps(p, s)))
  475.  Case Else
  476.    v_ps = CVErr(xlErrValue)
  477.  End Select
  478. End Function
  479.  
  480. '***********************************************************************************************************
  481. '*1.6 Density (rho)
  482. ' Density is calculated as 1/v
  483. Function rhoV_p(ByVal p As Double) As Double
  484.   rhoV_p = 1 / vV_p(p)
  485. End Function
  486. Function rhoL_p(ByVal p As Double) As Double
  487.   rhoL_p = 1 / vL_p(p)
  488. End Function
  489. Function rhoL_T(ByVal t As Double) As Double
  490.   rhoL_T = 1 / vL_T(t)
  491. End Function
  492. Function rhoV_T(ByVal t As Double) As Double
  493.   rhoV_T = 1 / vV_T(t)
  494. End Function
  495. Function rho_pT(ByVal p As Double, ByVal t As Double) As Double
  496.   rho_pT = 1 / v_pT(p, t)
  497. End Function
  498. Function rho_ph(ByVal p As Double, ByVal h As Double) As Double
  499.   rho_ph = 1 / v_ph(p, h)
  500. End Function
  501. Function rho_ps(ByVal p As Double, ByVal s As Double) As Double
  502.   rho_ps = 1 / v_ps(p, s)
  503. End Function
  504.  
  505. '***********************************************************************************************************
  506. '*1.7 Specific entropy (s)
  507. Function sV_p(ByVal p As Double) As Double
  508.  p = toSIunit_p(p)
  509.  If p > 0.000611657 And p < 22.06395 Then
  510.   If p < 16.529 Then
  511.    sV_p = fromSIunit_s(s2_pT(p, T4_p(p)))
  512.   Else
  513.    sV_p = fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)))
  514.   End If
  515.  Else
  516.    sV_p = CVErr(xlErrValue)
  517.  End If
  518. End Function
  519. Function sL_p(ByVal p As Double) As Double
  520.  p = toSIunit_p(p)
  521.  If p > 0.000611657 And p < 22.06395 Then
  522.   If p < 16.529 Then
  523.    sL_p = fromSIunit_s(s1_pT(p, T4_p(p)))
  524.   Else
  525.    sL_p = fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)))
  526.   End If
  527.  Else
  528.    sL_p = CVErr(xlErrValue)
  529.  End If
  530. End Function
  531. Function sV_T(ByVal t As Double) As Double
  532.  t = toSIunit_T(t)
  533.  If t > 273.15 And t < 647.096 Then
  534.   If t <= 623.15 Then
  535.    sV_T = fromSIunit_s(s2_pT(p4_T(t), t))
  536.   Else
  537.    sV_T = fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(t), h4V_p(p4_T(t)))), t))
  538.   End If
  539.  Else
  540.    sV_T = CVErr(xlErrValue)
  541.  End If
  542. End Function
  543. Function sL_T(ByVal t As Double) As Double
  544.  t = toSIunit_T(t)
  545.  If t > 273.15 And t < 647.096 Then
  546.   If t <= 623.15 Then
  547.    sL_T = fromSIunit_s(s1_pT(p4_T(t), t))
  548.   Else
  549.    sL_T = fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(t), h4L_p(p4_T(t)))), t))
  550.   End If
  551.  Else
  552.    sL_T = CVErr(xlErrValue)
  553.  End If
  554. End Function
  555. Function s_pT(ByVal p As Double, ByVal t As Double) As Double
  556.  p = toSIunit_p(p)
  557.  t = toSIunit_T(t)
  558.  Select Case region_pT(p, t)
  559.  Case 1
  560.    s_pT = fromSIunit_s(s1_pT(p, t))
  561.  Case 2
  562.    s_pT = fromSIunit_s(s2_pT(p, t))
  563.  Case 3
  564.    s_pT = fromSIunit_s(s3_rhoT(1 / v3_ph(p, h3_pT(p, t)), t))
  565.  Case 4
  566.    s_pT = CVErr(xlErrValue)
  567.  Case 5
  568.    s_pT = fromSIunit_s(s5_pT(p, t))
  569.  Case Else
  570.   s_pT = CVErr(xlErrValue)
  571.  End Select
  572. End Function
  573. Function s_ph(ByVal p As Double, ByVal h As Double) As Double
  574.  Dim Ts As Double
  575.  Dim xs As Double
  576.  Dim s4V As Double
  577.  Dim s4L As Double
  578.  Dim v4V As Double
  579.  Dim v4L As Double
  580.  p = toSIunit_p(p)
  581.  h = toSIunit_h(h)
  582.  Select Case region_ph(p, h)
  583.  Case 1
  584.    s_ph = fromSIunit_s(s1_pT(p, T1_ph(p, h)))
  585.  Case 2
  586.    s_ph = fromSIunit_s(s2_pT(p, T2_ph(p, h)))
  587.  Case 3
  588.    s_ph = fromSIunit_s(s3_rhoT(1 / v3_ph(p, h), T3_ph(p, h)))
  589.  Case 4
  590.    Ts = T4_p(p)
  591.    xs = x4_ph(p, h)
  592.    If p < 16.529 Then
  593.      s4V = s2_pT(p, Ts)
  594.      s4L = s1_pT(p, Ts)
  595.    Else
  596.      v4V = v3_ph(p, h4V_p(p))
  597.      s4V = s3_rhoT(1 / v4V, Ts)
  598.      v4L = v3_ph(p, h4L_p(p))
  599.      s4L = s3_rhoT(1 / v4L, Ts)
  600.     End If
  601.    s_ph = fromSIunit_s((xs * s4V + (1 - xs) * s4L))
  602.  Case 5
  603.    s_ph = fromSIunit_s(s5_pT(p, T5_ph(p, h)))
  604.  Case Else
  605.   s_ph = CVErr(xlErrValue)
  606.  End Select
  607. End Function
  608. '***********************************************************************************************************
  609. '*1.8 Specific internal energy (u)
  610. Function uV_p(ByVal p As Double) As Double
  611.  p = toSIunit_p(p)
  612.  If p > 0.000611657 And p < 22.06395 Then
  613.   If p < 16.529 Then
  614.    uV_p = fromSIunit_u(u2_pT(p, T4_p(p)))
  615.   Else
  616.    uV_p = fromSIunit_u(u3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)))
  617.   End If
  618.  Else
  619.    uV_p = CVErr(xlErrValue)
  620.  End If
  621. End Function
  622. Function uL_p(ByVal p As Double) As Double
  623.  p = toSIunit_p(p)
  624.  If p > 0.000611657 And p < 22.06395 Then
  625.   If p < 16.529 Then
  626.    uL_p = fromSIunit_u(u1_pT(p, T4_p(p)))
  627.   Else
  628.    uL_p = fromSIunit_u(u3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)))
  629.   End If
  630.  Else
  631.    uL_p = CVErr(xlErrValue)
  632.  End If
  633. End Function
  634. Function uV_T(ByVal t As Double) As Double
  635.  t = toSIunit_T(t)
  636.  If t > 273.15 And t < 647.096 Then
  637.   If t <= 623.15 Then
  638.    uV_T = fromSIunit_u(u2_pT(p4_T(t), t))
  639.   Else
  640.    uV_T = fromSIunit_u(u3_rhoT(1 / (v3_ph(p4_T(t), h4V_p(p4_T(t)))), t))
  641.   End If
  642.  Else
  643.    uV_T = CVErr(xlErrValue)
  644.  End If
  645. End Function
  646. Function uL_T(ByVal t As Double) As Double
  647.  t = toSIunit_T(t)
  648.  If t > 273.15 And t < 647.096 Then
  649.   If t <= 623.15 Then
  650.    uL_T = fromSIunit_u(u1_pT(p4_T(t), t))
  651.   Else
  652.    uL_T = fromSIunit_u(u3_rhoT(1 / (v3_ph(p4_T(t), h4L_p(p4_T(t)))), t))
  653.   End If
  654.  Else
  655.    uL_T = CVErr(xlErrValue)
  656.  End If
  657. End Function
  658. Function u_pT(ByVal p As Double, ByVal t As Double) As Double
  659.  p = toSIunit_p(p)
  660.  t = toSIunit_T(t)
  661.  Select Case region_pT(p, t)
  662.  Case 1
  663.    u_pT = fromSIunit_u(u1_pT(p, t))
  664.  Case 2
  665.    u_pT = fromSIunit_u(u2_pT(p, t))
  666.  Case 3
  667.    u_pT = fromSIunit_u(u3_rhoT(1 / v3_ph(p, h3_pT(p, t)), t))
  668.  Case 4
  669.    u_pT = CVErr(xlErrValue)
  670.  Case 5
  671.    u_pT = fromSIunit_u(u5_pT(p, t))
  672.  Case Else
  673.   u_pT = CVErr(xlErrValue)
  674.  End Select
  675. End Function
  676. Function u_ph(ByVal p As Double, ByVal h As Double) As Double
  677.  Dim Ts As Double
  678.  Dim xs As Double
  679.  Dim u4v As Double
  680.  Dim u4L As Double
  681.  Dim v4V As Double
  682.  Dim v4L As Double
  683.  p = toSIunit_p(p)
  684.  h = toSIunit_h(h)
  685.  Select Case region_ph(p, h)
  686.  Case 1
  687.    u_ph = fromSIunit_u(u1_pT(p, T1_ph(p, h)))
  688.  Case 2
  689.    u_ph = fromSIunit_u(u2_pT(p, T2_ph(p, h)))
  690.  Case 3
  691.    u_ph = fromSIunit_u(u3_rhoT(1 / v3_ph(p, h), T3_ph(p, h)))
  692.  Case 4
  693.    Ts = T4_p(p)
  694.    xs = x4_ph(p, h)
  695.    If p < 16.529 Then
  696.      u4v = u2_pT(p, Ts)
  697.      u4L = u1_pT(p, Ts)
  698.    Else
  699.      v4V = v3_ph(p, h4V_p(p))
  700.      u4v = u3_rhoT(1 / v4V, Ts)
  701.      v4L = v3_ph(p, h4L_p(p))
  702.      u4L = u3_rhoT(1 / v4L, Ts)
  703.    End If
  704.    u_ph = fromSIunit_u((xs * u4v + (1 - xs) * u4L))
  705.  Case 5
  706.    Ts = T5_ph(p, h)
  707.    u_ph = fromSIunit_u(u5_pT(p, Ts))
  708.  Case Else
  709.   u_ph = CVErr(xlErrValue)
  710.  End Select
  711. End Function
  712. Function u_ps(ByVal p As Double, ByVal s As Double) As Double
  713.  Dim x As Double
  714.  Dim u4v, uLp, uVp, u4L As Double
  715.  p = toSIunit_p(p)
  716.  s = toSIunit_s(s)
  717.  Select Case region_ps(p, s)
  718.  Case 1
  719.    u_ps = fromSIunit_u(u1_pT(p, T1_ps(p, s)))
  720.  Case 2
  721.    u_ps = fromSIunit_u(u2_pT(p, T2_ps(p, s)))
  722.  Case 3
  723.    u_ps = fromSIunit_u(u3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)))
  724.  Case 4
  725.    If p < 16.529 Then
  726.      uLp = u1_pT(p, T4_p(p))
  727.      uVp = u2_pT(p, T4_p(p))
  728.    Else
  729.      uLp = u3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p))
  730.      uVp = u3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p))
  731.    End If
  732.    x = x4_ps(p, s)
  733.    u_ps = fromSIunit_u((x * uVp + (1 - x) * uLp))
  734.  Case 5
  735.    u_ps = fromSIunit_u(u5_pT(p, T5_ps(p, s)))
  736.  Case Else
  737.   u_ps = CVErr(xlErrValue)
  738.  End Select
  739. End Function
  740. '***********************************************************************************************************
  741. '*1.9 Specific isobaric heat capacity (Cp)
  742. Function CpV_p(ByVal p As Double) As Double
  743.  p = toSIunit_p(p)
  744.  If p > 0.000611657 And p < 22.06395 Then
  745.   If p < 16.529 Then
  746.    CpV_p = fromSIunit_Cp(Cp2_pT(p, T4_p(p)))
  747.   Else
  748.    CpV_p = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)))
  749.   End If
  750.  Else
  751.    CpV_p = CVErr(xlErrValue)
  752.  End If
  753. End Function
  754. Function CpL_p(ByVal p As Double) As Double
  755.  Dim t, h, v As Double
  756.  p = toSIunit_p(p)
  757.  If p > 0.000611657 And p < 22.06395 Then
  758.   If p < 16.529 Then
  759.    CpL_p = fromSIunit_Cp(Cp1_pT(p, T4_p(p)))
  760.   Else
  761.   t = T4_p(p)
  762.   h = h4L_p(p)
  763.   v = v3_ph(p, h4L_p(p))
  764.  
  765.    CpL_p = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)))
  766.   End If
  767.  Else
  768.    CpL_p = CVErr(xlErrValue)
  769.  End If
  770. End Function
  771. Function CpV_T(ByVal t As Double) As Double
  772.  t = toSIunit_T(t)
  773.  If t > 273.15 And t < 647.096 Then
  774.   If t <= 623.15 Then
  775.    CpV_T = fromSIunit_Cp(Cp2_pT(p4_T(t), t))
  776.   Else
  777.    CpV_T = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p4_T(t), h4V_p(p4_T(t)))), t))
  778.   End If
  779.  Else
  780.    CpV_T = CVErr(xlErrValue)
  781.  End If
  782. End Function
  783. Function CpL_T(ByVal t As Double) As Double
  784.  t = toSIunit_T(t)
  785.  If t > 273.15 And t < 647.096 Then
  786.   If t <= 623.15 Then
  787.    CpL_T = fromSIunit_Cp(Cp1_pT(p4_T(t), t))
  788.   Else
  789.    CpL_T = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p4_T(t), h4L_p(p4_T(t)))), t))
  790.   End If
  791.  Else
  792.    CpL_T = CVErr(xlErrValue)
  793.  End If
  794. End Function
  795. Function Cp_pT(ByVal p As Double, ByVal t As Double) As Double
  796.  p = toSIunit_p(p)
  797.  t = toSIunit_T(t)
  798.  Select Case region_pT(p, t)
  799.  Case 1
  800.    Cp_pT = fromSIunit_Cp(Cp1_pT(p, t))
  801.  Case 2
  802.    Cp_pT = fromSIunit_Cp(Cp2_pT(p, t))
  803.  Case 3
  804.    Cp_pT = fromSIunit_Cp(Cp3_rhoT(1 / v3_ph(p, h3_pT(p, t)), t))
  805.  Case 4
  806.    Cp_pT = CVErr(xlErrValue)
  807.  Case 5
  808.    Cp_pT = fromSIunit_Cp(Cp5_pT(p, t))
  809.  Case Else
  810.   Cp_pT = CVErr(xlErrValue)
  811.  End Select
  812. End Function
  813. Function Cp_ph(ByVal p As Double, ByVal h As Double) As Double
  814.  p = toSIunit_p(p)
  815.  h = toSIunit_h(h)
  816.  Select Case region_ph(p, h)
  817.  Case 1
  818.    Cp_ph = fromSIunit_Cp(Cp1_pT(p, T1_ph(p, h)))
  819.  Case 2
  820.    Cp_ph = fromSIunit_Cp(Cp2_pT(p, T2_ph(p, h)))
  821.  Case 3
  822.    Cp_ph = fromSIunit_Cp(Cp3_rhoT(1 / v3_ph(p, h), T3_ph(p, h)))
  823.  Case 4
  824.    Cp_ph = CVErr(xlErrValue) '#Not def. for mixture"
  825. Case 5
  826.    Cp_ph = fromSIunit_Cp(Cp5_pT(p, T5_ph(p, h)))
  827.  Case Else
  828.   Cp_ph = CVErr(xlErrValue)
  829.  End Select
  830. End Function
  831. Function Cp_ps(ByVal p As Double, ByVal s As Double) As Double
  832.  p = toSIunit_p(p)
  833.  s = toSIunit_s(s)
  834.  Select Case region_ps(p, s)
  835.  Case 1
  836.    Cp_ps = fromSIunit_Cp(Cp1_pT(p, T1_ps(p, s)))
  837.  Case 2
  838.    Cp_ps = fromSIunit_Cp(Cp2_pT(p, T2_ps(p, s)))
  839.  Case 3
  840.    Cp_ps = fromSIunit_Cp(Cp3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)))
  841.  Case 4
  842.    Cp_ps = CVErr(xlErrValue) '#Not def. for mixture"
  843. Case 5
  844.    Cp_ps = fromSIunit_Cp(Cp5_pT(p, T5_ps(p, s)))
  845.  Case Else
  846.   Cp_ps = CVErr(xlErrValue)
  847.  End Select
  848. End Function
  849. '***********************************************************************************************************
  850. '*1.10 Specific isochoric heat capacity (Cv)
  851. Function CvV_p(ByVal p As Double) As Double
  852.  p = toSIunit_p(p)
  853.  If p > 0.000611657 And p < 22.06395 Then
  854.   If p < 16.529 Then
  855.    CvV_p = fromSIunit_Cv(Cv2_pT(p, T4_p(p)))
  856.   Else
  857.    CvV_p = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)))
  858.   End If
  859.  Else
  860.    CvV_p = CVErr(xlErrValue)
  861.  End If
  862. End Function
  863. Function CvL_p(ByVal p As Double) As Double
  864.  p = toSIunit_p(p)
  865.  If p > 0.000611657 And p < 22.06395 Then
  866.   If p < 16.529 Then
  867.    CvL_p = fromSIunit_Cv(Cv1_pT(p, T4_p(p)))
  868.   Else
  869.    CvL_p = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)))
  870.   End If
  871.  Else
  872.    CvL_p = CVErr(xlErrValue)
  873.  End If
  874. End Function
  875. Function CvV_T(ByVal t As Double) As Double
  876.  t = toSIunit_T(t)
  877.  If t > 273.15 And t < 647.096 Then
  878.   If t <= 623.15 Then
  879.    CvV_T = fromSIunit_Cv(Cv2_pT(p4_T(t), t))
  880.   Else
  881.    CvV_T = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p4_T(t), h4V_p(p4_T(t)))), t))
  882.   End If
  883.  Else
  884.    CvV_T = CVErr(xlErrValue)
  885.  End If
  886. End Function
  887. Function CvL_T(ByVal t As Double) As Double
  888.  t = toSIunit_T(t)
  889.  If t > 273.15 And t < 647.096 Then
  890.   If t <= 623.15 Then
  891.    CvL_T = fromSIunit_Cv(Cv1_pT(p4_T(t), t))
  892.   Else
  893.    CvL_T = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p4_T(t), h4L_p(p4_T(t)))), t))
  894.   End If
  895.  Else
  896.    CvL_T = CVErr(xlErrValue)
  897.  End If
  898. End Function
  899. Function Cv_pT(ByVal p As Double, ByVal t As Double) As Double
  900.  p = toSIunit_p(p)
  901.  t = toSIunit_T(t)
  902.  Select Case region_pT(p, t)
  903.  Case 1
  904.    Cv_pT = fromSIunit_Cv(Cv1_pT(p, t))
  905.  Case 2
  906.    Cv_pT = fromSIunit_Cv(Cv2_pT(p, t))
  907.  Case 3
  908.    Cv_pT = fromSIunit_Cv(Cv3_rhoT(1 / v3_ph(p, h3_pT(p, t)), t))
  909.  Case 4
  910.    Cv_pT = CVErr(xlErrValue)
  911.  Case 5
  912.    Cv_pT = fromSIunit_Cv(Cv5_pT(p, t))
  913.  Case Else
  914.   Cv_pT = CVErr(xlErrValue)
  915.  End Select
  916. End Function
  917. Function Cv_ph(ByVal p, ByVal h) As Double
  918.  p = toSIunit_p(p)
  919.  h = toSIunit_h(h)
  920.  Select Case region_ph(p, h)
  921.  Case 1
  922.    Cv_ph = fromSIunit_Cv(Cv1_pT(p, T1_ph(p, h)))
  923.  Case 2
  924.    Cv_ph = fromSIunit_Cv(Cv2_pT(p, T2_ph(p, h)))
  925.  Case 3
  926.    Cv_ph = fromSIunit_Cv(Cv3_rhoT(1 / v3_ph(p, h), T3_ph(p, h)))
  927.  Case 4
  928.    Cv_ph = CVErr(xlErrValue) '#Not def. for mixture"
  929. Case 5
  930.    Cv_ph = fromSIunit_Cv(Cv5_pT(p, T5_ph(p, h)))
  931.  Case Else
  932.   Cv_ph = CVErr(xlErrValue)
  933.  End Select
  934. End Function
  935.  
  936. Function Cv_ps(ByVal p As Double, ByVal s As Double) As Double
  937.  p = toSIunit_p(p)
  938.  s = toSIunit_s(s)
  939.  Select Case region_ps(p, s)
  940.  Case 1
  941.    Cv_ps = fromSIunit_Cv(Cv1_pT(p, T1_ps(p, s)))
  942.  Case 2
  943.    Cv_ps = fromSIunit_Cv(Cv2_pT(p, T2_ps(p, s)))
  944.  Case 3
  945.    Cv_ps = fromSIunit_Cv(Cv3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)))
  946.  Case 4
  947.    Cv_ps = CVErr(xlErrValue) '#Not def. for mixture
  948. Case 5
  949.    Cv_ps = fromSIunit_Cv(Cv5_pT(p, T5_ps(p, s)))
  950.  Case Else
  951.   Cv_ps = CVErr(xlErrValue)
  952.  End Select
  953. End Function
  954.  
  955.  
  956. '***********************************************************************************************************
  957. '*1.11 Speed of sound
  958. Function wV_p(ByVal p As Double) As Double
  959.  p = toSIunit_p(p)
  960.  If p > 0.000611657 And p < 22.06395 Then
  961.   If p < 16.529 Then
  962.    wV_p = fromSIunit_w(w2_pT(p, T4_p(p)))
  963.   Else
  964.    wV_p = fromSIunit_w(w3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)))
  965.   End If
  966.  Else
  967.    wV_p = CVErr(xlErrValue)
  968.  End If
  969. End Function
  970. Function wL_p(ByVal p As Double) As Double
  971.  p = toSIunit_p(p)
  972.  If p > 0.000611657 And p < 22.06395 Then
  973.   If p < 16.529 Then
  974.    wL_p = fromSIunit_w(w1_pT(p, T4_p(p)))
  975.   Else
  976.    wL_p = fromSIunit_w(w3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)))
  977.   End If
  978.  Else
  979.    wL_p = CVErr(xlErrValue)
  980.  End If
  981. End Function
  982. Function wV_T(ByVal t As Double) As Double
  983.  t = toSIunit_T(t)
  984.  If t > 273.15 And t < 647.096 Then
  985.   If t <= 623.15 Then
  986.    wV_T = fromSIunit_w(w2_pT(p4_T(t), t))
  987.   Else
  988.    wV_T = fromSIunit_w(w3_rhoT(1 / (v3_ph(p4_T(t), h4V_p(p4_T(t)))), t))
  989.   End If
  990.  Else
  991.    wV_T = CVErr(xlErrValue)
  992.  End If
  993. End Function
  994. Function wL_T(ByVal t As Double) As Double
  995.  t = toSIunit_T(t)
  996.  If t > 273.15 And t < 647.096 Then
  997.   If t <= 623.15 Then
  998.    wL_T = fromSIunit_w(w1_pT(p4_T(t), t))
  999.   Else
  1000.    wL_T = fromSIunit_w(w3_rhoT(1 / (v3_ph(p4_T(t), h4L_p(p4_T(t)))), t))
  1001.   End If
  1002.  Else
  1003.    wL_T = CVErr(xlErrValue)
  1004.  End If
  1005. End Function
  1006. Function w_pT(ByVal p As Double, ByVal t As Double) As Double
  1007.  p = toSIunit_p(p)
  1008.  t = toSIunit_T(t)
  1009.  Select Case region_pT(p, t)
  1010.  Case 1
  1011.    w_pT = fromSIunit_w(w1_pT(p, t))
  1012.  Case 2
  1013.    w_pT = fromSIunit_w(w2_pT(p, t))
  1014.  Case 3
  1015.    w_pT = fromSIunit_w(w3_rhoT(1 / v3_ph(p, h3_pT(p, t)), t))
  1016.  Case 4
  1017.    w_pT = CVErr(xlErrValue)
  1018.  Case 5
  1019.    w_pT = fromSIunit_w(w5_pT(p, t))
  1020.  Case Else
  1021.   w_pT = CVErr(xlErrValue)
  1022.  End Select
  1023. End Function
  1024.  
  1025. Function w_ph(ByVal p As Double, ByVal h As Double) As Double
  1026.  p = toSIunit_p(p)
  1027.  h = toSIunit_h(h)
  1028.  Select Case region_ph(p, h)
  1029.  Case 1
  1030.    w_ph = fromSIunit_w(w1_pT(p, T1_ph(p, h)))
  1031.  Case 2
  1032.    w_ph = fromSIunit_w(w2_pT(p, T2_ph(p, h)))
  1033.  Case 3
  1034.    w_ph = fromSIunit_w(w3_rhoT(1 / v3_ph(p, h), T3_ph(p, h)))
  1035.  Case 4
  1036.    w_ph = CVErr(xlErrValue) '#Not def. for mixture
  1037. Case 5
  1038.    w_ph = fromSIunit_w(w5_pT(p, T5_ph(p, h)))
  1039.  Case Else
  1040.   w_ph = CVErr(xlErrValue)
  1041.  End Select
  1042. End Function
  1043.  
  1044. Function w_ps(ByVal p As Double, ByVal s As Double) As Double
  1045.  p = toSIunit_p(p)
  1046.  s = toSIunit_s(s)
  1047.  Select Case region_ps(p, s)
  1048.  Case 1
  1049.    w_ps = fromSIunit_w(w1_pT(p, T1_ps(p, s)))
  1050.  Case 2
  1051.    w_ps = fromSIunit_w(w2_pT(p, T2_ps(p, s)))
  1052.  Case 3
  1053.    w_ps = fromSIunit_w(w3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)))
  1054.  Case 4
  1055.    w_ps = CVErr(xlErrValue) '#Not def. for mixture
  1056. Case 5
  1057.    w_ps = fromSIunit_w(w5_pT(p, T5_ps(p, s)))
  1058.  Case Else
  1059.   w_ps = CVErr(xlErrValue)
  1060.  End Select
  1061. End Function
  1062. '***********************************************************************************************************
  1063. '*1.12 Viscosity
  1064. Function my_pT(ByVal p As Double, ByVal t As Double) As Double
  1065.  p = toSIunit_p(p)
  1066.  t = toSIunit_T(t)
  1067.  Select Case region_pT(p, t)
  1068.  Case 4
  1069.    my_pT = CVErr(xlErrValue)
  1070.  Case 1, 2, 3, 5
  1071.    my_pT = fromSIunit_my(my_AllRegions_pT(p, t))
  1072.  Case Else
  1073.    my_pT = CVErr(xlErrValue)
  1074.  End Select
  1075. End Function
  1076. Function my_ph(ByVal p As Double, ByVal h As Double) As Double
  1077.  p = toSIunit_p(p)
  1078.  h = toSIunit_h(h)
  1079.  Select Case region_ph(p, h)
  1080.  Case 1, 2, 3, 5
  1081.    my_ph = fromSIunit_my(my_AllRegions_ph(p, h))
  1082.  Case 4
  1083.   my_ph = CVErr(xlErrValue)
  1084.  Case Else
  1085.   my_ph = CVErr(xlErrValue)
  1086.  End Select
  1087. End Function
  1088. Function my_ps(ByVal p As Double, ByVal s As Double) As Double
  1089.  my_ps = my_ph(p, h_ps(p, s))
  1090. End Function
  1091. '***********************************************************************************************************
  1092. '*1.13 Prandtl
  1093. Function Pr_pT(ByVal p As Double, ByVal t As Double) As Double
  1094.   Dim Cp As Double
  1095.   Dim my As Double
  1096.   Dim tc As Double
  1097.   Cp = toSIunit_Cp(Cp_pT(p, t))
  1098.   my = toSIunit_my(my_pT(p, t))
  1099.   tc = toSIunit_tc(tc_pT(p, t))
  1100.   Pr_pT = Cp * 1000 * my / tc
  1101. End Function
  1102. Function Pr_ph(ByVal p As Double, ByVal h As Double) As Double
  1103.   Dim Cp As Double
  1104.   Dim my As Double
  1105.   Dim tc As Double
  1106.   Cp = toSIunit_Cp(Cp_ph(p, h))
  1107.   my = toSIunit_my(my_ph(p, h))
  1108.   tc = toSIunit_tc(tc_ph(p, h))
  1109.   Pr_ph = Cp * 1000 * my / tc
  1110. End Function
  1111. '***********************************************************************************************************
  1112. '*1.14 Kappa
  1113. Function Kappa_pT(ByVal p As Double, ByVal t As Double) As Double
  1114.   Dim Cp As Double
  1115.   Dim Cv As Double
  1116.   Cp = Cp_pT(p, t)
  1117.   Cv = Cv_pT(p, t)
  1118.   Kappa_pT = Cp / Cv
  1119. End Function
  1120.  
  1121. Function Kappa_ph(ByVal p As Double, ByVal h As Double) As Double
  1122.   Dim Cp As Double
  1123.   Dim Cv As Double
  1124.   Cv = Cv_ph(p, h)
  1125.   Cp = Cp_ph(p, h)
  1126.   Kappa_ph = Cp / Cv
  1127. End Function
  1128. '***********************************************************************************************************
  1129. '*1.15 Surface tension
  1130. Function st_t(ByVal t As Double) As Double
  1131.   t = toSIunit_T(t)
  1132.   st_t = fromSIunit_st(Surface_Tension_T(t))
  1133. End Function
  1134. Function st_p(ByVal p As Double) As Double
  1135.    Dim t As Double
  1136.    t = Tsat_p(p)
  1137.    t = toSIunit_T(t)
  1138.    st_p = fromSIunit_st(Surface_Tension_T(t))
  1139. End Function
  1140. '***********************************************************************************************************
  1141. '*1.16 Thermal conductivity
  1142. Function tcL_p(ByVal p As Double) As Double
  1143.   Dim t As Double
  1144.   Dim v As Double
  1145.   t = Tsat_p(p)
  1146.   v = vL_p(p)
  1147.   p = toSIunit_p(p)
  1148.   t = toSIunit_T(t)
  1149.   v = toSIunit_v(v)
  1150.   tcL_p = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1151. End Function
  1152. Function tcV_p(ByVal p As Double) As Double
  1153.   Dim t As Double
  1154.   Dim v As Double
  1155.   t = Tsat_p(p)
  1156.   v = vV_p(p)
  1157.   p = toSIunit_p(p)
  1158.   t = toSIunit_T(t)
  1159.   v = toSIunit_v(v)
  1160.   tcV_p = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1161. End Function
  1162. Function tcL_T(ByVal t As Double) As Double
  1163.   Dim p, v As Double
  1164.   p = psat_T(t)
  1165.   v = vL_T(t)
  1166.   p = toSIunit_p(p)
  1167.   t = toSIunit_T(t)
  1168.   v = toSIunit_v(v)
  1169.   tcL_T = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1170. End Function
  1171. Function tcV_T(ByVal t As Double) As Double
  1172.   Dim p, v As Double
  1173.   p = psat_T(t)
  1174.   v = vV_T(t)
  1175.   p = toSIunit_p(p)
  1176.   t = toSIunit_T(t)
  1177.   v = toSIunit_v(v)
  1178.   tcV_T = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1179. End Function
  1180. Function tc_pT(ByVal p As Double, ByVal t As Double) As Double
  1181.   Dim v As Double
  1182.   v = v_pT(p, t)
  1183.   p = toSIunit_p(p)
  1184.   t = toSIunit_T(t)
  1185.   v = toSIunit_v(v)
  1186.   tc_pT = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1187. End Function
  1188. Function tc_ph(ByVal p As Double, ByVal h As Double) As Double
  1189.   Dim v As Double
  1190.   Dim t As Double
  1191.   v = v_ph(p, h)
  1192.   t = T_ph(p, h)
  1193.   p = toSIunit_p(p)
  1194.   t = toSIunit_T(t)
  1195.   v = toSIunit_v(v)
  1196.   tc_ph = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1197. End Function
  1198. Function tc_hs(ByVal h As Double, ByVal s As Double) As Double
  1199.   Dim p As Double
  1200.   Dim v As Double
  1201.   Dim t As Double
  1202.   p = p_hs(h, s)
  1203.   v = v_ph(p, h)
  1204.   t = T_ph(p, h)
  1205.   p = toSIunit_p(p)
  1206.   t = toSIunit_T(t)
  1207.   v = toSIunit_v(v)
  1208.   tc_hs = fromSIunit_tc(tc_ptrho(p, t, 1 / v))
  1209. End Function
  1210. '***********************************************************************************************************
  1211. '*1.17 Vapour fraction
  1212. Function x_ph(ByVal p As Double, ByVal h As Double) As Double
  1213.  p = toSIunit_p(p)
  1214.  h = toSIunit_h(h)
  1215.   If p > 0.000611657 And p < 22.06395 Then
  1216.     x_ph = fromSIunit_x(x4_ph(p, h))
  1217.   Else
  1218.     x_ph = CVErr(xlErrValue)
  1219.   End If
  1220. End Function
  1221. Function x_ps(ByVal p As Double, ByVal s As Double) As Double
  1222.  p = toSIunit_p(p)
  1223.  s = toSIunit_s(s)
  1224.   If p > 0.000611657 And p < 22.06395 Then
  1225.     x_ps = fromSIunit_x(x4_ps(p, s))
  1226.   Else
  1227.     x_ps = CVErr(xlErrValue)
  1228.   End If
  1229. End Function
  1230. '***********************************************************************************************************
  1231. '*1.18 Vapour Volume Fraction
  1232. Function vx_ph(ByVal p As Double, ByVal h As Double) As Double
  1233.  Dim vL As Double
  1234.  Dim vV As Double
  1235.  Dim xs As Double
  1236.  p = toSIunit_p(p)
  1237.  h = toSIunit_h(h)
  1238.  If p > 0.000611657 And p < 22.06395 Then
  1239.     If p < 16.529 Then
  1240.       vL = v1_pT(p, T4_p(p))
  1241.       vV = v2_pT(p, T4_p(p))
  1242.     Else
  1243.       vL = v3_ph(p, h4L_p(p))
  1244.       vV = v3_ph(p, h4V_p(p))
  1245.     End If
  1246.     xs = x4_ph(p, h)
  1247.     vx_ph = fromSIunit_vx((xs * vV / (xs * vV + (1 - xs) * vL)))
  1248.   Else
  1249.     vx_ph = CVErr(xlErrValue)
  1250.   End If
  1251. End Function
  1252. Function vx_ps(ByVal p As Double, ByVal s As Double) As Double
  1253.  Dim vL As Double
  1254.  Dim vV As Double
  1255.  Dim xs As Double
  1256.  p = toSIunit_p(p)
  1257.  s = toSIunit_s(s)
  1258.  If p > 0.000611657 And p < 22.06395 Then
  1259.     If p < 16.529 Then
  1260.       vL = v1_pT(p, T4_p(p))
  1261.       vV = v2_pT(p, T4_p(p))
  1262.     Else
  1263.       vL = v3_ph(p, h4L_p(p))
  1264.       vV = v3_ph(p, h4V_p(p))
  1265.     End If
  1266.     xs = x4_ps(p, s)
  1267.     vx_ps = fromSIunit_vx((xs * vV / (xs * vV + (1 - xs) * vL)))
  1268.   Else
  1269.     vx_ps = CVErr(xlErrValue)
  1270.   End If
  1271. End Function
  1272.  
  1273. '***********************************************************************************************************
  1274. '*2 IAPWS IF 97 Calling functions                                                                          *
  1275. '***********************************************************************************************************
  1276. '
  1277. '***********************************************************************************************************
  1278. '*2.1 Functions for region 1
  1279. Function v1_pT(ByVal p As Double, ByVal t As Double) As Double
  1280. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1281. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1282. 'Eqution 7, Table 3, Page 6
  1283.  Dim i As Integer
  1284.   Dim ps, tau, g_p As Double
  1285.   Dim I1, J1, n1 As Variant
  1286.   Const R As Double = 0.461526 'kJ/(kg K)
  1287.  I1 = Array(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)
  1288.   J1 = Array(-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)
  1289.   n1 = Array(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)
  1290.   ps = p / 16.53
  1291.   tau = 1386 / t
  1292.   g_p = 0#
  1293.   For i = 0 To 33
  1294.    g_p = g_p - n1(i) * I1(i) * (7.1 - ps) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i)
  1295.   Next i
  1296.  v1_pT = R * t / p * ps * g_p / 1000
  1297. End Function
  1298. Function h1_pT(ByVal p As Double, ByVal t As Double) As Double
  1299. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1300. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1301. 'Eqution 7, Table 3, Page 6
  1302.  Dim i As Integer
  1303.   Dim ps, tau, g_t As Double
  1304.   Dim I1, J1, n1 As Variant
  1305.   Const R As Double = 0.461526 'kJ/(kg K)
  1306.  I1 = Array(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)
  1307.   J1 = Array(-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)
  1308.   n1 = Array(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)
  1309.   p = p / 16.53
  1310.   tau = 1386 / t
  1311.   g_t = 0#
  1312.   For i = 0 To 33
  1313.    g_t = g_t + (n1(i) * (7.1 - p) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1))
  1314.   Next i
  1315.  h1_pT = R * t * tau * g_t
  1316. End Function
  1317. Function u1_pT(ByVal p As Double, ByVal t As Double) As Double
  1318. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1319. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1320. 'Eqution 7, Table 3, Page 6
  1321.  Dim i As Integer
  1322.   Dim tau, g_t, g_p As Double
  1323.   Dim I1, J1, n1 As Variant
  1324.   Const R As Double = 0.461526 'kJ/(kg K)
  1325.  I1 = Array(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)
  1326.   J1 = Array(-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)
  1327.   n1 = Array(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)
  1328.   p = p / 16.53
  1329.   tau = 1386 / t
  1330.   g_t = 0#
  1331.   g_p = 0#
  1332.   For i = 0 To 33
  1333.    g_p = g_p - n1(i) * I1(i) * (7.1 - p) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i)
  1334.    g_t = g_t + (n1(i) * (7.1 - p) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1))
  1335.   Next i
  1336.   u1_pT = R * t * (tau * g_t - p * g_p)
  1337. End Function
  1338. Function s1_pT(ByVal p As Double, ByVal t As Double) As Double
  1339. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1340. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1341. 'Eqution 7, Table 3, Page 6
  1342.  Dim i As Integer
  1343.   Dim g, g_t As Double
  1344.   Dim I1, J1, n1 As Variant
  1345.   Const R As Double = 0.461526 'kJ/(kg K)
  1346.  I1 = Array(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)
  1347.   J1 = Array(-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)
  1348.   n1 = Array(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)
  1349.   p = p / 16.53
  1350.   t = 1386 / t
  1351.   g = 0#
  1352.   g_t = 0#
  1353.   For i = 0 To 33
  1354.    g_t = g_t + (n1(i) * (7.1 - p) ^ I1(i) * J1(i) * (t - 1.222) ^ (J1(i) - 1))
  1355.    g = g + n1(i) * (7.1 - p) ^ I1(i) * (t - 1.222) ^ J1(i)
  1356.   Next i
  1357.   s1_pT = R * t * g_t - R * g
  1358. End Function
  1359. Function Cp1_pT(ByVal p As Double, ByVal t As Double) As Double
  1360. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1361. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1362. 'Eqution 7, Table 3, Page 6
  1363.  Dim i As Integer
  1364.   Dim G_tt As Double
  1365.   Dim I1, J1, n1 As Variant
  1366.   Const R As Double = 0.461526 'kJ/(kg K)
  1367.  I1 = Array(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)
  1368.   J1 = Array(-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)
  1369.   n1 = Array(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)
  1370.   p = p / 16.53
  1371.   t = 1386 / t
  1372.   G_tt = 0#
  1373.   For i = 0 To 33
  1374.    G_tt = G_tt + (n1(i) * (7.1 - p) ^ I1(i) * J1(i) * (J1(i) - 1) * (t - 1.222) ^ (J1(i) - 2))
  1375.   Next i
  1376.   Cp1_pT = -R * t ^ 2 * G_tt
  1377. End Function
  1378. Function Cv1_pT(ByVal p As Double, ByVal t As Double) As Double
  1379. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1380. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1381. 'Eqution 7, Table 3, Page 6
  1382.  Dim i As Integer
  1383.   Dim g_p, g_pp, g_pt, G_tt As Double
  1384.   Dim I1, J1, n1 As Variant
  1385.   Const R As Double = 0.461526 'kJ/(kg K)
  1386.  I1 = Array(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)
  1387.   J1 = Array(-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)
  1388.   n1 = Array(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)
  1389.   p = p / 16.53
  1390.   t = 1386 / t
  1391.   g_p = 0#
  1392.   g_pp = 0#
  1393.   g_pt = 0#
  1394.   G_tt = 0#
  1395.   For i = 0 To 33
  1396.    g_p = g_p - n1(i) * I1(i) * (7.1 - p) ^ (I1(i) - 1) * (t - 1.222) ^ J1(i)
  1397.    g_pp = g_pp + n1(i) * I1(i) * (I1(i) - 1) * (7.1 - p) ^ (I1(i) - 2) * (t - 1.222) ^ J1(i)
  1398.    g_pt = g_pt - n1(i) * I1(i) * (7.1 - p) ^ (I1(i) - 1) * J1(i) * (t - 1.222) ^ (J1(i) - 1)
  1399.    G_tt = G_tt + n1(i) * (7.1 - p) ^ I1(i) * J1(i) * (J1(i) - 1) * (t - 1.222) ^ (J1(i) - 2)
  1400.   Next i
  1401.   Cv1_pT = R * (-(t ^ 2 * G_tt) + (g_p - t * g_pt) ^ 2 / g_pp)
  1402. End Function
  1403. Function w1_pT(ByVal p As Double, ByVal t As Double) As Double
  1404. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1405. '5 Equations for Region 1, Section. 5.1 Basic Equation
  1406. 'Eqution 7, Table 3, Page 6
  1407.  Dim i As Integer
  1408.   Dim g_p, g_pp, g_pt, G_tt, tau As Double
  1409.   Dim I1, J1, n1 As Variant
  1410.   Const R As Double = 0.461526 'kJ/(kg K)
  1411.  I1 = Array(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)
  1412.   J1 = Array(-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)
  1413.   n1 = Array(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)
  1414.   p = p / 16.53
  1415.   tau = 1386 / t
  1416.   g_p = 0#
  1417.   g_pp = 0#
  1418.   g_pt = 0#
  1419.   G_tt = 0#
  1420.   For i = 0 To 33
  1421.    g_p = g_p - n1(i) * I1(i) * (7.1 - p) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i)
  1422.    g_pp = g_pp + n1(i) * I1(i) * (I1(i) - 1) * (7.1 - p) ^ (I1(i) - 2) * (tau - 1.222) ^ J1(i)
  1423.    g_pt = g_pt - n1(i) * I1(i) * (7.1 - p) ^ (I1(i) - 1) * J1(i) * (tau - 1.222) ^ (J1(i) - 1)
  1424.    G_tt = G_tt + n1(i) * (7.1 - p) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2)
  1425.   Next i
  1426.   w1_pT = (1000 * R * t * g_p ^ 2 / ((g_p - tau * g_pt) ^ 2 / (tau ^ 2 * G_tt) - g_pp)) ^ 0.5
  1427. End Function
  1428. Function T1_ph(ByVal p As Double, ByVal h As Double) As Double
  1429. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1430. '5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.1 The Backward Equation T ( p,h )
  1431. 'Eqution 11, Table 6, Page 10
  1432.  Dim i As Integer
  1433.   Dim t As Double
  1434.   Dim I1, J1, n1 As Variant
  1435.   I1 = Array(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6)
  1436.   J1 = Array(0, 1, 2, 6, 22, 32, 0, 1, 2, 3, 4, 10, 32, 10, 32, 10, 32, 32, 32, 32)
  1437.   n1 = Array(-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)
  1438.   h = h / 2500
  1439.   t = 0#
  1440.   For i = 0 To 19
  1441.    t = t + n1(i) * p ^ I1(i) * (h + 1) ^ J1(i)
  1442.   Next i
  1443.   T1_ph = t
  1444. End Function
  1445. Function T1_ps(ByVal p As Double, ByVal s As Double) As Double
  1446. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1447. '5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.2 The Backward Equation T ( p, s )
  1448. 'Eqution 13, Table 8, Page 11
  1449.  Dim i As Integer
  1450.   Dim I1, J1, n1 As Variant
  1451.   I1 = Array(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4)
  1452.   J1 = Array(0, 1, 2, 3, 11, 31, 0, 1, 2, 3, 12, 31, 0, 1, 2, 9, 31, 10, 32, 32)
  1453.   n1 = Array(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)
  1454.   T1_ps = 0#
  1455.   For i = 0 To 19
  1456.    T1_ps = T1_ps + n1(i) * p ^ I1(i) * (s + 2) ^ J1(i)
  1457.   Next i
  1458. End Function
  1459. Function p1_hs(ByVal h As Double, ByVal s As Double) As Double
  1460. '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
  1461. '5 Backward Equation p(h,s) for Region 1
  1462. 'Eqution 1, Table 2, Page 5
  1463.  Dim i As Integer
  1464.   Dim p As Double
  1465.   Dim I1, J1, n1 As Variant
  1466.   I1 = Array(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 5)
  1467.   J1 = Array(0, 1, 2, 4, 5, 6, 8, 14, 0, 1, 4, 6, 0, 1, 10, 4, 1, 4, 0)
  1468.   n1 = Array(-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)
  1469.   h = h / 3400
  1470.   s = s / 7.6
  1471.   p = 0#
  1472.   For i = 0 To 18
  1473.    p = p + n1(i) * (h + 0.05) ^ I1(i) * (s + 0.05) ^ J1(i)
  1474.   Next i
  1475.   p1_hs = p * 100
  1476. End Function
  1477. Function T1_prho(ByVal p As Double, ByVal rho As Double) As Double
  1478.   'Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  1479.  'Solve with half interval method
  1480.  Dim i As Integer
  1481.   Dim Ts, Low_Bound, High_Bound, rhos As Double
  1482.   Low_Bound = 273.15
  1483.   High_Bound = T4_p(p)
  1484.   Do While Abs(rho - rhos) > 0.00001
  1485.     Ts = (Low_Bound + High_Bound) / 2
  1486.     rhos = 1 / v1_pT(p, Ts)
  1487.     If rhos < rho Then
  1488.       High_Bound = Ts
  1489.     Else
  1490.       Low_Bound = Ts
  1491.     End If
  1492.     Loop
  1493.     T1_prho = Ts
  1494. End Function
  1495. '***********************************************************************************************************
  1496. '*2.2 Functions for region 2
  1497.  
  1498. Function v2_pT(ByVal p As Double, ByVal t As Double) As Double
  1499. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1500. '6 Equations for Region 2, Section. 6.1 Basic Equation
  1501. 'Table 11 and 12, Page 14 and 15
  1502.  Dim i As Integer
  1503.   Dim tau, g0_pi, gr_pi As Double
  1504.   Dim Ir, Jr, nr, J0, n0 As Variant
  1505.   Const R As Double = 0.461526 'kJ/(kg K)
  1506.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1507.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1508.   Ir = Array(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)
  1509.   Jr = Array(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)
  1510.   nr = Array(-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)
  1511.   tau = 540 / t
  1512.   g0_pi = 1 / p
  1513.   gr_pi = 0#
  1514.   For i = 0 To 42
  1515.    gr_pi = gr_pi + nr(i) * Ir(i) * p ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i)
  1516.   Next i
  1517.   v2_pT = R * t / p * p * (g0_pi + gr_pi) / 1000
  1518. End Function
  1519. Function h2_pT(ByVal p As Double, ByVal t As Double) As Double
  1520. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1521. '6 Equations for Region 2, Section. 6.1 Basic Equation
  1522. 'Table 11 and 12, Page 14 and 15
  1523.  Dim i As Integer
  1524.   Dim tau, g0_tau, gr_tau As Double
  1525.   Dim Ir, Jr, nr, J0, n0 As Variant
  1526.   Const R As Double = 0.461526 'kJ/(kg K)
  1527.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1528.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1529.   Ir = Array(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)
  1530.   Jr = Array(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)
  1531.   nr = Array(-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)
  1532.   tau = 540 / t
  1533.   g0_tau = 0#
  1534.   For i = 0 To 8
  1535.     g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1)
  1536.   Next i
  1537.   gr_tau = 0#
  1538.   For i = 0 To 42
  1539.    gr_tau = gr_tau + nr(i) * p ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1)
  1540.   Next i
  1541.   h2_pT = R * t * tau * (g0_tau + gr_tau)
  1542. End Function
  1543. Function u2_pT(ByVal p As Double, ByVal t As Double) As Double
  1544. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1545. '6 Equations for Region 2, Section. 6.1 Basic Equation
  1546. 'Table 11 and 12, Page 14 and 15
  1547.  Dim i As Integer
  1548.   Dim tau, g0_tau, g0_pi, gr_pi, gr_tau As Double
  1549.   Dim Ir, Jr, nr, J0, n0 As Variant
  1550.   Const R As Double = 0.461526 'kJ/(kg K)
  1551.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1552.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1553.   Ir = Array(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)
  1554.   Jr = Array(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)
  1555.   nr = Array(-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)
  1556.   tau = 540 / t
  1557.   g0_pi = 1 / p
  1558.   g0_tau = 0#
  1559.   For i = 0 To 8
  1560.     g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1)
  1561.   Next i
  1562.   gr_pi = 0#
  1563.   gr_tau = 0#
  1564.   For i = 0 To 42
  1565.    gr_pi = gr_pi + nr(i) * Ir(i) * p ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i)
  1566.    gr_tau = gr_tau + nr(i) * p ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1)
  1567.   Next i
  1568.   u2_pT = R * t * (tau * (g0_tau + gr_tau) - p * (g0_pi + gr_pi))
  1569. End Function
  1570.  
  1571. Function s2_pT(ByVal p As Double, ByVal t As Double) As Double
  1572. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1573. '6 Equations for Region 2, Section. 6.1 Basic Equation
  1574. 'Table 11 and 12, Page 14 and 15
  1575.  Dim i As Integer
  1576.   Dim tau, g0, g0_tau, gr, gr_tau As Double
  1577.   Dim Ir, Jr, nr, J0, n0 As Variant
  1578.   Const R As Double = 0.461526 'kJ/(kg K)
  1579.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1580.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1581.   Ir = Array(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)
  1582.   Jr = Array(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)
  1583.   nr = Array(-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)
  1584.   tau = 540 / t
  1585.   g0 = Log(p)
  1586.   g0_tau = 0#
  1587.   For i = 0 To 8
  1588.     g0 = g0 + n0(i) * tau ^ J0(i)
  1589.     g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1)
  1590.   Next i
  1591.   gr = 0#
  1592.   gr_tau = 0#
  1593.   For i = 0 To 42
  1594.    gr = gr + nr(i) * p ^ Ir(i) * (tau - 0.5) ^ Jr(i)
  1595.    gr_tau = gr_tau + nr(i) * p ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1)
  1596.   Next i
  1597.   s2_pT = R * (tau * (g0_tau + gr_tau) - (g0 + gr))
  1598. End Function
  1599. Function Cp2_pT(ByVal p As Double, ByVal t As Double) As Double
  1600. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1601. '6 Equations for Region 2, Section. 6.1 Basic Equation
  1602. 'Table 11 and 12, Page 14 and 15
  1603.  Dim i As Integer
  1604.   Dim tau, g0_tautau, gr_tautau As Double
  1605.   Dim Ir, Jr, nr, J0, n0 As Variant
  1606.   Const R As Double = 0.461526 'kJ/(kg K)
  1607.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1608.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1609.   Ir = Array(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)
  1610.   Jr = Array(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)
  1611.   nr = Array(-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)
  1612.   tau = 540 / t
  1613.   g0_tautau = 0#
  1614.   For i = 0 To 8
  1615.     g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2)
  1616.   Next i
  1617.   gr_tautau = 0#
  1618.   For i = 0 To 42
  1619.    gr_tautau = gr_tautau + nr(i) * p ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2)
  1620.   Next i
  1621.   Cp2_pT = -R * tau ^ 2 * (g0_tautau + gr_tautau)
  1622. End Function
  1623. Function Cv2_pT(ByVal p As Double, ByVal t As Double) As Double
  1624.   Dim i As Integer
  1625.   Dim tau, g0_tautau, gr_pi, gr_pitau, gr_pipi, gr_tautau As Double
  1626.   Dim Ir, Jr, nr, J0, n0 As Variant
  1627.   Const R As Double = 0.461526 'kJ/(kg K)
  1628.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1629.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1630.   Ir = Array(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)
  1631.   Jr = Array(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)
  1632.   nr = Array(-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)
  1633.   tau = 540 / t
  1634.   g0_tautau = 0#
  1635.   For i = 0 To 8
  1636.     g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2)
  1637.   Next i
  1638.   gr_pi = 0#
  1639.   gr_pitau = 0#
  1640.   gr_pipi = 0#
  1641.   gr_tautau = 0#
  1642.   For i = 0 To 42
  1643.    gr_pi = gr_pi + nr(i) * Ir(i) * p ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i)
  1644.    gr_pipi = gr_pipi + nr(i) * Ir(i) * (Ir(i) - 1) * p ^ (Ir(i) - 2) * (tau - 0.5) ^ Jr(i)
  1645.    gr_pitau = gr_pitau + nr(i) * Ir(i) * p ^ (Ir(i) - 1) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1)
  1646.    gr_tautau = gr_tautau + nr(i) * p ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2)
  1647.   Next i
  1648.   Cv2_pT = R * (-(tau ^ 2 * (g0_tautau + gr_tautau)) - ((1 + p * gr_pi - tau * p * gr_pitau) ^ 2) / (1 - p ^ 2 * gr_pipi))
  1649. End Function
  1650. Function w2_pT(ByVal p As Double, ByVal t As Double) As Double
  1651. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1652. '6 Equations for Region 2, Section. 6.1 Basic Equation
  1653. 'Table 11 and 12, Page 14 and 15
  1654.  Dim i As Integer
  1655.   Dim tau, g0_tautau, gr_pi, gr_pitau, gr_pipi, gr_tautau As Double
  1656.   Dim Ir, Jr, nr, J0, n0 As Variant
  1657.   Const R As Double = 0.461526 'kJ/(kg K)
  1658.  J0 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)
  1659.   n0 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)
  1660.   Ir = Array(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)
  1661.   Jr = Array(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)
  1662.   nr = Array(-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)
  1663.   tau = 540 / t
  1664.   g0_tautau = 0#
  1665.   For i = 0 To 8
  1666.     g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2)
  1667.   Next i
  1668.   gr_pi = 0#
  1669.   gr_pitau = 0#
  1670.   gr_pipi = 0#
  1671.   gr_tautau = 0#
  1672.   For i = 0 To 42
  1673.    gr_pi = gr_pi + nr(i) * Ir(i) * p ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i)
  1674.    gr_pipi = gr_pipi + nr(i) * Ir(i) * (Ir(i) - 1) * p ^ (Ir(i) - 2) * (tau - 0.5) ^ Jr(i)
  1675.    gr_pitau = gr_pitau + nr(i) * Ir(i) * p ^ (Ir(i) - 1) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1)
  1676.    gr_tautau = gr_tautau + nr(i) * p ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2)
  1677.   Next i
  1678.   w2_pT = (1000 * R * t * (1 + 2 * p * gr_pi + p ^ 2 * gr_pi ^ 2) / ((1 - p ^ 2 * gr_pipi) + (1 + p * gr_pi - tau * p * gr_pitau) ^ 2 / (tau ^ 2 * (g0_tautau + gr_tautau)))) ^ 0.5
  1679. End Function
  1680. Function T2_ph(ByVal p As Double, ByVal h As Double) As Double
  1681.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1682.  '6 Equations for Region 2,6.3.1 The Backward Equations T( p, h ) for Subregions 2a, 2b, and 2c
  1683.  Dim sub_reg As Integer
  1684.   Dim i As Integer
  1685.   Dim Ts, hs As Double
  1686.   Dim Ji, Ii, ni As Variant
  1687.  
  1688.   If p < 4 Then
  1689.    sub_reg = 1
  1690.   Else
  1691.    If p < (905.84278514723 - 0.67955786399241 * h + 1.2809002730136E-04 * h ^ 2) Then
  1692.      sub_reg = 2
  1693.    Else
  1694.      sub_reg = 3
  1695.    End If
  1696.   End If
  1697.  
  1698.   Select Case sub_reg
  1699.   Case 1
  1700.     'Subregion A
  1701.    'Table 20, Eq 22, page 22
  1702.    Ji = Array(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)
  1703.     Ii = Array(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)
  1704.     ni = Array(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)
  1705.     Ts = 0
  1706.     hs = h / 2000
  1707.     For i = 0 To 33
  1708.       Ts = Ts + ni(i) * p ^ (Ii(i)) * (hs - 2.1) ^ Ji(i)
  1709.     Next i
  1710.     T2_ph = Ts
  1711.   Case 2
  1712.     'Subregion B
  1713.    'Table 21, Eq 23, page 23
  1714.    Ji = Array(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)
  1715.     Ii = Array(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)
  1716.     ni = Array(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)
  1717.     Ts = 0
  1718.     hs = h / 2000
  1719.     For i = 0 To 37
  1720.       Ts = Ts + ni(i) * (p - 2) ^ (Ii(i)) * (hs - 2.6) ^ Ji(i)
  1721.     Next i
  1722.     T2_ph = Ts
  1723.   Case Else
  1724.     'Subregion C
  1725.    'Table 22, Eq 24, page 24
  1726.    Ji = Array(0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16, 20, 22)
  1727.     Ii = Array(-7, -7, -6, -6, -5, -5, -2, -2, -1, -1, 0, 0, 1, 1, 2, 6, 6, 6, 6, 6, 6, 6, 6)
  1728.     ni = Array(-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)
  1729.     Ts = 0
  1730.     hs = h / 2000
  1731.     For i = 0 To 22
  1732.       Ts = Ts + ni(i) * (p + 25) ^ (Ii(i)) * (hs - 1.8) ^ Ji(i)
  1733.     Next i
  1734.     T2_ph = Ts
  1735.   End Select
  1736. End Function
  1737. Function T2_ps(ByVal p As Double, ByVal s As Double) As Double
  1738.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1739.  '6 Equations for Region 2,6.3.2 The Backward Equations T( p, s ) for Subregions 2a, 2b, and 2c
  1740.  'Page 26
  1741.  Dim sub_reg As Integer
  1742.   Dim i As Integer
  1743.   Dim teta, sigma As Double
  1744.   Dim Ji, Ii, ni As Variant
  1745.  
  1746.   If p < 4 Then
  1747.      sub_reg = 1
  1748.   Else
  1749.      If s < 5.85 Then
  1750.        sub_reg = 3
  1751.      Else
  1752.        sub_reg = 2
  1753.      End If
  1754.   End If
  1755.   Select Case sub_reg
  1756.   Case 1
  1757.    'Subregion A
  1758.   'Table 25, Eq 25, page 26
  1759.   Ii = Array(-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)
  1760.    Ji = Array(-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)
  1761.    ni = Array(-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)
  1762.    sigma = s / 2
  1763.    teta = 0
  1764.    For i = 0 To 45
  1765.      teta = teta + ni(i) * p ^ Ii(i) * (sigma - 2) ^ Ji(i)
  1766.    Next i
  1767.    T2_ps = teta
  1768.   Case 2
  1769.     'Subregion B
  1770.    'Table 26, Eq 26, page 27
  1771.    Ii = Array(-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)
  1772.     Ji = Array(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)
  1773.     ni = Array(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)
  1774.     sigma = s / 0.7853
  1775.     teta = 0
  1776.     For i = 0 To 43
  1777.       teta = teta + ni(i) * p ^ Ii(i) * (10 - sigma) ^ Ji(i)
  1778.     Next i
  1779.     T2_ps = teta
  1780.   Case Else
  1781.     'Subregion C
  1782.    'Table 27, Eq 27, page 28
  1783.    Ii = Array(-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)
  1784.     Ji = Array(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)
  1785.     ni = Array(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)
  1786.     sigma = s / 2.9251
  1787.     teta = 0
  1788.     For i = 0 To 29
  1789.       teta = teta + ni(i) * p ^ Ii(i) * (2 - sigma) ^ Ji(i)
  1790.     Next i
  1791.     T2_ps = teta
  1792.     End Select
  1793. End Function
  1794. Function p2_hs(ByVal h As Double, ByVal s As Double) As Double
  1795. '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
  1796. 'Chapter 6:Backward Equations p(h,s) for Region 2
  1797.  Dim sub_reg As Integer
  1798.   Dim i As Integer
  1799.   Dim eta, teta, sigma, p As Double
  1800.   Dim Ji, Ii, ni As Variant
  1801.  
  1802.   If h < (-3498.98083432139 + 2575.60716905876 * s - 421.073558227969 * s ^ 2 + 27.6349063799944 * s ^ 3) Then
  1803.     sub_reg = 1
  1804.   Else
  1805.    If s < 5.85 Then
  1806.      sub_reg = 3
  1807.    Else
  1808.      sub_reg = 2
  1809.    End If
  1810.   End If
  1811.   Select Case sub_reg
  1812.   Case 1
  1813.     'Subregion A
  1814.    'Table 6, Eq 3, page 8
  1815.    Ii = Array(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)
  1816.     Ji = Array(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)
  1817.     ni = Array(-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)
  1818.     eta = h / 4200
  1819.     sigma = s / 12
  1820.     p = 0
  1821.     For i = 0 To 28
  1822.       p = p + ni(i) * (eta - 0.5) ^ Ii(i) * (sigma - 1.2) ^ Ji(i)
  1823.     Next i
  1824.     p2_hs = p ^ 4 * 4
  1825.   Case 2
  1826.     'Subregion B
  1827.    'Table 7, Eq 4, page 9
  1828.    Ii = Array(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)
  1829.     Ji = Array(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)
  1830.     ni = Array(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#)
  1831.     eta = h / 4100
  1832.     sigma = s / 7.9
  1833.     For i = 0 To 32
  1834.       p = p + ni(i) * (eta - 0.6) ^ Ii(i) * (sigma - 1.01) ^ Ji(i)
  1835.     Next i
  1836.     p2_hs = p ^ 4 * 100
  1837.   Case Else
  1838.     'Subregion C
  1839.    'Table 8, Eq 5, page 10
  1840.    Ii = Array(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)
  1841.     Ji = Array(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)
  1842.     ni = Array(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)
  1843.     eta = h / 3500
  1844.     sigma = s / 5.9
  1845.     For i = 0 To 30
  1846.       p = p + ni(i) * (eta - 0.7) ^ Ii(i) * (sigma - 1.1) ^ Ji(i)
  1847.     Next i
  1848.     p2_hs = p ^ 4 * 100
  1849.   End Select
  1850. End Function
  1851. Function T2_prho(ByVal p As Double, ByVal rho As Double) As Double
  1852.   'Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  1853.  'Solve with half interval method
  1854.  Dim Low_Bound, High_Bound, rhos, Ts As Double
  1855.  
  1856.   If p < 16.5292 Then
  1857.     Low_Bound = T4_p(p)
  1858.   Else
  1859.     Low_Bound = B23T_p(p)
  1860.   End If
  1861.   High_Bound = 1073.15
  1862.   Do While Abs(rho - rhos) > 0.000001
  1863.     Ts = (Low_Bound + High_Bound) / 2
  1864.     rhos = 1 / v2_pT(p, Ts)
  1865.     If rhos < rho Then
  1866.       High_Bound = Ts
  1867.     Else
  1868.       Low_Bound = Ts
  1869.     End If
  1870.     Loop
  1871.     T2_prho = Ts
  1872. End Function
  1873. '***********************************************************************************************************
  1874. '*2.3 Functions for region 3
  1875.  
  1876. Function p3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1877.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1878.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1879.  'Table 30 and 31, Page 30 and 31
  1880.  Dim i As Integer, Ji, Ii, ni As Variant, delta, tau, fidelta As Double
  1881.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  1882.   Ii = Array(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)
  1883.   Ji = Array(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)
  1884.   ni = Array(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)
  1885.   delta = rho / rhoc
  1886.   tau = tc / t
  1887.   fidelta = 0
  1888.   For i = 1 To 39
  1889.     fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i)
  1890.   Next i
  1891.   fidelta = fidelta + ni(0) / delta
  1892.   p3_rhoT = rho * R * t * delta * fidelta / 1000
  1893. End Function
  1894. Function u3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1895.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1896.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1897.  'Table 30 and 31, Page 30 and 31
  1898.  Dim i As Integer, Ji, Ii, ni As Variant, delta, tau, fitau As Double
  1899.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  1900.   Ii = Array(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)
  1901.   Ji = Array(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)
  1902.   ni = Array(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)
  1903.   delta = rho / rhoc
  1904.   tau = tc / t
  1905.   fitau = 0
  1906.   For i = 1 To 39
  1907.     fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1)
  1908.   Next i
  1909.   u3_rhoT = R * t * (tau * fitau)
  1910. End Function
  1911.  
  1912. Function h3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1913.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1914.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1915.  'Table 30 and 31, Page 30 and 31
  1916.  Dim i As Integer, Ji, Ii, ni As Variant, delta, tau, fidelta, fitau As Double
  1917.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  1918.   Ii = Array(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)
  1919.   Ji = Array(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)
  1920.   ni = Array(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)
  1921.   delta = rho / rhoc
  1922.   tau = tc / t
  1923.   fidelta = 0
  1924.   fitau = 0
  1925.   For i = 1 To 39
  1926.     fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i)
  1927.     fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1)
  1928.   Next i
  1929.   fidelta = fidelta + ni(0) / delta
  1930.   h3_rhoT = R * t * (tau * fitau + delta * fidelta)
  1931. End Function
  1932. Function s3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1933.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1934.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1935.  'Table 30 and 31, Page 30 and 31
  1936.  Dim i As Integer, Ji, Ii, ni As Variant, fi, delta, tau, fitau As Double
  1937.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  1938.   Ii = Array(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)
  1939.   Ji = Array(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)
  1940.   ni = Array(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)
  1941.   delta = rho / rhoc
  1942.   tau = tc / t
  1943.   fi = 0
  1944.   fitau = 0
  1945.   For i = 1 To 39
  1946.     fi = fi + ni(i) * delta ^ Ii(i) * tau ^ Ji(i)
  1947.     fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1)
  1948.   Next i
  1949.   fi = fi + ni(0) * Log(delta)
  1950.   s3_rhoT = R * (tau * fitau - fi)
  1951. End Function
  1952. Function Cp3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1953.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1954.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1955.  'Table 30 and 31, Page 30 and 31
  1956.  Dim i As Integer, Ji, Ii, ni As Variant, fideltatau, fi, delta, tau, fitautau, fidelta, fideltatautau, fideltadelta As Double
  1957.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  1958.   Ii = Array(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)
  1959.   Ji = Array(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)
  1960.   ni = Array(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)
  1961.   delta = rho / rhoc
  1962.   tau = tc / t
  1963.   fitautau = 0
  1964.   fidelta = 0
  1965.   fideltatau = 0
  1966.   fideltadelta = 0
  1967.   For i = 1 To 39
  1968.     fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2)
  1969.     fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i)
  1970.     fideltatau = fideltatau + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * Ji(i) * tau ^ (Ji(i) - 1)
  1971.     fideltadelta = fideltadelta + ni(i) * Ii(i) * (Ii(i) - 1) * delta ^ (Ii(i) - 2) * tau ^ Ji(i)
  1972.   Next i
  1973.   fidelta = fidelta + ni(0) / delta
  1974.   fideltadelta = fideltadelta - ni(0) / (delta ^ 2)
  1975.   Cp3_rhoT = R * (-(tau ^ 2 * fitautau) + (delta * fidelta - delta * tau * fideltatau) ^ 2 / (2 * delta * fidelta + delta ^ 2 * fideltadelta))
  1976. End Function
  1977. Function Cv3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1978.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1979.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1980.  'Table 30 and 31, Page 30 and 31
  1981.  Dim i As Integer, Ji, Ii, ni As Variant, fi, delta, tau, fitautau As Double
  1982.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  1983.   Ii = Array(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)
  1984.   Ji = Array(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)
  1985.   ni = Array(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)
  1986.   delta = rho / rhoc
  1987.   tau = tc / t
  1988.   fitautau = 0
  1989.   For i = 1 To 39
  1990.     fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2)
  1991.   Next i
  1992.   Cv3_rhoT = R * -(tau * tau * fitautau)
  1993. End Function
  1994. Function w3_rhoT(ByVal rho As Double, ByVal t As Double) As Double
  1995.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  1996.  '7 Basic Equation for Region 3, Section. 6.1 Basic Equation
  1997.  'Table 30 and 31, Page 30 and 31
  1998.  Dim i As Integer, Ji, Ii, ni As Variant, fi, delta, tau, fitautau, fidelta, fideltatau, fideltadelta As Double
  1999.   Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  2000.   Ii = Array(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)
  2001.   Ji = Array(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)
  2002.   ni = Array(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)
  2003.   delta = rho / rhoc
  2004.   tau = tc / t
  2005.   fitautau = 0
  2006.   fidelta = 0
  2007.   fideltatau = 0
  2008.   fideltadelta = 0
  2009.   For i = 1 To 39
  2010.     fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2)
  2011.     fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i)
  2012.     fideltatau = fideltatau + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * Ji(i) * tau ^ (Ji(i) - 1)
  2013.     fideltadelta = fideltadelta + ni(i) * Ii(i) * (Ii(i) - 1) * delta ^ (Ii(i) - 2) * tau ^ Ji(i)
  2014.   Next i
  2015.   fidelta = fidelta + ni(0) / delta
  2016.   fideltadelta = fideltadelta - ni(0) / (delta ^ 2)
  2017.   w3_rhoT = (1000 * R * t * (2 * delta * fidelta + delta ^ 2 * fideltadelta - (delta * fidelta - delta * tau * fideltatau) ^ 2 / (tau ^ 2 * fitautau))) ^ 0.5
  2018. End Function
  2019. Function T3_ph(ByVal p As Double, ByVal h As Double) As Double
  2020. '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
  2021. '2004
  2022. 'Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
  2023. 'Boundary equation, Eq 1 Page 5
  2024. Dim i As Integer, Ji, Ii, ni As Variant, h3ab, ps, hs, Ts As Double
  2025.  Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  2026.  h3ab = (2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3)
  2027.   If h < h3ab Then
  2028.     'Subregion 3a
  2029.    'Eq 2, Table 3, Page 7
  2030.    Ii = Array(-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)
  2031.     Ji = Array(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)
  2032.     ni = Array(-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)
  2033.     ps = p / 100
  2034.     hs = h / 2300
  2035.     Ts = 0
  2036.     For i = 0 To 30
  2037.       Ts = Ts + ni(i) * (ps + 0.24) ^ Ii(i) * (hs - 0.615) ^ Ji(i)
  2038.     Next i
  2039.     T3_ph = Ts * 760
  2040.   Else
  2041.     'Subregion 3b
  2042.    'Eq 3, Table 4, Page 7,8
  2043.    Ii = Array(-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)
  2044.     Ji = Array(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)
  2045.     ni = Array(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)
  2046.     hs = h / 2800
  2047.     ps = p / 100
  2048.     Ts = 0
  2049.     For i = 0 To 32
  2050.       Ts = Ts + ni(i) * (ps + 0.298) ^ Ii(i) * (hs - 0.72) ^ Ji(i)
  2051.     Next i
  2052.     T3_ph = Ts * 860
  2053.   End If
  2054. End Function
  2055. Function v3_ph(ByVal p As Double, ByVal h As Double) As Double
  2056. '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
  2057. '2004
  2058. 'Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
  2059. 'Boundary equation, Eq 1 Page 5
  2060. Dim i As Integer, Ji, Ii, ni As Variant, h3ab, ps, hs, vs As Double
  2061.  Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  2062.   h3ab = (2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3)
  2063.   If h < h3ab Then
  2064.     'Subregion 3a
  2065.    'Eq 4, Table 6, Page 9
  2066.    Ii = Array(-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)
  2067.     Ji = Array(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)
  2068.     ni = Array(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)
  2069.     ps = p / 100
  2070.     hs = h / 2100
  2071.     vs = 0
  2072.     For i = 0 To 31
  2073.       vs = vs + ni(i) * (ps + 0.128) ^ Ii(i) * (hs - 0.727) ^ Ji(i)
  2074.     Next i
  2075.     v3_ph = vs * 0.0028
  2076.   Else
  2077.     'Subregion 3b
  2078.    'Eq 5, Table 7, Page 9
  2079.    Ii = Array(-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)
  2080.     Ji = Array(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)
  2081.     ni = Array(-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)
  2082.     ps = p / 100
  2083.     hs = h / 2800
  2084.     vs = 0
  2085.     For i = 0 To 29
  2086.       vs = vs + ni(i) * (ps + 0.0661) ^ Ii(i) * (hs - 0.72) ^ Ji(i)
  2087.     Next i
  2088.     v3_ph = vs * 0.0088
  2089.   End If
  2090. End Function
  2091. Function T3_ps(ByVal p As Double, ByVal s As Double) As Double
  2092. '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
  2093. '2004
  2094. '3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
  2095. 'Boundary equation, Eq 6 Page 11
  2096. Dim i As Integer, Ji, Ii, ni As Variant, ps, sigma, teta As Double
  2097.  Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  2098.  
  2099.   If s <= 4.41202148223476 Then
  2100.     'Subregion 3a
  2101.    'Eq 6, Table 10, Page 11
  2102.    Ii = Array(-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)
  2103.     Ji = Array(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)
  2104.     ni = Array(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)
  2105.     sigma = s / 4.4
  2106.     ps = p / 100
  2107.     teta = 0
  2108.     For i = 0 To 32
  2109.       teta = teta + ni(i) * (ps + 0.24) ^ Ii(i) * (sigma - 0.703) ^ Ji(i)
  2110.     Next i
  2111.     T3_ps = teta * 760
  2112.   Else
  2113.     'Subregion 3b
  2114.    'Eq 7, Table 11, Page 11
  2115.    Ii = Array(-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)
  2116.     Ji = Array(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)
  2117.     ni = Array(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)
  2118.     sigma = s / 5.3
  2119.     ps = p / 100
  2120.     teta = 0
  2121.     For i = 0 To 27
  2122.       teta = teta + ni(i) * (ps + 0.76) ^ Ii(i) * (sigma - 0.818) ^ Ji(i)
  2123.     Next i
  2124.     T3_ps = teta * 860
  2125.   End If
  2126. End Function
  2127. Function v3_ps(ByVal p As Double, ByVal s As Double) As Double
  2128. '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
  2129. '2004
  2130. '3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
  2131. 'Boundary equation, Eq 6 Page 11
  2132. Dim i As Integer, Ji, Ii, ni As Variant, ps, sigma, omega As Double
  2133.  Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  2134.   If s <= 4.41202148223476 Then
  2135.     'Subregion 3a
  2136.    'Eq 8, Table 13, Page 14
  2137.    Ii = Array(-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)
  2138.     Ji = Array(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)
  2139.     ni = Array(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)
  2140.     ps = p / 100
  2141.     sigma = s / 4.4
  2142.     omega = 0
  2143.     For i = 0 To 27
  2144.       omega = omega + ni(i) * (ps + 0.187) ^ Ii(i) * (sigma - 0.755) ^ Ji(i)
  2145.     Next i
  2146.     v3_ps = omega * 0.0028
  2147.   Else
  2148.     'Subregion 3b
  2149.    'Eq 9, Table 14, Page 14
  2150.    Ii = Array(-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)
  2151.     Ji = Array(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)
  2152.     ni = Array(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)
  2153.     ps = p / 100
  2154.     sigma = s / 5.3
  2155.     omega = 0
  2156.     For i = 0 To 30
  2157.       omega = omega + ni(i) * (ps + 0.298) ^ Ii(i) * (sigma - 0.816) ^ Ji(i)
  2158.    Next i
  2159.    v3_ps = omega * 0.0088
  2160.   End If
  2161. End Function
  2162.  
  2163. Function p3_hs(ByVal h As Double, ByVal s As Double) As Double
  2164. 'Supplementary Release on Backward Equations ( ) , p h s for Region 3,
  2165. 'Equations as a Function of h and s for the Region Boundaries, and an Equation
  2166. '( ) sat , T hs for Region 4 of the IAPWS Industrial Formulation 1997 for the
  2167. 'Thermodynamic Properties of Water and Steam
  2168. '2004
  2169. 'Section 3 Backward Functions p(h,s), T(h,s), and v(h,s) for Region 3
  2170. Dim i As Integer, Ji, Ii, ni As Variant, ps, sigma, eta As Double
  2171.  Const R As Double = 0.461526, tc As Double = 647.096, pc As Double = 22.064, rhoc As Double = 322
  2172.   If s < 4.41202148223476 Then
  2173.     'Subregion 3a
  2174.    'Eq 1, Table 3, Page 8
  2175.    Ii = Array(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)
  2176.     Ji = Array(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)
  2177.     ni = Array(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)
  2178.     sigma = s / 4.4
  2179.     eta = h / 2300
  2180.     ps = 0
  2181.     For i = 0 To 32
  2182.       ps = ps + ni(i) * (eta - 1.01) ^ Ii(i) * (sigma - 0.75) ^ Ji(i)
  2183.     Next i
  2184.     p3_hs = ps * 99
  2185.   Else
  2186.     'Subregion 3b
  2187.    'Eq 2, Table 4, Page 8
  2188.    Ii = Array(-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)
  2189.     Ji = Array(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)
  2190.     ni = Array(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)
  2191.     sigma = s / 5.3
  2192.     eta = h / 2800
  2193.     ps = 0
  2194.     For i = 0 To 34
  2195.       ps = ps + ni(i) * (eta - 0.681) ^ Ii(i) * (sigma - 0.792) ^ Ji(i)
  2196.     Next i
  2197.     p3_hs = 16.6 / ps
  2198.   End If
  2199. End Function
  2200. Function h3_pT(ByVal p As Double, ByVal t As Double) As Double
  2201.    'Not avalible with IF 97
  2202.   'Solve function T3_ph-T=0 with half interval method.
  2203.   Dim Ts, Low_Bound, High_Bound, hs As Double
  2204.    'ver2.6 Start corrected bug
  2205.   If p < 22.06395 Then   'Bellow tripple point
  2206.     Ts = T4_p(p)    'Saturation temperature
  2207.     If t <= Ts Then   'Liquid side
  2208.       High_Bound = h4L_p(p) 'Max h är liauid h.
  2209.       Low_Bound = h1_pT(p, 623.15)
  2210.      Else
  2211.        Low_Bound = h4V_p(p)  'Min h är Vapour h.
  2212.       High_Bound = h2_pT(p, B23T_p(p))
  2213.      End If
  2214.    Else                  'Above tripple point. R3 from R2 till R3.
  2215.     Low_Bound = h1_pT(p, 623.15)
  2216.      High_Bound = h2_pT(p, B23T_p(p))
  2217.    End If
  2218.    'ver2.6 End corrected bug
  2219.   Ts = t + 1
  2220.    Do While Abs(t - Ts) > 0.000001
  2221.     hs = (Low_Bound + High_Bound) / 2
  2222.     Ts = T3_ph(p, hs)
  2223.     If Ts > t Then
  2224.       High_Bound = hs
  2225.     Else
  2226.       Low_Bound = hs
  2227.     End If
  2228.    Loop
  2229.    h3_pT = hs
  2230. End Function
  2231. Function T3_prho(ByVal p As Double, ByVal rho As Double)
  2232.   'Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  2233.  'Solve with half interval method
  2234.  Dim ps, Low_Bound, High_Bound, Ts As Double
  2235.   Low_Bound = 623.15
  2236.   High_Bound = 1073.15
  2237.   Do While Abs(p - ps) > 0.00000001
  2238.     Ts = (Low_Bound + High_Bound) / 2
  2239.     ps = p3_rhoT(rho, Ts)
  2240.     If ps > p Then
  2241.       High_Bound = Ts
  2242.     Else
  2243.       Low_Bound = Ts
  2244.     End If
  2245.     Loop
  2246.     T3_prho = Ts
  2247. End Function
  2248.  
  2249. '***********************************************************************************************************
  2250. '*2.4 Functions for region 4
  2251. Function p4_T(ByVal t As Double) As Double
  2252.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2253.  'Section 8.1 The Saturation-Pressure Equation
  2254.  'Eq 30, Page 33
  2255.  Dim teta, a, b, c As Double
  2256.   teta = t - 0.23855557567849 / (t - 650.17534844798)
  2257.   a = teta ^ 2 + 1167.0521452767 * teta - 724213.16703206
  2258.   b = -17.073846940092 * teta ^ 2 + 12020.82470247 * teta - 3232555.0322333
  2259.   c = 14.91510861353 * teta ^ 2 - 4823.2657361591 * teta + 405113.40542057
  2260.   p4_T = (2 * c / (-b + (b ^ 2 - 4 * a * c) ^ 0.5)) ^ 4
  2261. End Function
  2262.  
  2263. Function T4_p(ByVal p As Double) As Double
  2264.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2265.  'Section 8.2 The Saturation-Temperature Equation
  2266.  'Eq 31, Page 34
  2267.  Dim beta, e, f, g, d As Double
  2268.   beta = p ^ 0.25
  2269.   e = beta ^ 2 - 17.073846940092 * beta + 14.91510861353
  2270.   f = 1167.0521452767 * beta ^ 2 + 12020.82470247 * beta - 4823.2657361591
  2271.   g = -724213.16703206 * beta ^ 2 - 3232555.0322333 * beta + 405113.40542057
  2272.   d = 2 * g / (-f - (f ^ 2 - 4 * e * g) ^ 0.5)
  2273.   T4_p = (650.17534844798 + d - ((650.17534844798 + d) ^ 2 - 4 * (-0.23855557567849 + 650.17534844798 * d)) ^ 0.5) / 2
  2274. End Function
  2275. Function h4_s(ByVal s As Double) As Double
  2276. '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
  2277. '4 Equations for Region Boundaries Given Enthalpy and Entropy
  2278. ' Se picture page 14
  2279. Dim Ii, Ji, ni As Variant, eta, sigma, sigma1, sigma2 As Double, i As Integer
  2280.   If s > -0.0001545495919 And s <= 3.77828134 Then
  2281.     'hL1_s
  2282.    'Eq 3,Table 9,Page 16
  2283.    Ii = Array(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)
  2284.     Ji = Array(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)
  2285.     ni = Array(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)
  2286.     sigma = s / 3.8
  2287.     eta = 0
  2288.     For i = 0 To 26
  2289.       eta = eta + ni(i) * (sigma - 1.09) ^ Ii(i) * (sigma + 0.0000366) ^ Ji(i)
  2290.     Next i
  2291.     h4_s = eta * 1700
  2292.   ElseIf s > 3.77828134 And s <= 4.41202148223476 Then
  2293.     'hL3_s
  2294.    'Eq 4,Table 10,Page 16
  2295.    Ii = Array(0, 0, 0, 0, 2, 3, 4, 4, 5, 5, 6, 7, 7, 7, 10, 10, 10, 32, 32)
  2296.     Ji = Array(1, 4, 10, 16, 1, 36, 3, 16, 20, 36, 4, 2, 28, 32, 14, 32, 36, 0, 6)
  2297.     ni = Array(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)
  2298.     sigma = s / 3.8
  2299.     eta = 0
  2300.     For i = 0 To 18
  2301.       eta = eta + ni(i) * (sigma - 1.09) ^ Ii(i) * (sigma + 0.0000366) ^ Ji(i)
  2302.     Next i
  2303.     h4_s = eta * 1700
  2304.   ElseIf s > 4.41202148223476 And s <= 5.85 Then
  2305.     'Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
  2306.    'Page 19, Eq 5
  2307.    'hV2c3b_s(s)
  2308.    Ii = Array(0, 0, 0, 1, 1, 5, 6, 7, 8, 8, 12, 16, 22, 22, 24, 36)
  2309.     Ji = Array(0, 3, 4, 0, 12, 36, 12, 16, 2, 20, 32, 36, 2, 32, 7, 20)
  2310.     ni = Array(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)
  2311.     sigma = s / 5.9
  2312.     eta = 0
  2313.     For i = 0 To 15
  2314.       eta = eta + ni(i) * (sigma - 1.02) ^ Ii(i) * (sigma - 0.726) ^ Ji(i)
  2315.     Next i
  2316.     h4_s = eta ^ 4 * 2800
  2317.   ElseIf s > 5.85 And s < 9.155759395 Then
  2318.     'Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
  2319.    'Page 20, Eq 6
  2320.    Ii = Array(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)
  2321.     Ji = Array(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)
  2322.     ni = Array(-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)
  2323.     sigma1 = s / 5.21
  2324.     sigma2 = s / 9.2
  2325.     eta = 0
  2326.     For i = 0 To 29
  2327.       eta = eta + ni(i) * (1 / sigma1 - 0.513) ^ Ii(i) * (sigma2 - 0.524) ^ Ji(i)
  2328.     Next i
  2329.     h4_s = Exp(eta) * 2800
  2330.   Else
  2331.     h4_s = CVErr(xlErrValue)
  2332.   End If
  2333. End Function
  2334. Function p4_s(ByVal s As Double) As Double
  2335.   'Uses h4_s and p_hs for the diffrent regions to determine p4_s
  2336.  Dim hsat As Double
  2337.   hsat = h4_s(s)
  2338.   If s > -0.0001545495919 And s <= 3.77828134 Then
  2339.     p4_s = p1_hs(hsat, s)
  2340.   ElseIf s > 3.77828134 And s <= 5.210887663 Then
  2341.     p4_s = p3_hs(hsat, s)
  2342.   ElseIf s > 5.210887663 And s < 9.155759395 Then
  2343.     p4_s = p2_hs(hsat, s)
  2344.   Else
  2345.     p4_s = CVErr(xlErrValue)
  2346.   End If
  2347. End Function
  2348. Function h4L_p(ByVal p As Double) As Double
  2349.  Dim Low_Bound, High_Bound, hs, ps, Ts As Double
  2350.  If p > 0.000611657 And p < 22.06395 Then
  2351.   Ts = T4_p(p)
  2352.   If p < 16.529 Then
  2353.     h4L_p = h1_pT(p, Ts)
  2354.   Else
  2355.     'Iterate to find the the backward solution of p3sat_h
  2356.    Low_Bound = 1670.858218
  2357.     High_Bound = 2087.23500164864
  2358.     Do While Abs(p - ps) > 0.00001
  2359.       hs = (Low_Bound + High_Bound) / 2
  2360.       ps = p3sat_h(hs)
  2361.       If ps > p Then
  2362.         High_Bound = hs
  2363.       Else
  2364.         Low_Bound = hs
  2365.       End If
  2366.     Loop
  2367.    
  2368.     h4L_p = hs
  2369.   End If
  2370.  Else
  2371.   h4L_p = CVErr(xlErrValue)
  2372.  End If
  2373. End Function
  2374. Function h4V_p(ByVal p As Double) As Double
  2375.  Dim Low_Bound, High_Bound, hs, ps, Ts As Double
  2376.  If p > 0.000611657 And p < 22.06395 Then
  2377.   Ts = T4_p(p)
  2378.   If p < 16.529 Then
  2379.     h4V_p = h2_pT(p, Ts)
  2380.   Else
  2381.     'Iterate to find the the backward solution of p3sat_h
  2382.    Low_Bound = 2087.23500164864
  2383.     High_Bound = 2563.592004 + 5 '5 added to extrapolate to ensure even the border ==350°C solved.
  2384.    Do While Abs(p - ps) > 0.000001
  2385.       hs = (Low_Bound + High_Bound) / 2
  2386.       ps = p3sat_h(hs)
  2387.       If ps < p Then
  2388.         High_Bound = hs
  2389.       Else
  2390.         Low_Bound = hs
  2391.       End If
  2392.     Loop
  2393.     h4V_p = hs
  2394.   End If
  2395.  Else
  2396.   h4V_p = CVErr(xlErrValue)
  2397.  End If
  2398. End Function
  2399. Function x4_ph(ByVal p As Double, ByVal h As Double) As Double
  2400. 'Calculate vapour fraction from hL and hV for given p
  2401. Dim h4v, h4l As Double
  2402.   h4v = h4V_p(p)
  2403.   h4l = h4L_p(p)
  2404.   If h > h4v Then
  2405.     x4_ph = 1
  2406.   ElseIf h < h4l Then
  2407.     x4_ph = 0
  2408.   Else
  2409.     x4_ph = (h - h4l) / (h4v - h4l)
  2410.   End If
  2411. End Function
  2412. Function x4_ps(ByVal p As Double, ByVal s As Double) As Double
  2413.  Dim ssV, ssL As Double
  2414.   If p < 16.529 Then
  2415.    ssV = s2_pT(p, T4_p(p))
  2416.    ssL = s1_pT(p, T4_p(p))
  2417.   Else
  2418.    ssV = s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p))
  2419.    ssL = s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p))
  2420.   End If
  2421.   If s < ssL Then
  2422.     x4_ps = 0
  2423.   ElseIf s > ssV Then
  2424.     x4_ps = 1
  2425.   Else
  2426.     x4_ps = (s - ssL) / (ssV - ssL)
  2427.   End If
  2428. End Function
  2429. Function T4_hs(ByVal h As Double, ByVal s As Double) As Double
  2430. 'Supplementary Release on Backward Equations ( ) , p h s for Region 3,
  2431. 'Chapter 5.3 page 30.
  2432. 'The if 97 function is only valid for part of region4. Use iteration outsida.
  2433. Dim Ii, Ji, ni As Variant, hL, Ts, ss, p, sigma, eta, teta, High_Bound, Low_Bound, PL, s4V, v4V, s4L, v4L, xs As Double, i As Integer
  2434.    Ii = Array(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)
  2435.    Ji = Array(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)
  2436.    ni = Array(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)
  2437.   If (s > 5.210887825 And s < 9.15546555571324) Then
  2438.     sigma = s / 9.2
  2439.     eta = h / 2800
  2440.     teta = 0
  2441.     For i = 0 To 35
  2442.       teta = teta + ni(i) * (eta - 0.119) ^ Ii(i) * (sigma - 1.07) ^ Ji(i)
  2443.     Next i
  2444.     T4_hs = teta * 550
  2445. Else
  2446.     'Function psat_h
  2447.    If s > -0.0001545495919 And s <= 3.77828134 Then
  2448.       Low_Bound = 0.000611
  2449.       High_Bound = 165.291642526045
  2450.       Do While Abs(hL - h) > 0.00001 And Abs(High_Bound - Low_Bound) > 0.0001
  2451.        PL = (Low_Bound + High_Bound) / 2
  2452.        Ts = T4_p(PL)
  2453.        hL = h1_pT(PL, Ts)
  2454.        If hL > h Then
  2455.          High_Bound = PL
  2456.        Else
  2457.          Low_Bound = PL
  2458.        End If
  2459.       Loop
  2460.     ElseIf s > 3.77828134 And s <= 4.41202148223476 Then
  2461.       PL = p3sat_h(h)
  2462.     ElseIf s > 4.41202148223476 And s <= 5.210887663 Then
  2463.       PL = p3sat_h(h)
  2464.     End If
  2465.     Low_Bound = 0.000611
  2466.     High_Bound = PL
  2467.     Do While Abs(s - ss) > 0.000001 And Abs(High_Bound - Low_Bound) > 0.0000001
  2468.       p = (Low_Bound + High_Bound) / 2
  2469.      
  2470.       'Calculate s4_ph
  2471.      Ts = T4_p(p)
  2472.       xs = x4_ph(p, h)
  2473.       If p < 16.529 Then
  2474.         s4V = s2_pT(p, Ts)
  2475.         s4L = s1_pT(p, Ts)
  2476.       Else
  2477.         v4V = v3_ph(p, h4V_p(p))
  2478.         s4V = s3_rhoT(1 / v4V, Ts)
  2479.         v4L = v3_ph(p, h4L_p(p))
  2480.         s4L = s3_rhoT(1 / v4L, Ts)
  2481.       End If
  2482.       ss = (xs * s4V + (1 - xs) * s4L)
  2483.      
  2484.       If ss < s Then
  2485.         High_Bound = p
  2486.       Else
  2487.         Low_Bound = p
  2488.       End If
  2489.     Loop
  2490.     T4_hs = T4_p(p)
  2491. End If
  2492. End Function
  2493. '***********************************************************************************************************
  2494. '*2.5 Functions for region 5
  2495. Function h5_pT(ByVal p As Double, ByVal t As Double) As Double
  2496.   'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2497.  'Basic Equation for Region 5
  2498.  'Eq 32,33, Page 36, Tables 37-41
  2499.   Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0_tau, gammar_tau As Double, i As Integer
  2500.    Const R As Double = 0.461526   'kJ/(kg K)
  2501.  Ji0 = Array(0, 1, -3, -2, -1, 2)
  2502.   ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2503.   Iir = Array(1, 1, 1, 2, 3)
  2504.   Jir = Array(0, 1, 3, 9, 3)
  2505.   nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2506.   tau = 1000 / t
  2507.   gamma0_tau = 0
  2508.   For i = 0 To 5
  2509.     gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1)
  2510.   Next i
  2511.   gammar_tau = 0
  2512.   For i = 0 To 4
  2513.     gammar_tau = gammar_tau + nir(i) * p ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1)
  2514.   Next i
  2515.   h5_pT = R * t * tau * (gamma0_tau + gammar_tau)
  2516. End Function
  2517.  
  2518.  
  2519. Function v5_pT(ByVal p As Double, ByVal t As Double) As Double
  2520. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2521. 'Basic Equation for Region 5
  2522. 'Eq 32,33, Page 36, Tables 37-41
  2523. Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0_pi, gammar_pi As Double, i As Integer
  2524. Const R As Double = 0.461526   'kJ/(kg K)
  2525. Ji0 = Array(0, 1, -3, -2, -1, 2)
  2526. ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2527. Iir = Array(1, 1, 1, 2, 3)
  2528. Jir = Array(0, 1, 3, 9, 3)
  2529. nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2530. tau = 1000 / t
  2531. gamma0_pi = 1 / p
  2532. gammar_pi = 0
  2533. For i = 0 To 4
  2534.   gammar_pi = gammar_pi + nir(i) * Iir(i) * p ^ (Iir(i) - 1) * tau ^ Jir(i)
  2535. Next i
  2536. v5_pT = R * t / p * p * (gamma0_pi + gammar_pi) / 1000
  2537. End Function
  2538.  
  2539. Function u5_pT(ByVal p As Double, ByVal t As Double) As Double
  2540. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2541. 'Basic Equation for Region 5
  2542. 'Eq 32,33, Page 36, Tables 37-41
  2543. Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0_pi, gammar_pi, gamma0_tau, gammar_tau As Double, i As Integer
  2544. Const R As Double = 0.461526   'kJ/(kg K)
  2545. Ji0 = Array(0, 1, -3, -2, -1, 2)
  2546. ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2547. Iir = Array(1, 1, 1, 2, 3)
  2548. Jir = Array(0, 1, 3, 9, 3)
  2549. nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2550. tau = 1000 / t
  2551. gamma0_pi = 1 / p
  2552. gamma0_tau = 0
  2553. For i = 0 To 5
  2554.   gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1)
  2555. Next i
  2556. gammar_pi = 0
  2557. gammar_tau = 0
  2558. For i = 0 To 4
  2559.   gammar_pi = gammar_pi + nir(i) * Iir(i) * p ^ (Iir(i) - 1) * tau ^ Jir(i)
  2560.   gammar_tau = gammar_tau + nir(i) * p ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1)
  2561. Next i
  2562. u5_pT = R * t * (tau * (gamma0_tau + gammar_tau) - p * (gamma0_pi + gammar_pi))
  2563. End Function
  2564. Function Cp5_pT(ByVal p As Double, ByVal t As Double) As Double
  2565. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2566. 'Basic Equation for Region 5
  2567. 'Eq 32,33, Page 36, Tables 37-41
  2568. Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0_tautau, gammar_tautau As Double, i As Integer
  2569. Const R As Double = 0.461526   'kJ/(kg K)
  2570. Ji0 = Array(0, 1, -3, -2, -1, 2)
  2571. ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2572. Iir = Array(1, 1, 1, 2, 3)
  2573. Jir = Array(0, 1, 3, 9, 3)
  2574. nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2575. tau = 1000 / t
  2576. gamma0_tautau = 0
  2577. For i = 0 To 5
  2578.   gamma0_tautau = gamma0_tautau + ni0(i) * Ji0(i) * (Ji0(i) - 1) * tau ^ (Ji0(i) - 2)
  2579. Next i
  2580. gammar_tautau = 0
  2581. For i = 0 To 4
  2582.   gammar_tautau = gammar_tautau + nir(i) * p ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2)
  2583. Next i
  2584. Cp5_pT = -R * tau ^ 2 * (gamma0_tautau + gammar_tautau)
  2585. End Function
  2586.  
  2587. Function s5_pT(ByVal p As Double, ByVal t As Double) As Double
  2588. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2589. 'Basic Equation for Region 5
  2590. 'Eq 32,33, Page 36, Tables 37-41
  2591. Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0, gamma0_tau, gammar, gammar_tau As Double, i As Integer
  2592. Const R As Double = 0.461526   'kJ/(kg K)
  2593. Ji0 = Array(0, 1, -3, -2, -1, 2)
  2594. ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2595. Iir = Array(1, 1, 1, 2, 3)
  2596. Jir = Array(0, 1, 3, 9, 3)
  2597. nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2598. tau = 1000 / t
  2599. gamma0 = Log(p)
  2600. gamma0_tau = 0
  2601. For i = 0 To 5
  2602.   gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1)
  2603.   gamma0 = gamma0 + ni0(i) * tau ^ Ji0(i)
  2604. Next i
  2605. gammar = 0
  2606. gammar_tau = 0
  2607. For i = 0 To 4
  2608.   gammar = gammar + nir(i) * p ^ Iir(i) * tau ^ Jir(i)
  2609.   gammar_tau = gammar_tau + nir(i) * p ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1)
  2610. Next i
  2611. s5_pT = R * (tau * (gamma0_tau + gammar_tau) - (gamma0 + gammar))
  2612. End Function
  2613. Function Cv5_pT(ByVal p As Double, ByVal t As Double) As Double
  2614. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2615. 'Basic Equation for Region 5
  2616. 'Eq 32,33, Page 36, Tables 37-41
  2617. Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0_tautau, gammar_pi, gammar_pitau, gammar_pipi, gammar_tautau As Double, i As Integer
  2618. Const R As Double = 0.461526   'kJ/(kg K)
  2619. Ji0 = Array(0, 1, -3, -2, -1, 2)
  2620. ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2621. Iir = Array(1, 1, 1, 2, 3)
  2622. Jir = Array(0, 1, 3, 9, 3)
  2623. nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2624. tau = 1000 / t
  2625. gamma0_tautau = 0
  2626. For i = 0 To 5
  2627.   gamma0_tautau = gamma0_tautau + ni0(i) * (Ji0(i) - 1) * Ji0(i) * tau ^ (Ji0(i) - 2)
  2628. Next i
  2629. gammar_pi = 0
  2630. gammar_pitau = 0
  2631. gammar_pipi = 0
  2632. gammar_tautau = 0
  2633. For i = 0 To 4
  2634.   gammar_pi = gammar_pi + nir(i) * Iir(i) * p ^ (Iir(i) - 1) * tau ^ Jir(i)
  2635.   gammar_pitau = gammar_pitau + nir(i) * Iir(i) * p ^ (Iir(i) - 1) * Jir(i) * tau ^ (Jir(i) - 1)
  2636.   gammar_pipi = gammar_pipi + nir(i) * Iir(i) * (Iir(i) - 1) * p ^ (Iir(i) - 2) * tau ^ Jir(i)
  2637.   gammar_tautau = gammar_tautau + nir(i) * p ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2)
  2638. Next i
  2639. Cv5_pT = R * (-(tau ^ 2 * (gamma0_tautau + gammar_tautau)) - (1 + p * gammar_pi - tau * p * gammar_pitau) ^ 2 / (1 - p ^ 2 * gammar_pipi))
  2640.  
  2641. End Function
  2642. Function w5_pT(ByVal p As Double, ByVal t As Double) As Double
  2643. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
  2644. 'Basic Equation for Region 5
  2645. 'Eq 32,33, Page 36, Tables 37-41
  2646. Dim Iir, Jir, nir, ni0, Ji0 As Variant, tau, gamma0_tautau, gammar_pi, gammar_pitau, gammar_pipi, gammar_tautau As Double, i As Integer
  2647. Const R As Double = 0.461526   'kJ/(kg K)
  2648. Ji0 = Array(0, 1, -3, -2, -1, 2)
  2649. ni0 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)
  2650. Iir = Array(1, 1, 1, 2, 3)
  2651. Jir = Array(0, 1, 3, 9, 3)
  2652. nir = Array(-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07)
  2653. tau = 1000 / t
  2654. gamma0_tautau = 0
  2655. For i = 0 To 5
  2656.   gamma0_tautau = gamma0_tautau + ni0(i) * (Ji0(i) - 1) * Ji0(i) * tau ^ (Ji0(i) - 2)
  2657. Next i
  2658. gammar_pi = 0
  2659. gammar_pitau = 0
  2660. gammar_pipi = 0
  2661. gammar_tautau = 0
  2662. For i = 0 To 4
  2663.   gammar_pi = gammar_pi + nir(i) * Iir(i) * p ^ (Iir(i) - 1) * tau ^ Jir(i)
  2664.   gammar_pitau = gammar_pitau + nir(i) * Iir(i) * p ^ (Iir(i) - 1) * Jir(i) * tau ^ (Jir(i) - 1)
  2665.   gammar_pipi = gammar_pipi + nir(i) * Iir(i) * (Iir(i) - 1) * p ^ (Iir(i) - 2) * tau ^ Jir(i)
  2666.   gammar_tautau = gammar_tautau + nir(i) * p ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2)
  2667. Next i
  2668. w5_pT = (1000 * R * t * (1 + 2 * p * gammar_pi + p ^ 2 * gammar_pi ^ 2) / ((1 - p ^ 2 * gammar_pipi) + (1 + p * gammar_pi - tau * p * gammar_pitau) ^ 2 / (tau ^ 2 * (gamma0_tautau + gammar_tautau)))) ^ 0.5
  2669. End Function
  2670.  
  2671. Function T5_ph(ByVal p As Double, ByVal h As Double) As Double
  2672.     'Solve with half interval method
  2673.    Dim Low_Bound, High_Bound, Ts, hs As Double
  2674.     Low_Bound = 1073.15
  2675.     High_Bound = 2273.15
  2676.     Do While Abs(h - hs) > 0.00001
  2677.       Ts = (Low_Bound + High_Bound) / 2
  2678.       hs = h5_pT(p, Ts)
  2679.       If hs > h Then
  2680.         High_Bound = Ts
  2681.       Else
  2682.         Low_Bound = Ts
  2683.       End If
  2684.     Loop
  2685.     T5_ph = Ts
  2686. End Function
  2687.  
  2688. Function T5_ps(ByVal p As Double, ByVal s As Double) As Double
  2689.     'Solve with half interval method
  2690.    Dim Low_Bound, High_Bound, Ts, ss As Double
  2691.     Low_Bound = 1073.15
  2692.     High_Bound = 2273.15
  2693.     Do While Abs(s - ss) > 0.00001
  2694.       Ts = (Low_Bound + High_Bound) / 2
  2695.       ss = s5_pT(p, Ts)
  2696.       If ss > s Then
  2697.         High_Bound = Ts
  2698.       Else
  2699.         Low_Bound = Ts
  2700.       End If
  2701.     Loop
  2702.     T5_ps = Ts
  2703. End Function
  2704. Function T5_prho(ByVal p As Double, ByVal rho As Double) As Double
  2705.   'Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
  2706.  'Solve with half interval method
  2707.  Dim Low_Bound, High_Bound, Ts, rhos As Double
  2708.     Low_Bound = 1073.15
  2709.     High_Bound = 2073.15
  2710.   Do While Abs(rho - rhos) > 0.000001
  2711.     Ts = (Low_Bound + High_Bound) / 2
  2712.     rhos = 1 / v2_pT(p, Ts)
  2713.     If rhos < rho Then
  2714.       High_Bound = Ts
  2715.     Else
  2716.       Low_Bound = Ts
  2717.     End If
  2718.     Loop
  2719.     T5_prho = Ts
  2720. End Function
  2721. '***********************************************************************************************************
  2722. '*3 Region Selection
  2723. '***********************************************************************************************************
  2724. '*3.1 Regions as a function of pT
  2725. Function region_pT(ByVal p As Double, ByVal t As Double) As Integer
  2726. Dim ps As Double
  2727. If t > 1073.15 And p < 10 And t < 2273.15 And p > 0.000611 Then
  2728.   region_pT = 5
  2729. ElseIf t <= 1073.15 And t > 273.15 And p <= 100 And p > 0.000611 Then
  2730.   If t > 623.15 Then
  2731.     If p > B23p_T(t) Then
  2732.      region_pT = 3
  2733.      If t < 647.096 Then
  2734.         ps = p4_T(t)
  2735.         If Abs(p - ps) < 0.00001 Then
  2736.            region_pT = 4
  2737.         End If
  2738.      End If
  2739.     Else
  2740.      region_pT = 2
  2741.     End If
  2742.    Else
  2743.     ps = p4_T(t)
  2744.     If Abs(p - ps) < 0.00001 Then
  2745.       region_pT = 4
  2746.     ElseIf p > ps Then
  2747.       region_pT = 1
  2748.     Else
  2749.       region_pT = 2
  2750.     End If
  2751.    End If
  2752.   Else
  2753.     region_pT = 0 '**Error, Outside valid area
  2754.  End If
  2755. End Function
  2756. '***********************************************************************************************************
  2757. '*3.2 Regions as a function of ph
  2758. Function region_ph(ByVal p, ByVal h) As Integer
  2759. Dim hL, hV, h_45, h_5u, Ts As Double
  2760.  'Check if outside pressure limits
  2761. If p < 0.000611657 Or p > 100 Then
  2762.      region_ph = 0
  2763.      Exit Function
  2764.  End If
  2765.  
  2766.  'Check if outside low h.
  2767. If h < 0.963 * p + 2.2 Then 'Linear adaption to h1_pt()+2 to speed up calcualations.
  2768.    If h < h1_pT(p, 273.15) Then
  2769.       region_ph = 0
  2770.       Exit Function
  2771.     End If
  2772.  End If
  2773.  
  2774.  If p < 16.5292 Then 'Bellow region 3,Check  region 1,4,2,5
  2775.   'Check Region 1
  2776.   Ts = T4_p(p)
  2777.    hL = 109.6635 * Log(p) + 40.3481 * p + 734.58 'Approximate function for hL_p
  2778.   If Abs(h - hL) < 100 Then 'If approximate is not god enough use real function
  2779.      hL = h1_pT(p, Ts)
  2780.    End If
  2781.    If h <= hL Then
  2782.      region_ph = 1
  2783.      Exit Function
  2784.    End If
  2785.    'Check Region 4
  2786.   hV = 45.1768 * Log(p) - 20.158 * p + 2804.4 'Approximate function for hV_p
  2787.   If Abs(h - hV) < 50 Then 'If approximate is not god enough use real function
  2788.      hV = h2_pT(p, Ts)
  2789.    End If
  2790.    If h < hV Then
  2791.      region_ph = 4
  2792.      Exit Function
  2793.    End If
  2794.    'Check upper limit of region 2 Quick Test
  2795.   If h < 4000 Then
  2796.      region_ph = 2
  2797.      Exit Function
  2798.    End If
  2799.   'Check region 2 (Real value)
  2800.   h_45 = h2_pT(p, 1073.15)
  2801.    If h <= h_45 Then
  2802.      region_ph = 2
  2803.      Exit Function
  2804.    End If
  2805.   'Check region 5
  2806.   If p > 10 Then
  2807.      region_ph = 0
  2808.      Exit Function
  2809.    End If
  2810.    h_5u = h5_pT(p, 2273.15)
  2811.    If h < h_5u Then
  2812.       region_ph = 5
  2813.       Exit Function
  2814.    End If
  2815.    region_ph = 0
  2816.    Exit Function
  2817.   Else 'For p>16.5292
  2818.   'Check if in region1
  2819.   If h < h1_pT(p, 623.15) Then
  2820.      region_ph = 1
  2821.      Exit Function
  2822.    End If
  2823.    'Check if in region 3 or 4 (Bellow Reg 2)
  2824.   If h < h2_pT(p, B23T_p(p)) Then
  2825.      'Region 3 or 4
  2826.     If p > p3sat_h(h) Then
  2827.        region_ph = 3
  2828.        Exit Function
  2829.      Else
  2830.        region_ph = 4
  2831.        Exit Function
  2832.      End If
  2833.   End If
  2834.   'Check if region 2
  2835.  If h < h2_pT(p, 1073.15) Then
  2836.     region_ph = 2
  2837.     Exit Function
  2838.   End If
  2839.  End If
  2840.  region_ph = 0
  2841. End Function
  2842. '***********************************************************************************************************
  2843. '*3.3 Regions as a function of ps
  2844. Function region_ps(ByVal p As Double, ByVal s) As Integer
  2845. Dim ss As Double
  2846.   If p < 0.000611657 Or p > 100 Or s < 0 Or s > s5_pT(p, 2273.15) Then
  2847.    region_ps = 0
  2848.    Exit Function
  2849.   End If
  2850.  
  2851.   'Check region 5
  2852.  If s > s2_pT(p, 1073.15) Then
  2853.     If p <= 10 Then
  2854.       region_ps = 5
  2855.       Exit Function
  2856.     Else
  2857.       region_ps = 0
  2858.       Exit Function
  2859.     End If
  2860.   End If
  2861.  
  2862.   'Check region 2
  2863.  If p > 16.529 Then
  2864.     ss = s2_pT(p, B23T_p(p)) 'Between 5.047 and 5.261. Use to speed up!
  2865.  Else
  2866.     ss = s2_pT(p, T4_p(p))
  2867.   End If
  2868.   If s > ss Then
  2869.       region_ps = 2
  2870.       Exit Function
  2871.   End If
  2872.  
  2873.   'Check region 3
  2874.  ss = s1_pT(p, 623.15)
  2875.   If p > 16.529 And s > ss Then
  2876.     If p > p3sat_s(s) Then
  2877.       region_ps = 3
  2878.       Exit Function
  2879.     Else
  2880.       region_ps = 4
  2881.       Exit Function
  2882.     End If
  2883.   End If
  2884.  
  2885.   'Check region 4 (Not inside region 3)
  2886.  If p < 16.529 And s > s1_pT(p, T4_p(p)) Then
  2887.     region_ps = 4
  2888.     Exit Function
  2889.   End If
  2890.  
  2891.   'Check region 1
  2892.  If p > 0.000611657 And s > s1_pT(p, 273.15) Then
  2893.     region_ps = 1
  2894.     Exit Function
  2895.   End If
  2896.   region_ps = 1
  2897. End Function
  2898. '***********************************************************************************************************
  2899. '*3.4 Regions as a function of hs
  2900. Function Region_hs(ByVal h, ByVal s) As Integer
  2901. Dim TMax, hMax, hB, hL, hV, vmax, Tmin, hMin As Double
  2902.   If s < -0.0001545495919 Then
  2903.     Region_hs = 0
  2904.     Exit Function
  2905.   End If
  2906.   'Check linear adaption to p=0.000611. If bellow region 4.
  2907.  hMin = (((-0.0415878 - 2500.89262) / (-0.00015455 - 9.155759)) * s)
  2908.   If s < 9.155759395 And h < hMin Then
  2909.     Region_hs = 0
  2910.     Exit Function
  2911.   End If
  2912.  
  2913.   '******Kolla 1 eller 4. (+liten bit över B13)
  2914.  If s >= -0.0001545495919 And s <= 3.77828134 Then
  2915.     If h < h4_s(s) Then
  2916.       Region_hs = 4
  2917.       Exit Function
  2918.     ElseIf s < 3.397782955 Then '100MPa line is limiting
  2919.      TMax = T1_ps(100, s)
  2920.       hMax = h1_pT(100, TMax)
  2921.       If h < hMax Then
  2922.        Region_hs = 1
  2923.        Exit Function
  2924.       Else
  2925.        Region_hs = 0
  2926.        Exit Function
  2927.       End If
  2928.      Else 'The point is either in region 4,1,3. Check B23
  2929.      hB = hB13_s(s)
  2930.       If h < hB Then
  2931.         Region_hs = 1
  2932.         Exit Function
  2933.       End If
  2934.       TMax = T3_ps(100, s)
  2935.       vmax = v3_ps(100, s)
  2936.       hMax = h3_rhoT(1 / vmax, TMax)
  2937.       If h < hMax Then
  2938.         Region_hs = 3
  2939.         Exit Function
  2940.       Else
  2941.         Region_hs = 0
  2942.         Exit Function
  2943.       End If
  2944.      End If
  2945.   End If
  2946.  
  2947.   '******Kolla region 2 eller 4. (Övre delen av omrĺde b23-> max)
  2948.  If s >= 5.260578707 And s <= 11.9212156897728 Then
  2949.     If s > 9.155759395 Then 'Above region 4
  2950.      Tmin = T2_ps(0.000611, s)
  2951.       hMin = h2_pT(0.000611, Tmin)
  2952.       'Function adapted to h(1073.15,s)
  2953.      hMax = -0.07554022 * s ^ 4 + 3.341571 * s ^ 3 - 55.42151 * s ^ 2 + 408.515 * s + 3031.338
  2954.       If h > hMin And h < hMax Then
  2955.         Region_hs = 2
  2956.         Exit Function
  2957.       Else
  2958.         Region_hs = 0
  2959.         Exit Function
  2960.       End If
  2961.     End If
  2962.    
  2963.    
  2964.       hV = h4_s(s)
  2965.    
  2966.     If h < hV Then  'Region 4. Under region 3.
  2967.        Region_hs = 4
  2968.         Exit Function
  2969.     End If
  2970.     If s < 6.04048367171238 Then
  2971.       TMax = T2_ps(100, s)
  2972.       hMax = h2_pT(100, TMax)
  2973.     Else
  2974.      'Function adapted to h(1073.15,s)
  2975.      hMax = -2.988734 * s ^ 4 + 121.4015 * s ^ 3 - 1805.15 * s ^ 2 + 11720.16 * s - 23998.33
  2976.     End If
  2977.      If h < hMax Then  'Region 2. Över region 4.
  2978.        Region_hs = 2
  2979.         Exit Function
  2980.     Else
  2981.         Region_hs = 0
  2982.         Exit Function
  2983.     End If
  2984.    End If
  2985.    
  2986.    'Kolla region 3 eller 4. Under kritiska punkten.
  2987.   If s >= 3.77828134 And s <= 4.41202148223476 Then
  2988.      hL = h4_s(s)
  2989.      If h < hL Then
  2990.        Region_hs = 4
  2991.        Exit Function
  2992.      End If
  2993.      TMax = T3_ps(100, s)
  2994.      vmax = v3_ps(100, s)
  2995.      hMax = h3_rhoT(1 / vmax, TMax)
  2996.      If h < hMax Then
  2997.         Region_hs = 3
  2998.         Exit Function
  2999.      Else
  3000.         Region_hs = 0
  3001.         Exit Function
  3002.     End If
  3003.    End If
  3004.    
  3005.    'Kolla region 3 eller 4 frĺn kritiska punkten till övre delen av b23
  3006.   If s >= 4.41202148223476 And s <= 5.260578707 Then
  3007.      hV = h4_s(s)
  3008.      If h < hV Then
  3009.         Region_hs = 4
  3010.         Exit Function
  3011.      End If
  3012.      'Kolla om vi är under b23 giltighetsomrĺde.
  3013.     If s <= 5.048096828 Then
  3014.        TMax = T3_ps(100, s)
  3015.        vmax = v3_ps(100, s)
  3016.        hMax = h3_rhoT(1 / vmax, TMax)
  3017.        If h < hMax Then
  3018.          Region_hs = 3
  3019.          Exit Function
  3020.        Else
  3021.          Region_hs = 0
  3022.          Exit Function
  3023.        End If
  3024.      Else 'Inom omrĺdet för B23 i s led.
  3025.       If (h > 2812.942061) Then 'Ovanför B23 i h_led
  3026.         If s > 5.09796573397125 Then
  3027.            TMax = T2_ps(100, s)
  3028.            hMax = h2_pT(100, TMax)
  3029.            If h < hMax Then
  3030.              Region_hs = 2
  3031.              Exit Function
  3032.            Else
  3033.              Region_hs = 0
  3034.              Exit Function
  3035.            End If
  3036.          Else
  3037.            Region_hs = 0
  3038.            Exit Function
  3039.          End If
  3040.        End If
  3041.        If (h < 2563.592004) Then   'Nedanför B23 i h_led men vi har redan kollat ovanför hV2c3b
  3042.          Region_hs = 3
  3043.           Exit Function
  3044.        End If
  3045.        'Vi är inom b23 omrĺdet i bĺde s och h led.
  3046.       If p2_hs(h, s) > B23p_T(TB23_hs(h, s)) Then
  3047.          Region_hs = 3
  3048.          Exit Function
  3049.        Else
  3050.          Region_hs = 2
  3051.          Exit Function
  3052.        End If
  3053.      End If
  3054.    End If
  3055.    Region_hs = CVErr(xlErrValue)
  3056. End Function
  3057. '***********************************************************************************************************
  3058. '*3.5 Regions as a function of p and rho
  3059. Function Region_prho(ByVal p As Double, ByVal rho) As Integer
  3060.   Dim v As Double
  3061.   v = 1 / rho
  3062.   If p < 0.000611657 Or p > 100 Then
  3063.     Region_prho = 0
  3064.     Exit Function
  3065.   End If
  3066.   If p < 16.5292 Then 'Bellow region 3, Check region 1,4,2
  3067.    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.
  3068.      Region_prho = 0
  3069.       Exit Function
  3070.     End If
  3071.     If v <= v1_pT(p, T4_p(p)) Then
  3072.       Region_prho = 1
  3073.       Exit Function
  3074.     End If
  3075.     If v < v2_pT(p, T4_p(p)) Then
  3076.       Region_prho = 4
  3077.       Exit Function
  3078.     End If
  3079.     If v <= v2_pT(p, 1073.15) Then
  3080.       Region_prho = 2
  3081.       Exit Function
  3082.     End If
  3083.     If p > 10 Then 'Above region 5
  3084.      Region_prho = 0
  3085.       Exit Function
  3086.     End If
  3087.     If v <= v5_pT(p, 2073.15) Then
  3088.       Region_prho = 5
  3089.       Exit Function
  3090.     End If
  3091.   Else 'Check region 1,3,4,3,2 (Above the lowest point of region 3.)
  3092.    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.
  3093.      Region_prho = 0
  3094.       Exit Function
  3095.     End If
  3096.     If v < v1_pT(p, 623.15) Then
  3097.       Region_prho = 1
  3098.       Exit Function
  3099.     End If
  3100.     'Check if in region 3 or 4 (Bellow Reg 2)
  3101.    If v < v2_pT(p, B23T_p(p)) Then
  3102.      'Region 3 or 4
  3103.      If p > 22.064 Then 'Above region 4
  3104.        Region_prho = 3
  3105.         Exit Function
  3106.       End If
  3107.       If v < v3_ph(p, h4L_p(p)) Or v > v3_ph(p, h4V_p(p)) Then 'Uses iteration!!
  3108.        Region_prho = 3
  3109.         Exit Function
  3110.       Else
  3111.         Region_prho = 4
  3112.         Exit Function
  3113.       End If
  3114.     End If
  3115.     'Check if region 2
  3116.    If v < v2_pT(p, 1073.15) Then
  3117.       Region_prho = 2
  3118.       Exit Function
  3119.     End If
  3120.   End If
  3121.  
  3122.   Region_prho = 0
  3123. End Function
  3124.  
  3125.  
  3126. '***********************************************************************************************************
  3127. '*4 Region Borders
  3128. '***********************************************************************************************************
  3129. '***********************************************************************************************************
  3130. '*4.1 Boundary between region 2 and 3.
  3131. Function B23p_T(ByVal t As Double) As Double
  3132. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam
  3133. '1997
  3134. 'Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
  3135. 'Eq 5, Page 5
  3136. B23p_T = 348.05185628969 - 1.1671859879975 * t + 1.0192970039326E-03 * t ^ 2
  3137. End Function
  3138. Function B23T_p(ByVal p As Double) As Double
  3139. 'Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam
  3140. '1997
  3141. 'Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
  3142. 'Eq 6, Page 6
  3143. B23T_p = 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ^ 0.5
  3144. End Function
  3145. '***********************************************************************************************************
  3146. '*4.2 Region 3. pSat_h and pSat_s
  3147. Function p3sat_h(ByVal h As Double) As Double
  3148. '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
  3149. '2004
  3150. 'Section 4 Boundary Equations psat(h) and psat(s) for the Saturation Lines of Region 3
  3151. 'Se pictures Page 17, Eq 10, Table 17, Page 18
  3152. Dim Ii, Ji, ni As Variant, ps As Double, i As Integer
  3153. Ii = Array(0, 1, 1, 1, 1, 5, 7, 8, 14, 20, 22, 24, 28, 36)
  3154. Ji = Array(0, 1, 3, 4, 36, 3, 0, 24, 16, 16, 3, 18, 8, 24)
  3155. ni = Array(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)
  3156. h = h / 2600
  3157. ps = 0
  3158. For i = 0 To 13
  3159.   ps = ps + ni(i) * (h - 1.02) ^ Ii(i) * (h - 0.608) ^ Ji(i)
  3160. Next i
  3161. p3sat_h = ps * 22
  3162. End Function
  3163. Function p3sat_s(ByVal s As Double) As Double
  3164. Dim Ii, Ji, ni As Variant, sigma, p As Double, i As Integer
  3165. Ii = Array(0, 1, 1, 4, 12, 12, 16, 24, 28, 32)
  3166. Ji = Array(0, 1, 32, 7, 4, 14, 36, 10, 0, 18)
  3167. ni = Array(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)
  3168. sigma = s / 5.2
  3169. p = 0
  3170. For i = 0 To 9
  3171.   p = p + ni(i) * (sigma - 1.03) ^ Ii(i) * (sigma - 0.699) ^ Ji(i)
  3172. Next i
  3173. p3sat_s = p * 22
  3174. End Function
  3175. '***********************************************************************************************************
  3176. '4.3 Region boundary 1to3 and 3to2 as a functions of s
  3177. Function hB13_s(ByVal s As Double) As Double
  3178. 'Supplementary Release on Backward Equations ( ) , p h s for Region 3,
  3179. 'Chapter 4.5 page 23.
  3180. Dim Ii, Ji, ni As Variant, sigma, eta As Double, i As Integer
  3181.   Ii = Array(0, 1, 1, 3, 5, 6)
  3182.   Ji = Array(0, -2, 2, -12, -4, -3)
  3183.   ni = Array(0.913965547600543, -4.30944856041991E-05, 60.3235694765419, 1.17518273082168E-18, 0.220000904781292, -69.0815545851641)
  3184.   sigma = s / 3.8
  3185.   eta = 0
  3186.   For i = 0 To 5
  3187.     eta = eta + ni(i) * (sigma - 0.884) ^ Ii(i) * (sigma - 0.864) ^ Ji(i)
  3188.   Next i
  3189.   hB13_s = eta * 1700
  3190. End Function
  3191. Function TB23_hs(ByVal h As Double, ByVal s As Double) As Double
  3192. 'Supplementary Release on Backward Equations ( ) , p h s for Region 3,
  3193. 'Chapter 4.6 page 25.
  3194. Dim Ii, Ji, ni As Variant, sigma, eta, teta As Double, i As Integer
  3195.    Ii = Array(-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)
  3196.    Ji = Array(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)
  3197.    ni = Array(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)
  3198.   sigma = s / 5.3
  3199.   eta = h / 3000
  3200.   teta = 0
  3201.   For i = 0 To 24
  3202.     teta = teta + ni(i) * (eta - 0.727) ^ Ii(i) * (sigma - 0.864) ^ Ji(i)
  3203.   Next i
  3204.   TB23_hs = teta * 900
  3205. End Function
  3206.  
  3207. '***********************************************************************************************************
  3208. '*5 Transport properties
  3209. '***********************************************************************************************************
  3210. '*5.1 Viscosity (IAPWS formulation 1985, Revised 2003)
  3211. '***********************************************************************************************************
  3212. Function my_AllRegions_pT(ByVal p As Double, ByVal t As Double) As Double
  3213. Dim h0, h1, h2, h3, h4, h5, h6 As Variant, rho, Ts, ps, my0, sum, my1, rhos As Double, i As Integer
  3214.   h0 = Array(0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447)
  3215.   h1 = Array(0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0)
  3216.   h2 = Array(-0.2818107, -1.070786, -1.263184, 0, 0, 0)
  3217.   h3 = Array(0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0)
  3218.   h4 = Array(-0.0417661, 0, 0, 0.1600435, 0, 0)
  3219.   h5 = Array(0, -0.01578386, 0, 0, 0, 0)
  3220.   h6 = Array(0, 0, 0, -0.003629481, 0, 0)
  3221.  
  3222.   'Calcualte density.
  3223. Select Case region_pT(p, t)
  3224.  Case 1
  3225.    rho = 1 / v1_pT(p, t)
  3226.  Case 2
  3227.    rho = 1 / v2_pT(p, t)
  3228.  Case 3
  3229.    rho = 1 / v3_ph(p, h3_pT(p, t))
  3230.  Case 4
  3231.    rho = CVErr(xlErrValue)
  3232.  Case 5
  3233.    rho = 1 / v5_pT(p, t)
  3234.  Case Else
  3235.   my_AllRegions_pT = CVErr(xlErrValue)
  3236.   Exit Function
  3237.  End Select
  3238.  
  3239.   rhos = rho / 317.763
  3240.   Ts = t / 647.226
  3241.   ps = p / 22.115
  3242.  
  3243.   'Check valid area
  3244.  If t > 900 + 273.15 Or (t > 600 + 273.15 And p > 300) Or (t > 150 + 273.15 And p > 350) Or p > 500 Then
  3245.     my_AllRegions_pT = CVErr(xlErrValue)
  3246.     Exit Function
  3247.   End If
  3248.   my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3))
  3249.   sum = 0
  3250.   For i = 0 To 5
  3251.       sum = sum + h0(i) * (1 / Ts - 1) ^ i + h1(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6
  3252.   Next i
  3253.   my1 = Exp(rhos * sum)
  3254.   my_AllRegions_pT = my0 * my1 * 0.000055071
  3255. End Function
  3256.  
  3257. Function my_AllRegions_ph(ByVal p As Double, ByVal h As Double) As Double
  3258.   Dim h0, h1, h2, h3, h4, h5, h6 As Variant, rho, t, Ts, ps, my0, sum, my1, rhos, v4V, v4L, xs As Double, i As Integer
  3259.   h0 = Array(0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447)
  3260.   h1 = Array(0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0)
  3261.   h2 = Array(-0.2818107, -1.070786, -1.263184, 0, 0, 0)
  3262.   h3 = Array(0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0)
  3263.   h4 = Array(-0.0417661, 0, 0, 0.1600435, 0, 0)
  3264.   h5 = Array(0, -0.01578386, 0, 0, 0, 0)
  3265.   h6 = Array(0, 0, 0, -0.003629481, 0, 0)
  3266.  
  3267.   'Calcualte density.
  3268. Select Case region_ph(p, h)
  3269.  Case 1
  3270.    Ts = T1_ph(p, h)
  3271.    t = Ts
  3272.    rho = 1 / v1_pT(p, Ts)
  3273.  Case 2
  3274.    Ts = T2_ph(p, h)
  3275.    t = Ts
  3276.    rho = 1 / v2_pT(p, Ts)
  3277.  Case 3
  3278.    rho = 1 / v3_ph(p, h)
  3279.    t = T3_ph(p, h)
  3280.  Case 4
  3281.    xs = x4_ph(p, h)
  3282.    If p < 16.529 Then
  3283.      v4V = v2_pT(p, T4_p(p))
  3284.      v4L = v1_pT(p, T4_p(p))
  3285.    Else
  3286.      v4V = v3_ph(p, h4V_p(p))
  3287.      v4L = v3_ph(p, h4L_p(p))
  3288.     End If
  3289.     rho = 1 / (xs * v4V + (1 - xs) * v4L)
  3290.     t = T4_p(p)
  3291.  Case 5
  3292.    Ts = T5_ph(p, h)
  3293.    t = Ts
  3294.    rho = 1 / v5_pT(p, Ts)
  3295.  Case Else
  3296.   my_AllRegions_ph = CVErr(xlErrValue)
  3297.   Exit Function
  3298.  End Select
  3299.   rhos = rho / 317.763
  3300.   Ts = t / 647.226
  3301.   ps = p / 22.115
  3302.   'Check valid area
  3303.  If t > 900 + 273.15 Or (t > 600 + 273.15 And p > 300) Or (t > 150 + 273.15 And p > 350) Or p > 500 Then
  3304.     my_AllRegions_ph = CVErr(xlErrValue)
  3305.     Exit Function
  3306.   End If
  3307.   my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3))
  3308.  
  3309.   sum = 0
  3310.   For i = 0 To 5
  3311.       sum = sum + h0(i) * (1 / Ts - 1) ^ i + h1(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6
  3312.   Next i
  3313.   my1 = Exp(rhos * sum)
  3314.   my_AllRegions_ph = my0 * my1 * 0.000055071
  3315. End Function
  3316. '***********************************************************************************************************
  3317. '*5.2 Thermal Conductivity (IAPWS formulation 1985)
  3318. Function tc_ptrho(ByVal p As Double, ByVal t As Double, ByVal rho As Double) As Double
  3319. 'Revised release on the IAPS Formulation 1985 for the Thermal Conductivity of ordinary water
  3320. 'IAPWS September 1998
  3321. 'Page 8
  3322. 'ver2.6 Start corrected bug
  3323. Dim tc0, tc1, dt, Q, s, tc2, tc As Double
  3324.  If t < 273.15 Then
  3325.    tc_ptrho = CVErr(xlErrValue) 'Out of range of validity (para. B4)
  3326.   Exit Function
  3327.  ElseIf t < 500 + 273.15 Then
  3328.    If p > 100 Then
  3329.      tc_ptrho = CVErr(xlErrValue) 'Out of range of validity (para. B4)
  3330.     Exit Function
  3331.    End If
  3332.  ElseIf t <= 650 + 273.15 Then
  3333.    If p > 70 Then
  3334.    tc_ptrho = CVErr(xlErrValue) 'Out of range of validity (para. B4)
  3335.   Exit Function
  3336.  End If
  3337.  ElseIf t <= 800 + 273.15 Then
  3338.    If p > 40 Then
  3339.       tc_ptrho = CVErr(xlErrValue) 'Out of range of validity (para. B4)
  3340.      Exit Function
  3341.    End If
  3342.  End If
  3343.  'ver2.6 End corrected bug
  3344.  
  3345.   t = t / 647.26
  3346.   rho = rho / 317.7
  3347.   tc0 = t ^ 0.5 * (0.0102811 + 0.0299621 * t + 0.0156146 * t ^ 2 - 0.00422464 * t ^ 3)
  3348.   tc1 = -0.39707 + 0.400302 * rho + 1.06 * Exp(-0.171587 * (rho + 2.39219) ^ 2)
  3349.   dt = Abs(t - 1) + 0.00308976
  3350.   Q = 2 + 0.0822994 / dt ^ (3 / 5)
  3351.   If t >= 1 Then
  3352.    s = 1 / dt
  3353.   Else
  3354.    s = 10.0932 / dt ^ (3 / 5)
  3355.   End If
  3356.   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)
  3357.   tc_ptrho = tc0 + tc1 + tc2
  3358. End Function
  3359. '***********************************************************************************************************
  3360. '5.3 Surface Tension
  3361. Function Surface_Tension_T(ByVal t As Double)
  3362. 'IAPWS Release on Surface Tension of Ordinary Water Substance,
  3363. 'September 1994
  3364. Dim tau As Double
  3365. Const tc As Double = 647.096, b As Double = 0.2358, bb As Double = -0.625, my As Double = 1.256
  3366. If t < 0.01 Or t > tc Then
  3367.  Surface_Tension_T = "Out of valid region"
  3368.  Exit Function
  3369. End If
  3370. tau = 1 - t / tc
  3371. Surface_Tension_T = b * tau ^ my * (1 + bb * tau)
  3372. End Function
  3373. '***********************************************************************************************************
  3374. '*6 Units                                                                                      *
  3375. '***********************************************************************************************************
  3376.  
  3377. Function toSIunit_p(ByVal Ins As Double) As Double
  3378. 'Translate bar to MPa
  3379.  toSIunit_p = Ins / 10
  3380. End Function
  3381. Function fromSIunit_p(ByVal Ins As Double) As Double
  3382. 'Translate bar to MPa
  3383.  fromSIunit_p = Ins * 10
  3384. End Function
  3385. Function toSIunit_T(ByVal Ins As Double) As Double
  3386. 'Translate degC to Kelvon
  3387.  toSIunit_T = Ins + 273.15
  3388. End Function
  3389. Function fromSIunit_T(ByVal Ins As Double) As Double
  3390. 'Translate Kelvin to degC
  3391.  fromSIunit_T = Ins - 273.15
  3392. End Function
  3393. Function toSIunit_h(ByVal Ins As Double) As Double
  3394.   toSIunit_h = Ins
  3395. End Function
  3396. Function fromSIunit_h(ByVal Ins As Double) As Double
  3397.   fromSIunit_h = Ins
  3398. End Function
  3399. Function toSIunit_v(ByVal Ins As Double) As Double
  3400.   toSIunit_v = Ins
  3401. End Function
  3402. Function fromSIunit_v(ByVal Ins As Double) As Double
  3403.   fromSIunit_v = Ins
  3404. End Function
  3405. Function toSIunit_s(ByVal Ins As Double) As Double
  3406.   toSIunit_s = Ins
  3407. End Function
  3408. Function fromSIunit_s(ByVal Ins As Double) As Double
  3409.   fromSIunit_s = Ins
  3410. End Function
  3411. Function toSIunit_u(ByVal Ins As Double) As Double
  3412.   toSIunit_u = Ins
  3413. End Function
  3414. Function fromSIunit_u(ByVal Ins As Double) As Double
  3415.   fromSIunit_u = Ins
  3416. End Function
  3417. Function toSIunit_Cp(ByVal Ins As Double) As Double
  3418.   toSIunit_Cp = Ins
  3419. End Function
  3420. Function fromSIunit_Cp(ByVal Ins As Double) As Double
  3421.   fromSIunit_Cp = Ins
  3422. End Function
  3423. Function toSIunit_Cv(ByVal Ins As Double) As Double
  3424.   toSIunit_Cv = Ins
  3425. End Function
  3426. Function fromSIunit_Cv(ByVal Ins As Double) As Double
  3427.   fromSIunit_Cv = Ins
  3428. End Function
  3429. Function toSIunit_w(ByVal Ins As Double) As Double
  3430.   toSIunit_w = Ins
  3431. End Function
  3432. Function fromSIunit_w(ByVal Ins As Double) As Double
  3433.   fromSIunit_w = Ins
  3434. End Function
  3435. Function toSIunit_tc(ByVal Ins As Double) As Double
  3436.   toSIunit_tc = Ins
  3437. End Function
  3438. Function fromSIunit_tc(ByVal Ins As Double) As Double
  3439.   fromSIunit_tc = Ins
  3440. End Function
  3441. Function toSIunit_st(ByVal Ins As Double) As Double
  3442.   toSIunit_st = Ins
  3443. End Function
  3444. Function fromSIunit_st(ByVal Ins As Double) As Double
  3445.   fromSIunit_st = Ins
  3446. End Function
  3447. Function toSIunit_x(ByVal Ins As Double) As Double
  3448.   toSIunit_x = Ins
  3449. End Function
  3450. Function fromSIunit_x(ByVal Ins As Double) As Double
  3451.   fromSIunit_x = Ins
  3452. End Function
  3453. Function toSIunit_vx(ByVal Ins As Double) As Double
  3454.   toSIunit_vx = Ins
  3455. End Function
  3456. Function fromSIunit_vx(ByVal Ins As Double) As Double
  3457.   fromSIunit_vx = Ins
  3458. End Function
  3459. Function toSIunit_my(ByVal Ins As Double) As Double
  3460.   toSIunit_my = Ins
  3461. End Function
  3462. Function fromSIunit_my(ByVal Ins As Double) As Double
  3463.   fromSIunit_my = Ins
  3464. End Function
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top