https://mooseframework.inl.gov
Water97FluidProperties.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #pragma once
11 
13 #include "NewtonInversion.h"
14 #include <array>
15 
16 #pragma GCC diagnostic push
17 #pragma GCC diagnostic ignored "-Woverloaded-virtual"
18 
39 {
40 public:
42 
44  virtual ~Water97FluidProperties();
45 
46  virtual std::string fluidName() const override;
47 
48  virtual Real molarMass() const override;
49 
50  virtual Real criticalPressure() const override;
51 
52  virtual Real criticalTemperature() const override;
53 
54  virtual Real criticalDensity() const override;
55 
56  virtual Real triplePointPressure() const override;
57 
58  virtual Real triplePointTemperature() const override;
59 
60  virtual Real rho_from_p_T(Real pressure, Real temperature) const override;
61  virtual ADReal rho_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
62 
63  virtual Real v_from_p_T(Real pressure, Real temperature) const override;
64  virtual ADReal v_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
65  template <typename T>
66  T v_from_p_T_template(const T & pressure, const T & temperature) const;
67  virtual void
68  v_from_p_T(Real pressure, Real temperature, Real & v, Real & dv_dp, Real & dv_dT) const override;
69  virtual void v_from_p_T(const ADReal & pressure,
70  const ADReal & temperature,
71  ADReal & v,
72  ADReal & dv_dp,
73  ADReal & dv_dT) const override;
74  template <typename T>
75  void
76  v_from_p_T_template(const T & pressure, const T & temperature, T & v, T & dv_dp, T & dv_dt) const;
77 
78  template <typename T>
79  T rho_from_p_T_template(const T & pressure, const T & temperature) const;
80  template <typename T>
82  const T & pressure, const T & temperature, T & rho, T & drho_dp, T & drho_dt) const;
83 
84  virtual void rho_from_p_T(
85  Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const override;
86  virtual void rho_from_p_T(const ADReal & pressure,
87  const ADReal & temperature,
88  ADReal & rho,
89  ADReal & drho_dp,
90  ADReal & drho_dT) const override;
91 
92  Real p_from_v_e(Real v, Real e) const override;
93  ADReal p_from_v_e(const ADReal & v, const ADReal & e) const override;
94  template <typename T>
95  T p_from_v_e_template(const T & v, const T & e) const;
96 
97  virtual Real e_from_p_rho(Real p, Real rho) const override;
98  virtual ADReal e_from_p_rho(const ADReal & p, const ADReal & rho) const override;
99  template <typename T>
100  T e_from_p_rho_template(const T & p, const T & rho) const;
101  void e_from_p_rho(Real p, Real rho, Real & e, Real & de_dp, Real & de_drho) const override;
102  void e_from_p_rho(const ADReal & p,
103  const ADReal & rho,
104  ADReal & e,
105  ADReal & de_dp,
106  ADReal & de_drho) const override;
107  template <typename T>
108  void e_from_p_rho_template(const T & p, const T & rho, T & e, T & de_dp, T & de_drho) const;
109 
110  virtual Real e_from_p_T(Real pressure, Real temperature) const override;
111  virtual ADReal e_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
112  template <typename T>
113  T e_from_p_T_template(const T & pressure, const T & temperature) const;
114 
115  virtual void
116  e_from_p_T(Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const override;
117  virtual void e_from_p_T(const ADReal & pressure,
118  const ADReal & temperature,
119  ADReal & e,
120  ADReal & de_dp,
121  ADReal & de_dT) const override;
122  template <typename T>
123  void
124  e_from_p_T_template(const T & pressure, const T & temperature, T & e, T & de_dp, T & de_dT) const;
125 
126  ADReal e_from_v_h(const ADReal & v, const ADReal & h) const override;
127 
128  Real T_from_v_e(Real v, Real e) const override;
129  ADReal T_from_v_e(const ADReal & v, const ADReal & e) const override;
130 
131  virtual Real c_from_p_T(Real pressure, Real temperature) const override;
132  virtual ADReal c_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
133  template <typename T>
134  T c_from_p_T_template(const T & pressure, const T & temperature) const;
135  virtual ADReal c_from_v_e(const ADReal & v, const ADReal & e) const override;
136 
137  virtual Real cp_from_p_T(Real pressure, Real temperature) const override;
138  virtual ADReal cp_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
139  template <typename T>
140  T cp_from_p_T_template(const T & pressure, const T & temperature) const;
141 
142  using SinglePhaseFluidProperties::cp_from_p_T;
143 
144  virtual ADReal cp_from_v_e(const ADReal & v, const ADReal & e) const override;
145 
146  virtual Real cv_from_p_T(Real pressure, Real temperature) const override;
147  virtual ADReal cv_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
148  template <typename T>
149  T cv_from_p_T_template(const T & pressure, const T & temperature) const;
150 
151  virtual ADReal cv_from_v_e(const ADReal & v, const ADReal & e) const override;
152 
153  virtual Real mu_from_p_T(Real pressure, Real temperature) const override;
154  virtual ADReal mu_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
155  template <typename T>
156  T mu_from_p_T_template(const T & pressure, const T & temperature) const;
157 
158  virtual void mu_from_p_T(
159  Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const override;
160 
161  virtual Real mu_from_rho_T(Real density, Real temperature) const override;
162  template <typename T>
163  T mu_from_rho_T_template(const T & density, const T & temperature) const;
164 
165  void mu_from_rho_T(
166  Real rho, Real temperature, Real drho_dT, Real & mu, Real & dmu_drho, Real & dmu_dT) const;
167 
168  ADReal mu_from_v_e(const ADReal & v, const ADReal & e) const override;
169 
170  virtual void
171  rho_mu_from_p_T(Real pressure, Real temperature, Real & rho, Real & mu) const override;
172 
173  virtual void rho_mu_from_p_T(Real pressure,
175  Real & rho,
176  Real & drho_dp,
177  Real & drho_dT,
178  Real & mu,
179  Real & dmu_dp,
180  Real & dmu_dT) const override;
181 
182  virtual Real k_from_p_T(Real pressure, Real temperature) const override;
183  virtual ADReal k_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
184  template <typename T>
185  T k_from_p_T_template(const T & pressure, const T & temperature) const;
186 
187  virtual void
188  k_from_p_T(Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const override;
189 
190  virtual Real k_from_rho_T(Real density, Real temperature) const override;
191  template <typename T>
192  T k_from_rho_T_template(const T & density, const T & temperature) const;
193 
194  Real k_from_v_e(Real v, Real e) const override;
195  ADReal k_from_v_e(const ADReal & v, const ADReal & e) const override;
196  template <typename T>
197  T k_from_v_e_template(const T & v, const T & e) const;
198 
200 
201  virtual Real h_from_p_T(Real pressure, Real temperature) const override;
202  virtual ADReal h_from_p_T(const ADReal & pressure, const ADReal & temperature) const override;
203  template <typename T>
204  T h_from_p_T_template(const T & pressure, const T & temperature) const;
205 
206  virtual void
207  h_from_p_T(Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const override;
208  virtual void h_from_p_T(const ADReal & pressure,
209  const ADReal & temperature,
210  ADReal & h,
211  ADReal & dh_dp,
212  ADReal & dh_dT) const override;
213  template <typename T>
214  void
215  h_from_p_T_template(const T & pressure, const T & temperature, T & h, T & dh_dp, T & dh_dT) const;
216 
217  virtual Real s_from_h_p(Real enthalpy, Real pressure) const override;
218  virtual ADReal s_from_h_p(const ADReal & enthalpy, const ADReal & pressure) const override;
219  void virtual s_from_h_p(
220  Real enthalpy, Real pressure, Real & s, Real & ds_dh, Real & ds_dp) const override;
221 
222  virtual Real vaporPressure(Real temperature) const override;
223 
224  virtual void vaporPressure(Real temperature, Real & psat, Real & dpsat_dT) const override;
225 
226  template <typename T>
227  void vaporPressureTemplate(const T & temperature, T & psat, T & dpsat_dT) const;
228 
241  Real vaporTemperature(Real pressure) const override;
242  virtual void vaporTemperature(Real pressure, Real & Tsat, Real & dTsat_dp) const override;
243 
254  Real b23p(Real temperature) const;
255 
266  Real b23T(Real pressure) const;
267 
276  unsigned int inRegion(Real pressure, Real temperature) const;
277 
289  unsigned int subregion3(Real pressure, Real temperature) const;
290 
300  template <typename T>
301  T subregionVolume(const T & pi,
302  const T & theta,
303  Real a,
304  Real b,
305  Real c,
306  Real d,
307  Real e,
308  unsigned int sid) const;
309 
323  template <typename T>
324  T densityRegion3(const T & pressure, const T & temperature) const;
325 
335  virtual Real T_from_p_h(Real pressure, Real enthalpy) const override;
336  virtual void T_from_p_h(Real p, Real h, Real & T, Real & dT_dp, Real & dT_dh) const override;
337  virtual ADReal T_from_p_h(const ADReal & pressure, const ADReal & enthalpy) const override;
338 
348  Real b2bc(Real pressure) const;
349 
360  Real b3ab(Real pressure) const;
361 
372  Real henryConstant(Real temperature, const std::vector<Real> & coeffs) const;
373  void
374  henryConstant(Real temperature, const std::vector<Real> & coeffs, Real & Kh, Real & dKh_dT) const;
375  ADReal henryConstant(const ADReal & temperature, const std::vector<Real> & coeffs) const;
376 
381  template <typename T>
382  std::pair<T, T> p_T_from_v_e(const T & v, const T & e) const;
383 
388  template <typename T>
389  std::pair<T, T> rho_T_from_v_e(const T & v, const T & e) const;
390 
395  template <typename T>
396  std::pair<T, T> p_T_from_v_h(const T & v, const T & h) const;
398 
399 protected:
411  template <typename T>
412  T gamma1(const T & pi, const T & tau) const;
413 
421  template <typename T>
422  T dgamma1_dpi(const T & pi, const T & tau) const;
423 
431  template <typename T>
432  T d2gamma1_dpi2(const T & pi, const T & tau) const;
433 
441  template <typename T>
442  T dgamma1_dtau(const T & pi, const T & tau) const;
443 
451  template <typename T>
452  T d2gamma1_dtau2(const T & pi, const T & tau) const;
453 
461  template <typename T>
462  T d2gamma1_dpitau(const T & pi, const T & tau) const;
463 
475  template <typename T>
476  T gamma2(const T & pi, const T & tau) const;
477 
485  template <typename T>
486  T dgamma2_dpi(const T & pi, const T & tau) const;
487 
495  template <typename T>
496  T d2gamma2_dpi2(const T & pi, const T & tau) const;
497 
505  template <typename T>
506  T dgamma2_dtau(const T & pi, const T & tau) const;
507 
515  template <typename T>
516  T d2gamma2_dtau2(const T & pi, const T & tau) const;
517 
525  template <typename T>
526  T d2gamma2_dpitau(const T & pi, const T & tau) const;
527 
539  template <typename T>
540  T phi3(const T & delta, const T & tau) const;
541 
549  template <typename T>
550  T dphi3_ddelta(const T & delta, const T & tau) const;
551 
559  template <typename T>
560  T d2phi3_ddelta2(const T & delta, const T & tau) const;
561 
569  template <typename T>
570  T dphi3_dtau(const T & delta, const T & tau) const;
571 
579  template <typename T>
580  T d2phi3_dtau2(const T & delta, const T & tau) const;
581 
589  template <typename T>
590  T d2phi3_ddeltatau(const T & delta, const T & tau) const;
591 
603  template <typename T>
604  T gamma5(const T & pi, const T & tau) const;
605 
613  template <typename T>
614  T dgamma5_dpi(const T & pi, const T & tau) const;
615 
623  template <typename T>
624  T d2gamma5_dpi2(const T & pi, const T & tau) const;
625 
633  template <typename T>
634  T dgamma5_dtau(const T & pi, const T & tau) const;
635 
643  template <typename T>
644  T d2gamma5_dtau2(const T & pi, const T & tau) const;
645 
653  template <typename T>
654  T d2gamma5_dpitau(const T & pi, const T & tau) const;
655 
658  {
659  AB,
660  CD,
661  GH,
662  IJ,
663  JK,
664  MN,
665  OP,
666  QU,
667  RX,
668  UV,
669  WX,
671  };
672 
684  template <typename T>
685  T tempXY(const T & pressure, subregionEnum xy) const;
686 
695  unsigned int inRegionPH(Real pressure, Real enthalpy) const;
696 
707  unsigned int subregion2ph(Real pressure, Real enthalpy) const;
708 
720  unsigned int subregion3ph(Real pressure, Real enthalpy) const;
721 
731  ADReal T_from_p_h_ad(const ADReal & pressure, const ADReal & enthalpy) const;
732 
743  ADReal temperature_from_ph1(const ADReal & pressure, const ADReal & enthalpy) const;
744 
755  ADReal temperature_from_ph2a(const ADReal & pressure, const ADReal & enthalpy) const;
756 
767  ADReal temperature_from_ph2b(const ADReal & pressure, const ADReal & enthalpy) const;
768 
779  ADReal temperature_from_ph2c(const ADReal & pressure, const ADReal & enthalpy) const;
780 
792  ADReal temperature_from_ph3a(const ADReal & pressure, const ADReal & enthalpy) const;
793 
805  ADReal temperature_from_ph3b(const ADReal & pressure, const ADReal & enthalpy) const;
806 
819  ADReal vaporTemperature_ad(const ADReal & pressure) const;
820 
822  const Real _Mh2o;
824  const Real _Rw;
835 
844  const std::array<Real, 34> _n1{
846  {0.14632971213167e0, -0.84548187169114e0, -0.37563603672040e1, 0.33855169168385e1,
847  -0.95791963387872e0, 0.15772038513228e0, -0.16616417199501e-1, 0.81214629983568e-3,
848  0.28319080123804e-3, -0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1,
849  -0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3, -0.30001780793026e-3,
850  0.47661393906987e-4, -0.44141845330846e-5, -0.72694996297594e-15, -0.31679644845054e-4,
851  -0.28270797985312e-5, -0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6,
852  -0.14341729937924e-12, -0.40516996860117e-6, -0.12734301741641e-8, -0.17424871230634e-9,
853  -0.68762131295531e-18, 0.14478307828521e-19, 0.26335781662795e-22, -0.11947622640071e-22,
854  0.18228094581404e-23, -0.93537087292458e-25}};
855 
856  const std::array<int, 34> _I1{{0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2,
857  2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32}};
858 
859  const std::array<int, 34> _J1{{-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0,
860  1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2,
861  10, -8, -11, -6, -29, -31, -38, -39, -40, -41}};
862 
863  const std::array<Real, 20> _nph1{
864  {-0.23872489924521e3, 0.40421188637945e3, 0.11349746881718e3, -0.58457616048039e1,
865  -0.15285482413140e-3, -0.10866707695377e-5, -0.13391744872602e2, 0.43211039183559e2,
866  -0.54010067170506e2, 0.30535892203916e2, -0.65964749423638e1, 0.93965400878363e-2,
867  0.11573647505340e-6, -0.25858641282073e-4, -0.40644363084799e-8, 0.66456186191635e-7,
868  0.80670734103027e-10, -0.93477771213947e-12, 0.58265442020601e-14, -0.15020185953503e-16}};
869 
870  const std::array<int, 20> _Iph1{{0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6}};
871 
872  const std::array<int, 20> _Jph1{
873  {0, 1, 2, 6, 22, 32, 0, 1, 2, 3, 4, 10, 32, 10, 32, 10, 32, 32, 32, 32}};
874 
876  const std::array<Real, 9> _n02{{-0.96927686500217e1,
877  0.10086655968018e2,
878  -0.56087911283020e-2,
879  0.71452738081455e-1,
880  -0.40710498223928e0,
881  0.14240819171444e1,
882  -0.43839511319450e1,
883  -0.28408632460772e0,
884  0.21268463753307e-1}};
885 
886  const std::array<Real, 43> _n2{
887  {-0.17731742473213e-2, -0.17834862292358e-1, -0.45996013696365e-1, -0.57581259083432e-1,
888  -0.50325278727930e-1, -0.33032641670203e-4, -0.18948987516315e-3, -0.39392777243355e-2,
889  -0.43797295650573e-1, -0.26674547914087e-4, 0.20481737692309e-7, 0.43870667284435e-6,
890  -0.32277677238570e-4, -0.15033924542148e-2, -0.40668253562649e-1, -0.78847309559367e-9,
891  0.12790717852285e-7, 0.48225372718507e-6, 0.22922076337661e-5, -0.16714766451061e-10,
892  -0.21171472321355e-2, -0.23895741934104e2, -0.59059564324270e-17, -0.12621808899101e-5,
893  -0.38946842435739e-1, 0.11256211360459e-10, -0.82311340897998e1, 0.19809712802088e-7,
894  0.10406965210174e-18, -0.10234747095929e-12, -0.10018179379511e-8, -0.80882908646985e-10,
895  0.10693031879409e0, -0.33662250574171e0, 0.89185845355421e-24, 0.30629316876232e-12,
896  -0.42002467698208e-5, -0.59056029685639e-25, 0.37826947613457e-5, -0.12768608934681e-14,
897  0.73087610595061e-28, 0.55414715350778e-16, -0.94369707241210e-6}};
898 
899  const std::array<int, 9> _J02{{0, 1, -5, -4, -3, -2, -1, 2, 3}};
900 
901  const std::array<int, 43> _I2{{1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,
902  4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10,
903  10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24}};
904 
905  const std::array<int, 43> _J2{{0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35,
906  1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10,
907  14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58}};
908 
909  const std::array<Real, 36> _nph2a{
910  {0.10898952318288e4, 0.84951654495535e3, -0.10781748091826e3, 0.33153654801263e2,
911  -0.74232016790248e1, 0.11765048724356e2, 0.18445749355790e1, -0.41792700549624e1,
912  0.62478196935812e1, -0.17344563108114e2, -0.20058176862096e3, 0.27196065473796e3,
913  -0.45511318285818e3, 0.30919688604755e4, 0.25226640357872e6, -0.61707422868339e-2,
914  -0.31078046629583, 0.11670873077107e2, 0.12812798404046e9, -0.98554909623276e9,
915  0.28224546973002e10, -0.35948971410703e10, 0.17227349913197e10, -0.13551334240775e5,
916  0.12848734664650e8, 0.13865724283226e1, 0.23598832556514e6, -0.13105236545054e8,
917  0.73999835474766e4, -0.55196697030060e6, 0.37154085996233e7, 0.19127729239660e5,
918  -0.41535164835634e6, -0.62459855192507e2}};
919 
920  const std::array<int, 36> _Iph2a{{0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,
921  2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7}};
922 
923  const std::array<int, 36> _Jph2a{{0, 1, 2, 3, 7, 20, 0, 1, 2, 3, 7, 9,
924  11, 18, 44, 0, 2, 7, 36, 38, 40, 42, 44, 24,
925  44, 12, 32, 44, 32, 36, 42, 34, 44, 28}};
926 
927  const std::array<Real, 38> _nph2b{
928  {0.14895041079516e4, 0.74307798314034e3, -0.97708318797837e2, 0.24742464705674e1,
929  -0.63281320016026, 0.11385952129658e1, -0.47811863648625, 0.85208123431544e-2,
930  0.93747147377932, 0.33593118604916e1, 0.33809355601454e1, 0.16844539671904,
931  0.73875745236695, -0.47128737436186, 0.15020273139707, -0.21764114219750e-2,
932  -0.21810755324761e-1, -0.10829784403677, -0.46333324635812e-1, 0.71280351959551e-4,
933  0.11032831789999e-3, 0.18955248387902e-3, 0.30891541160537e-2, 0.13555504554949e-2,
934  0.28640237477456e-6, -0.10779857357512e-4, -0.76462712454814e-4, 0.14052392818316e-4,
935  -0.31083814331434e-4, -0.10302738212103e-5, 0.28217281635040e-6, 0.12704902271945e-5,
936  0.73803353468292e-7, -0.11030139238909e-7, -0.81456365207833e-13, -0.25180545682962e-10,
937  -0.17565233969407e-17, 0.86934156344163e-14}};
938 
939  const std::array<int, 38> _Iph2b{{0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2,
940  2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 7, 7, 9, 9}};
941 
942  const std::array<int, 38> _Jph2b{{0, 1, 2, 12, 18, 24, 28, 40, 0, 2, 6, 12, 18,
943  24, 28, 40, 2, 8, 18, 40, 1, 2, 12, 24, 2, 12,
944  18, 24, 28, 40, 18, 24, 40, 28, 2, 28, 1, 40}};
945 
946  const std::array<Real, 23> _nph2c{
947  {-0.32368398555242e13, 0.73263350902181e13, 0.35825089945447e12, -0.58340131851590e12,
948  -0.10783068217470e11, 0.20825544563171e11, 0.61074783564516e6, 0.85977722535580e6,
949  -0.25745723604170e5, 0.31081088422714e5, 0.12082315865936e4, 0.48219755109255e3,
950  0.37966001272486e1, -0.10842984880077e2, -0.45364172676660e-1, 0.14559115658698e-12,
951  0.11261597407230e-11, -0.17804982240686e-10, 0.12324579690832e-6, -0.11606921130984e-5,
952  0.27846367088554e-4, -0.59270038474176e-3, 0.12918582991878e-2}};
953 
954  const std::array<int, 23> _Iph2c{
955  {-7, -7, -6, -6, -5, -5, -2, -2, -1, -1, 0, 0, 1, 1, 2, 6, 6, 6, 6, 6, 6, 6, 6}};
956 
957  const std::array<int, 23> _Jph2c{
958  {0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16, 20, 22}};
959 
961  const std::array<Real, 5> _n23{{0.34805185628969e3,
962  -0.11671859879975e1,
963  0.10192970039326e-2,
964  0.57254459862746e3,
965  0.13918839778870e2}};
966 
968  const std::array<Real, 40> _n3{
969  {0.10658070028513e1, -0.15732845290239e2, 0.20944396974307e2, -0.76867707878716e1,
970  0.26185947787954e1, -0.28080781148620e1, 0.12053369696517e1, -0.84566812812502e-2,
971  -0.12654315477714e1, -0.11524407806681e1, 0.88521043984318e0, -0.64207765181607e0,
972  0.38493460186671e0, -0.85214708824206e0, 0.48972281541877e1, -0.30502617256965e1,
973  0.39420536879154e-1, 0.12558408424308e0, -0.27999329698710e0, 0.13899799569460e1,
974  -0.20189915023570e1, -0.82147637173963e-2, -0.47596035734923e0, 0.43984074473500e-1,
975  -0.44476435428739e0, 0.90572070719733e0, 0.70522450087967e0, 0.10770512626332e0,
976  -0.32913623258954e0, -0.50871062041158e0, -0.22175400873096e-1, 0.94260751665092e-1,
977  0.16436278447961e0, -0.13503372241348e-1, -0.14834345352472e-1, 0.57922953628084e-3,
978  0.32308904703711e-2, 0.80964802996215e-4, -0.16557679795037e-3, -0.44923899061815e-4}};
979 
980  const std::array<int, 40> _I3{{0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3,
981  3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11}};
982 
983  const std::array<int, 40> _J3{{0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2,
984  6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1,
985  3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26}};
986 
987  const std::array<std::array<Real, 8>, 26> _par3{
988  {{{0.0024, 100.0, 760.0, 0.085, 0.817, 1.0, 1.0, 1.0}},
989  {{0.0041, 100.0, 860.0, 0.280, 0.779, 1.0, 1.0, 1.0}},
990  {{0.0022, 40.0, 690.0, 0.259, 0.903, 1.0, 1.0, 1.0}},
991  {{0.0029, 40.0, 690.0, 0.559, 0.939, 1.0, 1.0, 4.0}},
992  {{0.0032, 40.0, 710.0, 0.587, 0.918, 1.0, 1.0, 1.0}},
993  {{0.0064, 40.0, 730.0, 0.587, 0.891, 0.5, 1.0, 4.0}},
994  {{0.0027, 25.0, 660.0, 0.872, 0.971, 1.0, 1.0, 4.0}},
995  {{0.0032, 25.0, 660.0, 0.898, 0.983, 1.0, 1.0, 4.0}},
996  {{0.0041, 25.0, 660.0, 0.910, 0.984, 0.5, 1.0, 4.0}},
997  {{0.0054, 25.0, 670.0, 0.875, 0.964, 0.5, 1.0, 4.0}},
998  {{0.0077, 25.0, 680.0, 0.802, 0.935, 1.0, 1.0, 1.0}},
999  {{0.0026, 24.0, 650.0, 0.908, 0.989, 1.0, 1.0, 4.0}},
1000  {{0.0028, 23.0, 650.0, 1.000, 0.997, 1.0, 0.25, 1.0}},
1001  {{0.0031, 23.0, 650.0, 0.976, 0.997, 0.0, 0.0, 0.0}},
1002  {{0.0034, 23.0, 650.0, 0.974, 0.996, 0.5, 1.0, 1.0}},
1003  {{0.0041, 23.0, 650.0, 0.972, 0.997, 0.5, 1.0, 1.0}},
1004  {{0.0022, 23.0, 650.0, 0.848, 0.983, 1.0, 1.0, 4.0}},
1005  {{0.0054, 23.0, 650.0, 0.874, 0.982, 1.0, 1.0, 1.0}},
1006  {{0.0022, 21.0, 640.0, 0.886, 0.990, 1.0, 1.0, 4.0}},
1007  {{0.0088, 20.0, 650.0, 0.803, 1.020, 1.0, 1.0, 1.0}},
1008  {{0.0026, 23.0, 650.0, 0.902, 0.988, 1.0, 1.0, 1.0}},
1009  {{0.0031, 23.0, 650.0, 0.960, 0.995, 1.0, 1.0, 1.0}},
1010  {{0.0039, 23.0, 650.0, 0.959, 0.995, 1.0, 1.0, 4.0}},
1011  {{0.0049, 23.0, 650.0, 0.910, 0.988, 1.0, 1.0, 1.0}},
1012  {{0.0031, 22.0, 650.0, 0.996, 0.994, 1.0, 1.0, 4.0}},
1013  {{0.0038, 22.0, 650.0, 0.993, 0.994, 1.0, 1.0, 4.0}}}};
1014 
1015  const std::array<unsigned int, 26> _par3N{{30, 32, 35, 38, 29, 42, 38, 29, 42, 29, 34, 43, 40,
1016  39, 24, 27, 24, 27, 29, 33, 38, 39, 35, 36, 20, 23}};
1017 
1019  const std::vector<std::vector<Real>> _n3s{
1020  {0.110879558823853e-2, 0.572616740810616e3, -0.767051948380852e5, -0.253321069529674e-1,
1021  0.628008049345689e4, 0.234105654131876e6, 0.216867826045856, -0.156237904341963e3,
1022  -0.269893956176613e5, -0.180407100085505e-3, 0.116732227668261e-2, 0.266987040856040e2,
1023  0.282776617243286e5, -0.242431520029523e4, 0.435217323022733e-3, -0.122494831387441e-1,
1024  0.179357604019989e1, 0.442729521058314e2, -0.593223489018342e-2, 0.453186261685774,
1025  0.135825703129140e1, 0.408748415856745e-1, 0.474686397863312, 0.118646814997915e1,
1026  0.546987265727549, 0.195266770452643, -0.502268790869663e-1, -0.369645308193377,
1027  0.633828037528420e-2, 0.797441793901017e-1},
1028  {-0.827670470003621e-1, 0.416887126010565e2, 0.483651982197059e-1, -0.291032084950276e5,
1029  -0.111422582236948e3, -0.202300083904014e-1, 0.294002509338515e3, 0.140244997609658e3,
1030  -0.344384158811459e3, 0.361182452612149e3, -0.140699677420738e4, -0.202023902676481e-2,
1031  0.171346792457471e3, -0.425597804058632e1, 0.691346085000334e-5, 0.151140509678925e-2,
1032  -0.416375290166236e-1, -0.413754957011042e2, -0.506673295721637e2, -0.572212965569023e-3,
1033  0.608817368401785e1, 0.239600660256161e2, 0.122261479925384e-1, 0.216356057692938e1,
1034  0.398198903368642, -0.116892827834085, -0.102845919373532, -0.492676637589284,
1035  0.655540456406790e-1, -0.240462535078530, -0.269798180310075e-1, 0.128369435967012},
1036  {0.311967788763030e1, 0.276713458847564e5, 0.322583103403269e8, -0.342416065095363e3,
1037  -0.899732529907377e6, -0.793892049821251e8, 0.953193003217388e2, 0.229784742345072e4,
1038  0.175336675322499e6, 0.791214365222792e7, 0.319933345844209e-4, -0.659508863555767e2,
1039  -0.833426563212851e6, 0.645734680583292e-1, -0.382031020570813e7, 0.406398848470079e-4,
1040  0.310327498492008e2, -0.892996718483724e-3, 0.234604891591616e3, 0.377515668966951e4,
1041  0.158646812591361e-1, 0.707906336241843, 0.126016225146570e2, 0.736143655772152,
1042  0.676544268999101, -0.178100588189137e2, -0.156531975531713, 0.117707430048158e2,
1043  0.840143653860447e-1, -0.186442467471949, -0.440170203949645e2, 0.123290423502494e7,
1044  -0.240650039730845e-1, -0.107077716660869e7, 0.438319858566475e-1},
1045  {-0.452484847171645e-9, 0.315210389538801e-4, -0.214991352047545e-2, 0.508058874808345e3,
1046  -0.127123036845932e8, 0.115371133120497e13, -0.197805728776273e-15, 0.241554806033972e-10,
1047  -0.156481703640525e-5, 0.277211346836625e-2, -0.203578994462286e2, 0.144369489909053e7,
1048  -0.411254217946539e11, 0.623449786243773e-5, -0.221774281146038e2, -0.689315087933158e5,
1049  -0.195419525060713e8, 0.316373510564015e4, 0.224040754426988e7, -0.436701347922356e-5,
1050  -0.404213852833996e-3, -0.348153203414663e3, -0.385294213555289e6, 0.135203700099403e-6,
1051  0.134648383271089e-3, 0.125031835351736e6, 0.968123678455841e-1, 0.225660517512438e3,
1052  -0.190102435341872e-3, -0.299628410819229e-1, 0.500833915372121e-2, 0.387842482998411,
1053  -0.138535367777182e4, 0.870745245971773, 0.171946252068742e1, -0.326650121426383e-1,
1054  0.498044171727877e4, 0.551478022765087e-2},
1055  {0.715815808404721e9, -0.114328360753449e12, 0.376531002015720e-11, -0.903983668691157e-4,
1056  0.665695908836252e6, 0.535364174960127e10, 0.794977402335603e11, 0.922230563421437e2,
1057  -0.142586073991215e6, -0.111796381424162e7, 0.896121629640760e4, -0.669989239070491e4,
1058  0.451242538486834e-2, -0.339731325977713e2, -0.120523111552278e1, 0.475992667717124e5,
1059  -0.266627750390341e6, -0.153314954386524e-3, 0.305638404828265, 0.123654999499486e3,
1060  -0.104390794213011e4, -0.157496516174308e-1, 0.685331118940253, 0.178373462873903e1,
1061  -0.544674124878910, 0.204529931318843e4, -0.228342359328752e5, 0.413197481515899,
1062  -0.341931835910405e2},
1063  {-0.251756547792325e-7, 0.601307193668763e-5, -0.100615977450049e-2,
1064  0.999969140252192, 0.214107759236486e1, -0.165175571959086e2,
1065  -0.141987303638727e-2, 0.269251915156554e1, 0.349741815858722e2,
1066  -0.300208695771783e2, -0.131546288252539e1, -0.839091277286169e1,
1067  0.181545608337015e-9, -0.591099206478909e-3, 0.152115067087106e1,
1068  0.252956470663225e-4, 0.100726265203786e-14, -0.14977453386065e1,
1069  -0.793940970562969e-9, -0.150290891264717e-3, 0.151205531275133e1,
1070  0.470942606221652e-5, 0.195049710391712e-12, -0.911627886266077e-8,
1071  0.604374640201265e-3, -0.225132933900136e-15, 0.610916973582981e-11,
1072  -0.303063908043404e-6, -0.137796070798409e-4, -0.919296736666106e-3,
1073  0.639288223132545e-9, 0.753259479898699e-6, -0.400321478682929e-12,
1074  0.756140294351614e-8, -0.912082054034891e-11, -0.237612381140539e-7,
1075  0.269586010591874e-4, -0.732828135157839e-10, 0.241995578306660e-9,
1076  -0.405735532730322e-3, 0.189424143498011e-9, -0.486632965074563e-9},
1077  {0.412209020652996e-4, -0.114987238280587e7, 0.948180885032080e10, -0.195788865718971e18,
1078  0.4962507048713e25, -0.105549884548496e29, -0.758642165988278e12, -0.922172769596101e23,
1079  0.725379072059348e30, -0.617718249205859e2, 0.107555033344858e5, -0.379545802336487e8,
1080  0.228646846221831e12, -0.499741093010619e7, -0.280214310054101e31, 0.104915406769586e7,
1081  0.613754229168619e28, 0.802056715528378e32, -0.298617819828065e8, -0.910782540134681e2,
1082  0.135033227281565e6, -0.712949383408211e19, -0.104578785289542e37, 0.304331584444093e2,
1083  0.593250797959445e10, -0.364174062110798e28, 0.921791403532461, -0.337693609657471,
1084  -0.724644143758508e2, -0.110480239272601, 0.536516031875059e1, -0.291441872156205e4,
1085  0.616338176535305e40, -0.120889175861180e39, 0.818396024524612e23, 0.940781944835829e9,
1086  -0.367279669545448e5, -0.837513931798655e16},
1087  {0.561379678887577e-1, 0.774135421587083e10, 0.111482975877938e-8, -0.143987128208183e-2,
1088  0.193696558764920e4, -0.605971823585005e9, 0.171951568124337e14, -0.185461154985145e17,
1089  0.38785116807801e-16, -0.395464327846105e-13, -0.170875935679023e3, -0.21201062070122e4,
1090  0.177683337348191e8, 0.110177443629575e2, -0.234396091693313e6, -0.656174421999594e7,
1091  0.156362212977396e-4, -0.212946257021400e1, 0.135249306374858e2, 0.177189164145813,
1092  0.139499167345464e4, -0.703670932036388e-2, -0.152011044389648, 0.981916922991113e-4,
1093  0.147199658618076e-2, 0.202618487025578e2, 0.899345518944240, -0.211346402240858,
1094  0.249971752957491e2},
1095  {0.106905684359136e1, -0.148620857922333e1, 0.259862256980408e15,
1096  -0.446352055678749e-11, -0.566620757170032e-6, -0.235302885736849e-2,
1097  -0.269226321968839, 0.922024992944392e1, 0.357633505503772e-11,
1098  -0.173942565562222e2, 0.700681785556229e-5, -0.267050351075768e-3,
1099  -0.231779669675624e1, -0.753533046979752e-12, 0.481337131452891e1,
1100  -0.223286270422356e22, -0.118746004987383e-4, 0.646412934136496e-2,
1101  -0.410588536330937e-9, 0.422739537057241e20, 0.313698180473812e-12,
1102  0.16439533434504e-23, -0.339823323754373e-5, -0.135268639905021e-1,
1103  -0.723252514211625e-14, 0.184386437538366e-8, -0.463959533752385e-1,
1104  -0.99226310037675e14, 0.688169154439335e-16, -0.222620998452197e-10,
1105  -0.540843018624083e-7, 0.345570606200257e-2, 0.422275800304086e11,
1106  -0.126974478770487e-14, 0.927237985153679e-9, 0.612670812016489e-13,
1107  -0.722693924063497e-11, -0.383669502636822e-3, 0.374684572410204e-3,
1108  -0.931976897511086e5, -0.247690616026922e-1, 0.658110546759474e2},
1109  {-0.111371317395540e-3, 0.100342892423685e1, 0.530615581928979e1,
1110  0.179058760078792e-5, -0.728541958464774e-3, -0.187576133371704e2,
1111  0.199060874071849e-2, 0.243574755377290e2, -0.177040785499444e-3,
1112  -0.25968038522713e-2, -0.198704578406823e3, 0.738627790224287e-4,
1113  -0.236264692844138e-2, -0.161023121314333e1, 0.622322971786473e4,
1114  -0.960754116701669e-8, -0.510572269720488e-10, 0.767373781404211e-2,
1115  0.663855469485254e-14, -0.717590735526745e-9, 0.146564542926508e-4,
1116  0.309029474277013e-11, -0.464216300971708e-15, -0.390499637961161e-13,
1117  -0.236716126781431e-9, 0.454652854268717e-11, -0.422271787482497e-2,
1118  0.283911742354706e-10, 0.270929002720228e1},
1119  {-0.401215699576099e9, 0.484501478318406e11, 0.394721471363678e-14,
1120  0.372629967374147e5, -0.369794374168666e-29, -0.380436407012452e-14,
1121  0.475361629970233e-6, -0.879148916140706e-3, 0.844317863844331,
1122  0.122433162656600e2, -0.104529634830279e3, 0.589702771277429e3,
1123  -0.291026851164444e14, 0.170343072841850e-5, -0.277617606975748e-3,
1124  -0.344709605486686e1, 0.221333862447095e2, -0.194646110037079e3,
1125  0.808354639772825e-15, -0.18084520914547e-10, -0.696664158132412e-5,
1126  -0.181057560300994e-2, 0.255830298579027e1, 0.328913873658481e4,
1127  -0.173270241249904e-18, -0.661876792558034e-6, -0.39568892342125e-2,
1128  0.604203299819132e-17, -0.400879935920517e-13, 0.160751107464958e-8,
1129  0.383719409025556e-4, -0.649565446702457e-14, -0.149095328506e-11,
1130  0.541449377329581e-8},
1131  {0.260702058647537e10, -0.188277213604704e15, 0.554923870289667e19, -0.758966946387758e23,
1132  0.413865186848908e27, -0.81503800073806e12, -0.381458260489955e33, -0.123239564600519e-1,
1133  0.226095631437174e8, -0.49501780950672e12, 0.529482996422863e16, -0.444359478746295e23,
1134  0.521635864527315e35, -0.487095672740742e55, -0.714430209937547e6, 0.127868634615495,
1135  -0.100752127917598e2, 0.777451437960990e7, -0.108105480796471e25, -0.357578581169659e-5,
1136  -0.212857169423484e1, 0.270706111085238e30, -0.695953622348829e33, 0.110609027472280,
1137  0.721559163361354e2, -0.306367307532219e15, 0.265839618885530e-4, 0.253392392889754e-1,
1138  -0.214443041836579e3, 0.937846601489667, 0.223184043101700e1, 0.338401222509191e2,
1139  0.494237237179718e21, -0.198068404154428, -0.141415349881140e31, -0.993862421613651e2,
1140  0.125070534142731e3, -0.996473529004439e3, 0.473137909872765e5, 0.116662121219322e33,
1141  -0.315874976271533e16, -0.445703369196945e33, 0.642794932373694e33},
1142  {0.811384363481847, -0.568199310990094e4, -0.178657198172556e11, 0.795537657613427e32,
1143  -0.814568209346872e5, -0.659774567602874e8, -0.152861148659302e11, -0.560165667510446e12,
1144  0.458384828593949e6, -0.385754000383848e14, 0.453735800004273e8, 0.939454935735563e12,
1145  0.266572856432938e28, -0.547578313899097e10, 0.200725701112386e15, 0.185007245563239e13,
1146  0.185135446828337e9, -0.170451090076385e12, 0.157890366037614e15, -0.202530509748774e16,
1147  0.36819392618357e60, 0.170215539458936e18, 0.639234909918741e42, -0.821698160721956e15,
1148  -0.795260241872306e24, 0.23341586947851e18, -0.600079934586803e23, 0.594584382273384e25,
1149  0.189461279349492e40, -0.810093428842645e46, 0.188813911076809e22, 0.111052244098768e36,
1150  0.291133958602503e46, -0.329421923951460e22, -0.137570282536696e26, 0.181508996303902e28,
1151  -0.346865122768353e30, -0.21196114877426e38, -0.128617899887675e49, 0.479817895699239e65},
1152  {0.280967799943151e-38, 0.614869006573609e-30, 0.582238667048942e-27,
1153  0.390628369238462e-22, 0.821445758255119e-20, 0.402137961842776e-14,
1154  0.651718171878301e-12, -0.211773355803058e-7, 0.264953354380072e-2,
1155  -0.135031446451331e-31, -0.607246643970893e-23, -0.402352115234494e-18,
1156  -0.744938506925544e-16, 0.189917206526237e-12, 0.364975183508473e-5,
1157  0.177274872361946e-25, -0.334952758812999e-18, -0.421537726098389e-8,
1158  -0.391048167929649e-1, 0.541276911564176e-13, 0.705412100773699e-11,
1159  0.258585887897486e-8, -0.493111362030162e-10, -0.158649699894543e-5,
1160  -0.525037427886100, 0.220019901729615e-2, -0.643064132636925e-2,
1161  0.629154149015048e2, 0.135147318617061e3, 0.240560808321713e-6,
1162  -0.890763306701305e-3, -0.440209599407714e4, -0.302807107747776e3,
1163  0.159158748314599e4, 0.232534272709876e6, -0.792681207132600e6,
1164  -0.869871364662769e11, 0.354542769185671e12, 0.400849240129329e15},
1165  {0.128746023979718e-34, -0.735234770382342e-11, 0.28907869214915e-2,
1166  0.244482731907223, 0.141733492030985e-23, -0.354533853059476e-28,
1167  -0.594539202901431e-17, -.585188401782779e-8, .201377325411803e-5,
1168  0.138647388209306e1, -0.173959365084772e-4, 0.137680878349369e-2,
1169  0.814897605805513e-14, 0.425596631351839e-25, -0.387449113787755e-17,
1170  0.13981474793024e-12, -0.171849638951521e-2, 0.641890529513296e-21,
1171  0.118960578072018e-10, -0.155282762571611e-17, 0.233907907347507e-7,
1172  -0.174093247766213e-12, 0.377682649089149e-8, -0.516720236575302e-10},
1173  {-0.982825342010366e-4, 0.105145700850612e1, 0.116033094095084e3, 0.324664750281543e4,
1174  -0.123592348610137e4, -0.561403450013495e-1, 0.856677401640869e-7, 0.236313425393924e3,
1175  0.972503292350109e-2, -0.103001994531927e1, -0.149653706199162e-8, -0.215743778861592e-4,
1176  -0.834452198291445e1, 0.586602660564988, 0.343480022104968e-25, 0.816256095947021e-5,
1177  0.294985697916798e-2, 0.711730466276584e-16, 0.400954763806941e-9, 0.107766027032853e2,
1178  -0.409449599138182e-6, -0.729121307758902e-5, 0.677107970938909e-8, 0.602745973022975e-7,
1179  -0.382323011855257e-10, 0.179946628317437e-2, -0.345042834640005e-3},
1180  {-0.820433843259950e5, 0.473271518461586e11, -0.805950021005413e-1, 0.328600025435980e2,
1181  -0.35661702998249e4, -0.172985781433335e10, 0.351769232729192e8, -0.775489259985144e6,
1182  0.710346691966018e-4, 0.993499883820274e5, -0.642094171904570, -0.612842816820083e4,
1183  0.232808472983776e3, -0.142808220416837e-4, -0.643596060678456e-2, -0.428577227475614e1,
1184  0.225689939161918e4, 0.100355651721510e-2, 0.333491455143516, 0.109697576888873e1,
1185  0.961917379376452, -0.838165632204598e-1, 0.247795908411492e1, -0.319114969006533e4},
1186  {0.144165955660863e-2, -0.701438599628258e13, -0.830946716459219e-16, 0.261975135368109,
1187  0.393097214706245e3, -0.104334030654021e5, 0.490112654154211e9, -0.147104222772069e-3,
1188  0.103602748043408e1, 0.305308890065089e1, -0.399745276971264e7, 0.569233719593750e-11,
1189  -0.464923504407778e-1, -0.535400396512906e-17, 0.399988795693162e-12, -0.536479560201811e-6,
1190  0.159536722411202e-1, 0.270303248860217e-14, 0.244247453858506e-7, -0.983430636716454e-5,
1191  0.663513144224454e-1, -0.993456957845006e1, 0.546491323528491e3, -0.143365406393758e5,
1192  0.150764974125511e6, -0.337209709340105e-9, 0.377501980025469e-8},
1193  {-0.532466612140254e23, 0.100415480000824e32, -0.191540001821367e30, 0.105618377808847e17,
1194  0.202281884477061e59, 0.884585472596134e8, 0.166540181638363e23, -0.313563197669111e6,
1195  -0.185662327545324e54, -0.624942093918942e-1, -0.50416072413259e10, 0.187514491833092e5,
1196  0.121399979993217e-2, 0.188317043049455e1, -0.167073503962060e4, 0.965961650599775,
1197  0.294885696802488e1, -0.653915627346115e5, 0.604012200163444e50, -0.198339358557937,
1198  -0.175984090163501e58, 0.356314881403987e1, -0.575991255144384e3, 0.456213415338071e5,
1199  -0.109174044987829e8, 0.437796099975134e34, -0.616552611135792e46, 0.193568768917797e10,
1200  0.950898170425042e54},
1201  {0.155287249586268e1, 0.664235115009031e1, -0.289366236727210e4, -0.385923202309848e13,
1202  -0.291002915783761e1, -0.829088246858083e12, 0.176814899675218e1, -0.534686695713469e9,
1203  0.160464608687834e18, 0.196435366560186e6, 0.156637427541729e13, -0.178154560260006e1,
1204  -0.229746237623692e16, 0.385659001648006e8, 0.110554446790543e10, -0.677073830687349e14,
1205  -0.327910592086523e31, -0.341552040860644e51, -0.527251339709047e21, 0.245375640937055e24,
1206  -0.168776617209269e27, 0.358958955867578e29, -0.656475280339411e36, 0.355286045512301e39,
1207  0.569021454413270e58, -0.700584546433113e48, -0.705772623326374e65, 0.166861176200148e53,
1208  -0.300475129680486e61, -0.668481295196808e51, 0.428432338620678e69, -0.444227367758304e72,
1209  -0.281396013562745e77},
1210  {0.122088349258355e18, 0.104216468608488e10, -.882666931564652e16, 0.259929510849499e20,
1211  0.222612779142211e15, -0.878473585050085e18, -0.314432577551552e22, -0.216934916996285e13,
1212  0.159079648196849e21, -0.339567617303423e3, 0.884387651337836e13, -0.843405926846418e21,
1213  0.114178193518022e2, -0.122708229235641e-3, -0.106201671767107e3, 0.903443213959313e25,
1214  -0.693996270370852e28, 0.648916718965575e-8, 0.718957567127851e4, 0.105581745346187e-2,
1215  -0.651903203602581e15, -0.160116813274676e25, -0.510254294237837e-8, -0.152355388953402,
1216  0.677143292290144e12, 0.276378438378930e15, 0.116862983141686e-1, -0.301426947980171e14,
1217  0.169719813884840e-7, 0.104674840020929e27, -0.10801690456014e5, -0.990623601934295e-12,
1218  0.536116483602738e7, 0.226145963747881e22, -0.488731565776210e-9, 0.151001548880670e-4,
1219  -0.227700464643920e5, -0.781754507698846e28},
1220  {-0.415652812061591e-54, 0.177441742924043e-60, -0.357078668203377e-54,
1221  0.359252213604114e-25, -0.259123736380269e2, 0.594619766193460e5,
1222  -0.624184007103158e11, 0.313080299915944e17, 0.105006446192036e-8,
1223  -0.192824336984852e-5, 0.654144373749937e6, 0.513117462865044e13,
1224  -0.697595750347391e19, -0.103977184454767e29, 0.119563135540666e-47,
1225  -0.436677034051655e-41, 0.926990036530639e-29, 0.587793105620748e21,
1226  0.280375725094731e-17, -0.192359972440634e23, 0.742705723302738e27,
1227  -0.517429682450605e2, 0.820612048645469e7, -0.188214882341448e-8,
1228  0.184587261114837e-1, -0.135830407782663e-5, -0.723681885626348e17,
1229  -0.223449194054124e27, -0.111526741826431e-34, 0.276032601145151e-28,
1230  0.134856491567853e15, 0.652440293345860e-9, 0.510655119774360e17,
1231  -0.468138358908732e32, -0.760667491183279e16, -0.417247986986821e-18,
1232  0.312545677756104e14, -0.100375333864186e15, 0.247761392329058e27},
1233  {-0.586219133817016e-7, -0.894460355005526e11, 0.531168037519774e-30,
1234  0.109892402329239, -0.575368389425212e-1, 0.228276853990249e5,
1235  -0.158548609655002e19, 0.329865748576503e-27, -0.634987981190669e-24,
1236  0.615762068640611e-8, -0.961109240985747e8, -0.406274286652625e-44,
1237  -0.471103725498077e-12, 0.725937724828145, 0.187768525763682e-38,
1238  -0.103308436323771e4, -0.662552816342168e-1, 0.579514041765710e3,
1239  0.237416732616644e-26, 0.271700235739893e-14, -0.9078862134836e2,
1240  -0.171242509570207e-36, 0.156792067854621e3, 0.923261357901470,
1241  -0.597865988422577e1, 0.321988767636389e7, -0.399441390042203e-29,
1242  0.493429086046981e-7, 0.812036983370565e-19, -0.207610284654137e-11,
1243  -0.340821291419719e-6, 0.542000573372233e-17, -0.856711586510214e-12,
1244  0.266170454405981e-13, 0.858133791857099e-5},
1245  {0.377373741298151e19, -0.507100883722913e13, -0.10336322559886e16, 0.184790814320773e-5,
1246  -0.924729378390945e-3, -0.425999562292738e24, -0.462307771873973e-12, 0.107319065855767e22,
1247  0.648662492280682e11, 0.244200600688281e1, -0.851535733484258e10, 0.169894481433592e22,
1248  0.215780222509020e-26, -0.320850551367334, -0.382642448458610e17, -0.275386077674421e-28,
1249  -0.563199253391666e6, -0.326068646279314e21, 0.397949001553184e14, 0.100824008584757e-6,
1250  0.162234569738433e5, -0.432355225319745e11, -0.59287424559861e12, 0.133061647281106e1,
1251  0.157338197797544e7, 0.258189614270853e14, 0.262413209706358e25, -0.920011937431142e-1,
1252  0.220213765905426e-2, -0.110433759109547e2, 0.847004870612087e7, -0.592910695762536e9,
1253  -0.183027173269660e-4, 0.181339603516302, -0.119228759669889e4, 0.430867658061468e7},
1254  {-0.525597995024633e-9, 0.583441305228407e4, -0.134778968457925e17, 0.118973500934212e26,
1255  -0.159096490904708e27, -0.315839902302021e-6, 0.496212197158239e3, 0.327777227273171e19,
1256  -0.527114657850696e22, 0.210017506281863e-16, 0.705106224399834e21, -0.266713136106469e31,
1257  -0.145370512554562e-7, 0.149333917053130e28, -0.149795620287641e8, -0.3818819062711e16,
1258  0.724660165585797e-4, -0.937808169550193e14, 0.514411468376383e10, -0.828198594040141e5},
1259  {0.24400789229065e-10, -0.463057430331242e7, 0.728803274777712e10, 0.327776302858856e16,
1260  -0.110598170118409e10, -0.323899915729957e13, 0.923814007023245e16, 0.842250080413712e-12,
1261  0.663221436245506e12, -0.167170186672139e15, 0.253749358701391e4, -0.819731559610523e-20,
1262  0.328380587890663e12, -0.625004791171543e8, 0.803197957462023e21, -0.204397011338353e-10,
1263  -0.378391047055938e4, 0.97287654593862e-2, 0.154355721681459e2, -0.373962862928643e4,
1264  -0.682859011374572e11, -0.248488015614543e-3, 0.394536049497068e7}};
1265 
1266  const std::vector<std::vector<int>> _I3s{
1267  {-12, -12, -12, -10, -10, -10, -8, -8, -8, -6, -5, -5, -5, -4, -3,
1268  -3, -3, -3, -2, -2, -2, -1, -1, -1, 0, 0, 1, 1, 2, 2},
1269  {-12, -12, -10, -10, -8, -6, -6, -6, -5, -5, -5, -4, -4, -4, -3, -3,
1270  -3, -3, -3, -2, -2, -2, -1, -1, 0, 0, 1, 1, 2, 3, 4, 4},
1271  {-12, -12, -12, -10, -10, -10, -8, -8, -8, -6, -5, -5, -5, -4, -4, -3, -3, -2,
1272  -2, -2, -1, -1, -1, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 8},
1273  {-12, -12, -12, -12, -12, -12, -10, -10, -10, -10, -10, -10, -10, -8, -8, -8, -8, -6, -6,
1274  -5, -5, -5, -5, -4, -4, -4, -3, -3, -2, -2, -1, -1, -1, 0, 0, 1, 1, 3},
1275  {-12, -12, -10, -10, -10, -10, -10, -8, -8, -8, -6, -5, -4, -4, -3,
1276  -3, -3, -2, -2, -2, -2, -1, 0, 0, 1, 1, 1, 2, 2},
1277  {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 5, 5, 6, 7, 7,
1278  10, 12, 12, 12, 14, 14, 14, 14, 14, 16, 16, 18, 18, 20, 20, 20, 22, 24, 24, 28, 32},
1279  {-12, -12, -12, -12, -12, -12, -10, -10, -10, -8, -8, -8, -8, -6, -6, -5, -5, -4, -3,
1280  -2, -2, -2, -2, -1, -1, -1, 0, 0, 0, 1, 1, 1, 3, 5, 6, 8, 10, 10},
1281  {-12, -12, -10, -10, -10, -10, -10, -10, -8, -8, -8, -8, -8, -6, -6,
1282  -6, -5, -5, -5, -4, -4, -3, -3, -2, -1, -1, 0, 1, 1},
1283  {0, 0, 0, 1, 1, 1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 7, 7, 8, 8, 10,
1284  12, 12, 12, 14, 14, 14, 14, 18, 18, 18, 18, 18, 20, 20, 22, 24, 24, 32, 32, 36, 36},
1285  {0, 0, 0, 1, 1, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6,
1286  10, 12, 12, 14, 14, 14, 16, 18, 20, 20, 24, 24, 28, 28},
1287  {-2, -2, -1, -1, 0, -0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
1288  1, 2, 2, 2, 2, 2, 2, 5, 5, 5, 6, 6, 6, 6, 8, 10, 12},
1289  {-12, -12, -12, -12, -12, -10, -10, -8, -8, -8, -8, -8, -8, -8, -6,
1290  -5, -5, -4, -4, -3, -3, -3, -3, -2, -2, -2, -1, -1, -1, 0,
1291  0, 0, 0, 1, 1, 2, 4, 5, 5, 6, 10, 10, 14},
1292  {0, 3, 8, 20, 1, 3, 4, 5, 1, 6, 2, 4, 14, 2, 5, 3, 0, 1, 1, 1,
1293  28, 2, 16, 0, 5, 0, 3, 4, 12, 16, 1, 8, 14, 0, 2, 3, 4, 8, 14, 24},
1294  {0, 3, 4, 6, 7, 10, 12, 14, 18, 0, 3, 5, 6, 8, 12, 0, 3, 7, 12, 2,
1295  3, 4, 2, 4, 7, 4, 3, 5, 6, 0, 0, 3, 1, 0, 1, 0, 1, 0, 1},
1296  {0, 0, 0, 2, 3, 4, 4, 4, 4, 4, 5, 5, 6, 7, 8, 8, 8, 10, 10, 14, 14, 20, 20, 24},
1297  {0, 0, 0, 0, 1, 2, 3, 3, 4, 6, 7, 7, 8, 10,
1298  12, 12, 12, 14, 14, 14, 16, 18, 20, 22, 24, 24, 36},
1299  {-12, -12, -10, -10, -10, -10, -8, -6, -5, -5, -4, -4,
1300  -3, -2, -2, -2, -2, -1, -1, -1, 0, 1, 1, 1},
1301  {-8, -8, -3, -3, -3, -3, -3, 0, 0, 0, 0, 3, 3, 8,
1302  8, 8, 8, 10, 10, 10, 10, 10, 10, 10, 10, 12, 14},
1303  {-12, -12, -10, -8, -6, -5, -5, -4, -4, -3, -3, -2, -1, -1, -1,
1304  0, 0, 0, 0, 1, 1, 3, 3, 3, 4, 4, 4, 5, 14},
1305  {0, 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 7, 7, 7, 7,
1306  7, 10, 10, 10, 10, 10, 18, 20, 22, 22, 24, 28, 32, 32, 32, 36},
1307  {-12, -10, -10, -10, -8, -8, -8, -6, -6, -5, -5, -5, -3, -1, -1, -1, -1, 0, 0,
1308  1, 2, 2, 3, 5, 5, 5, 6, 6, 8, 8, 10, 12, 12, 12, 14, 14, 14, 14},
1309  {-10, -8, -6, -6, -6, -6, -6, -6, -5, -5, -5, -5, -5, -5, -4, -4, -4, -4, -3, -3,
1310  -3, -2, -2, -1, -1, 0, 0, 0, 1, 1, 3, 4, 4, 4, 5, 8, 10, 12, 14},
1311  {-12, -12, -10, -10, -8, -8, -8, -6, -6, -6, -6, -5, -4, -4, -3, -3, -2, -2,
1312  -1, -1, -1, 0, 0, 1, 2, 2, 3, 3, 5, 5, 5, 8, 8, 10, 10},
1313  {-8, -6, -5, -4, -4, -4, -3, -3, -1, 0, 0, 0, 1, 1, 2, 3, 3, 3,
1314  4, 5, 5, 5, 6, 8, 8, 8, 8, 10, 12, 12, 12, 12, 14, 14, 14, 14},
1315  {0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 8, 8, 10, 12},
1316  {-8, -6, -5, -5, -4, -4, -4, -3, -3, -3, -2, -1, 0, 1, 2, 3, 3, 6, 6, 6, 6, 8, 8}};
1317 
1318  const std::vector<std::vector<int>> _J3s{
1319  {5, 10, 12, 5, 10, 12, 5, 8, 10, 1, 1, 5, 10, 8, 0,
1320  1, 3, 6, 0, 2, 3, 0, 1, 2, 0, 1, 0, 2, 0, 2},
1321  {10, 12, 8, 14, 8, 5, 6, 8, 5, 8, 10, 2, 4, 5, 0, 1,
1322  2, 3, 5, 0, 2, 5, 0, 2, 0, 1, 0, 2, 0, 2, 0, 1},
1323  {6, 8, 10, 6, 8, 10, 5, 6, 7, 8, 1, 4, 7, 2, 8, 0, 3, 0,
1324  4, 5, 0, 1, 2, 0, 1, 2, 0, 2, 0, 1, 3, 7, 0, 7, 1},
1325  {4, 6, 7, 10, 12, 16, 0, 2, 4, 6, 8, 10, 14, 3, 7, 8, 10, 6, 8,
1326  1, 2, 5, 7, 0, 1, 7, 2, 4, 0, 1, 0, 1, 5, 0, 2, 0, 6, 0},
1327  {14, 16, 3, 6, 10, 14, 16, 7, 8, 10, 6, 6, 2, 4, 2, 6, 7, 0, 1, 3, 4, 0, 0, 1, 0, 4, 6, 0, 2},
1328  {-3, -2, -1, 0, 1, 2, -1, 1, 2, 3, 0, 1, -5, -2,
1329  0, -3, -8, 1, -6, -4, 1, -6, -10, -8, -4, -12, -10, -8,
1330  -6, -4, -10, -8, -12, -10, -12, -10, -6, -12, -12, -4, -12, -12},
1331  {7, 12, 14, 18, 22, 24, 14, 20, 24, 7, 8, 10, 12, 8, 22, 7, 20, 22, 7,
1332  3, 5, 14, 24, 2, 8, 18, 0, 1, 2, 0, 1, 3, 24, 22, 12, 3, 0, 6},
1333  {8, 12, 4, 6, 8, 10, 14, 16, 0, 1, 6, 7, 8, 4, 6, 8, 2, 3, 4, 2, 4, 1, 2, 0, 0, 2, 0, 0, 2},
1334  {0, 1, 10, -4, -2, -1, 0, 0, -5, 0, -3, -2, -1, -6, -1, 12, -4, -3, -6, 10, -8,
1335  -12, -6, -4, -10, -8, -4, 5, -12, -10, -8, -6, 2, -12, -10, -12, -12, -8, -10, -5, -10, -8},
1336  {-1, 0, 1, -2, -1, 1, -1, 1, -2, -2, 2, -3, -2, 0, 3,
1337  -6, -8, -3, -10, -8, -5, -10, -12, -12, -10, -12, -6, -12, -5},
1338  {10, 12, -5, 6, -12, -6, -2, -1, 0, 1, 2, 3, 14, -3, -2, 0, 1,
1339  2, -8, -6, -3, -2, 0, 4, -12, -6, -3, -12, -10, -8, -5, -12, -12, -10},
1340  {14, 16, 18, 20, 22, 14, 24, 6, 10, 12, 14, 18, 24, 36, 8, 4, 5, 7, 16, 1, 3, 18,
1341  20, 2, 3, 10, 0, 1, 3, 0, 1, 2, 12, 0, 16, 1, 0, 0, 1, 14, 4, 12, 10},
1342  {0, 0, 0, 2, 5, 5, 5, 5, 6, 6, 7, 8, 8, 10, 10, 12, 14, 14, 18, 20,
1343  20, 22, 22, 24, 24, 28, 28, 28, 28, 28, 32, 32, 32, 36, 36, 36, 36, 36, 36, 36},
1344  {-12, -12, -12, -12, -12, -12, -12, -12, -12, -10, -10, -10, -10,
1345  -10, -10, -8, -8, -8, -8, -6, -6, -6, -5, -5, -5, -4,
1346  -3, -3, -3, -2, -1, -1, 0, 1, 1, 2, 4, 5, 6},
1347  {-12, -4, -1, -1, -10, -12, -8, -5, -4, -1, -4, -3,
1348  -8, -12, -10, -8, -4, -12, -8, -12, -8, -12, -10, -12},
1349  {-1, 0, 1, 2, 1, -1, -3, 0, -2, -2, -5, -4, -2, -3,
1350  -12, -6, -5, -10, -8, -3, -8, -8, -10, -10, -12, -8, -12},
1351  {10, 12, 6, 7, 8, 10, 8, 6, 2, 5, 3, 4, 3, 0, 1, 2, 4, 0, 1, 2, 0, 0, 1, 3},
1352  {6, 14, -3, 3, 4, 5, 8, -1, 0, 1, 5, -6, -2, -12,
1353  -10, -8, -5, -12, -10, -8, -6, -5, -4, -3, -2, -12, -12},
1354  {20, 24, 22, 14, 36, 8, 16, 6, 32, 3, 8, 4, 1, 2, 3,
1355  0, 1, 4, 28, 0, 32, 0, 1, 2, 3, 18, 24, 4, 24},
1356  {0, 1, 4, 12, 0, 10, 0, 6, 14, 3, 8, 0, 10, 3, 4, 7, 20,
1357  36, 10, 12, 14, 16, 22, 18, 32, 22, 36, 24, 28, 22, 32, 36, 36},
1358  {14, 10, 12, 14, 10, 12, 14, 8, 12, 4, 8, 12, 2, -1, 1, 12, 14, -3, 1,
1359  -2, 5, 10, -5, -4, 2, 3, -5, 2, -8, 8, -4, -12, -4, 4, -12, -10, -6, 6},
1360  {-8, -12, -12, -3, 5, 6, 8, 10, 1, 2, 6, 8, 10, 14, -12, -10, -6, 10, -3, 10,
1361  12, 2, 4, -2, 0, -2, 6, 10, -12, -10, 3, -6, 3, 10, 2, -12, -2, -3, 1},
1362  {8, 14, -1, 8, 6, 8, 14, -4, -3, 2, 8, -10, -1, 3, -10, 3, 1, 2,
1363  -8, -4, 1, -12, 1, -1, -1, 2, -12, -5, -10, -8, -6, -12, -10, -12, -8},
1364  {14, 10, 10, 1, 2, 14, -2, 12, 5, 0, 4, 10, -10, -1, 6, -12, 0, 8,
1365  3, -6, -2, 1, 1, -6, -3, 1, 8, -8, -10, -8, -5, -4, -12, -10, -8, -6},
1366  {-3, 1, 5, 8, 8, -4, -1, 4, 5, -8, 4, 8, -6, 6, -2, 1, -8, -2, -5, -8},
1367  {3, 6, 6, 8, 5, 6, 8, -2, 5, 6, 2, -6, 3, 1, 6, -6, -2, -6, -5, -4, -1, -8, -4}};
1368 
1369  const std::array<Real, 31> _nph3a{
1370  {-0.133645667811215e-6, 0.455912656802978e-5, -0.146294640700979e-4, 0.639341312970080e-2,
1371  0.372783927268847e3, -0.718654377460447e4, 0.573494752103400e6, -0.267569329111439e7,
1372  -0.334066283302614e-4, -0.245479214069597e-1, 0.478087847764996e2, 0.764664131818904e-5,
1373  0.128350627676972e-2, 0.171219081377331e-1, -0.851007304583213e1, -0.136513461629781e-1,
1374  -0.384460997596657e-5, 0.337423807911655e-2, -0.551624873066791, 0.729202277107470,
1375  -0.992522757376041e-2, -.119308831407288, .793929190615421, .454270731799386,
1376  .20999859125991, -0.642109823904738e-2, -0.235155868604540e-1, 0.252233108341612e-2,
1377  -0.764885133368119e-2, 0.136176427574291e-1, -0.133027883575669e-1}};
1378 
1379  const std::array<int, 31> _Iph3a{{-12, -12, -12, -12, -12, -12, -12, -12, -10, -10, -10,
1380  -8, -8, -8, -8, -5, -3, -2, -2, -2, -1, -1,
1381  0, 0, 1, 3, 3, 4, 4, 10, 12}};
1382 
1383  const std::array<int, 31> _Jph3a{{0, 1, 2, 6, 14, 16, 20, 22, 1, 5, 12, 0, 2, 4, 10, 2,
1384  0, 1, 3, 4, 0, 2, 0, 1, 1, 0, 1, 0, 3, 4, 5}};
1385 
1386  const std::array<Real, 33> _nph3b{
1387  {0.323254573644920e-4, -0.127575556587181e-3, -0.475851877356068e-3, 0.156183014181602e-2,
1388  0.105724860113781, -0.858514221132534e2, 0.724140095480911e3, 0.296475810273257e-2,
1389  -0.592721983365988e-2, -0.126305422818666e-1, -0.115716196364853, 0.849000969739595e2,
1390  -0.108602260086615e-1, 0.154304475328851e-1, 0.750455441524466e-1, 0.252520973612982e-1,
1391  -0.602507901232996e-1, -0.307622221350501e1, -0.574011959864879e-1, 0.503471360939849e1,
1392  -0.925081888584834, 0.391733882917546e1, -0.773146007130190e2, 0.949308762098587e4,
1393  -0.141043719679409e7, 0.849166230819026e7, 0.861095729446704, 0.323346442811720,
1394  0.873281936020439, -0.436653048526683, 0.286596714529479, -0.131778331276228,
1395  0.676682064330275e-2}};
1396 
1397  const std::array<int, 33> _Iph3b{{-12, -12, -10, -10, -10, -10, -10, -8, -8, -8, -8,
1398  -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1,
1399  -1, -1, -1, -1, 0, 0, 1, 3, 5, 6, 8}};
1400 
1401  const std::array<int, 33> _Jph3b{{0, 1, 0, 1, 5, 10, 12, 0, 1, 2, 4, 10, 0, 1, 2, 0, 1,
1402  5, 0, 4, 2, 4, 6, 10, 14, 16, 0, 2, 1, 1, 1, 1, 1}};
1403 
1405  const std::array<Real, 10> _n4{{0.11670521452767e4,
1406  -0.72421316703206e6,
1407  -0.17073846940092e2,
1408  0.12020824702470e5,
1409  -0.32325550322333e7,
1410  0.14915108613530e2,
1411  -0.48232657361591e4,
1412  0.40511340542057e6,
1413  -0.238555575678490,
1414  0.65017534844798e3}};
1415 
1417  const std::array<int, 6> _J05{{0, 1, -3, -2, -1, 2}};
1418  const std::array<Real, 6> _n05{{-0.13179983674201e2,
1419  0.68540841634434e1,
1420  -0.24805148933466e-1,
1421  0.36901534980333,
1422  -0.31161318213925e1,
1423  -0.32961626538917}};
1424 
1425  const std::array<int, 6> _I5{{1, 1, 1, 2, 2, 3}};
1426  const std::array<int, 6> _J5{{1, 2, 3, 3, 9, 7}};
1427  const std::array<Real, 6> _n5{{0.15736404855259e-2,
1428  0.90153761673944e-3,
1429  -0.50270077677648e-2,
1430  0.22440037409485e-5,
1431  -0.41163275453471e-5,
1432  0.37919454822955e-7}};
1433 
1435  const std::vector<std::vector<int>> _tempXY_I{{0, 1, 2, -1, -2},
1436  {0, 1, 2, 3},
1437  {0, 1, 2, 3, 4},
1438  {0, 1, 2, 3, 4},
1439  {0, 1, 2, 3, 4},
1440  {0, 1, 2, 3},
1441  {0, 1, 2, -1, -2},
1442  {0, 1, 2, 3},
1443  {0, 1, 2, 3},
1444  {0, 1, 2, 3, 4},
1445  {0, 1, 2, -1, -2}};
1446 
1447  const std::vector<std::vector<Real>> _tempXY_n{
1448  {0.154793642129415e4,
1449  -0.187661219490113e3,
1450  0.213144632222113e2,
1451  -0.191887498864292e4,
1452  0.918419702359447e3},
1453  {0.585276966696349e3, 0.278233532206915e1, -0.127283549295878e-1, 0.159090746562729e-3},
1454  {-0.249284240900418e5,
1455  0.428143584791546e4,
1456  -0.269029173140130e3,
1457  0.751608051114157e1,
1458  -0.787105249910383e-1},
1459  {0.584814781649163e3,
1460  -0.616179320924617,
1461  0.260763050899562,
1462  -0.587071076864459e-2,
1463  0.515308185433082e-4},
1464  {0.617229772068439e3,
1465  -0.770600270141675e1,
1466  0.697072596851896,
1467  -0.157391839848015e-1,
1468  0.137897492684194e-3},
1469  {0.535339483742384e3, 0.761978122720128e1, -0.158365725441648, 0.192871054508108e-2},
1470  {0.969461372400213e3,
1471  -0.332500170441278e3,
1472  0.642859598466067e2,
1473  0.773845935768222e3,
1474  -0.152313732937084e4},
1475  {0.565603648239126e3, 0.529062258221222e1, -0.102020639611016, 0.122240301070145e-2},
1476  {0.584561202520006e3, -0.102961025163669e1, 0.243293362700452, -0.294905044740799e-2},
1477  {0.528199646263062e3, 0.890579602135307e1, -0.222814134903755, 0.286791682263697e-2},
1478  {0.728052609145380e1,
1479  0.973505869861952e2,
1480  0.147370491183191e2,
1481  0.329196213998375e3,
1482  0.873371668682417e3}};
1483 
1486  const std::array<Real, 4> _mu_H0{{1.67752, 2.20462, 0.6366564, -0.241605}};
1487  const std::array<std::array<Real, 7>, 6> _mu_Hij{
1488  {{{5.20094e-1, 2.22531e-1, -2.81378e-1, 1.61913e-1, -3.25372e-2, 0.0, 0.0}},
1489  {{8.50895e-2, 9.99115e-1, -9.06851e-1, 2.57399e-1, 0.0, 0.0, 0.0}},
1490  {{-1.08374, 1.88797, -7.72479e-1, 0.0, 0.0, 0.0, 0.0}},
1491  {{-2.89555e-1, 1.26613, -4.89837e-1, 0.0, 6.98452e-2, 0.0, -4.35673e-3}},
1492  {{0.0, 0.0, -2.57040e-1, 0.0, 0.0, 8.72102e-3, 0.0}},
1493  {{0.0, 1.20573e-1, 0.0, 0.0, 0.0, 0.0, -5.93264e-4}}}};
1494 
1496  std::array<Real, 4> _k_a{{0.0102811, 0.0299621, 0.0156146, -0.00422464}};
1497 
1499  const std::array<Real, 5> _T_star{{1386.0, 540.0, _T_critical, 1.0, 1000.0}};
1501  const std::array<Real, 5> _p_star{{16.53e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6}};
1502 
1503 private:
1508  template <typename T>
1509  std::pair<T, T> T_drhodT_from_p_rho(const T & p, const T & rho) const;
1510 };
1511 
1512 #pragma GCC diagnostic pop
1513 
1514 template <typename T>
1515 T
1516 Water97FluidProperties::dgamma1_dpi(const T & pi, const T & tau) const
1517 {
1518  T sum = 0.0;
1519  for (std::size_t i = 0; i < _n1.size(); ++i)
1520  sum += -_n1[i] * _I1[i] * MathUtils::pow(7.1 - pi, _I1[i] - 1) *
1521  MathUtils::pow(tau - 1.222, _J1[i]);
1522 
1523  return sum;
1524 }
1525 
1526 template <typename T>
1527 T
1528 Water97FluidProperties::dgamma2_dpi(const T & pi, const T & tau) const
1529 {
1530  // Ideal gas part of the Gibbs free energy
1531  T dg0 = 1.0 / pi;
1532 
1533  // Residual part of the Gibbs free energy
1534  T dgr = 0.0;
1535  for (std::size_t i = 0; i < _n2.size(); ++i)
1536  dgr += _n2[i] * _I2[i] * MathUtils::pow(pi, _I2[i] - 1) * MathUtils::pow(tau - 0.5, _J2[i]);
1537 
1538  return dg0 + dgr;
1539 }
1540 
1541 template <typename T>
1542 T
1544 {
1545  using std::exp;
1546 
1547  // Region 3 is subdivided into 26 subregions, each with a given backwards equation
1548  // to directly calculate density from pressure and temperature without the need for
1549  // expensive iterations. Find the subregion id that the point is in:
1550  unsigned int sid =
1552 
1553  T vstar, pi, theta;
1554  Real a, b, c, d, e;
1555  unsigned int N;
1556 
1557  vstar = _par3[sid][0];
1558  pi = pressure / _par3[sid][1] / 1.0e6;
1559  theta = temperature / _par3[sid][2];
1560  a = _par3[sid][3];
1561  b = _par3[sid][4];
1562  c = _par3[sid][5];
1563  d = _par3[sid][6];
1564  e = _par3[sid][7];
1565  N = _par3N[sid];
1566 
1567  T sum = 0.0;
1568  T volume = 0.0;
1569 
1570  // Note that subregion 13 is the only different formulation
1571  if (sid == 13)
1572  {
1573  for (std::size_t i = 0; i < N; ++i)
1574  sum += _n3s[sid][i] * MathUtils::pow(pi - a, _I3s[sid][i]) *
1575  MathUtils::pow(theta - b, _J3s[sid][i]);
1576 
1577  volume = vstar * exp(sum);
1578  }
1579  else
1580  volume = vstar * subregionVolume(pi, theta, a, b, c, d, e, sid);
1581 
1582  // Density is the inverse of volume
1583  return 1.0 / volume;
1584 }
1585 
1586 template <typename T>
1587 T
1588 Water97FluidProperties::gamma1(const T & pi, const T & tau) const
1589 {
1590  T sum = 0.0;
1591  for (std::size_t i = 0; i < _n1.size(); ++i)
1592  sum += _n1[i] * MathUtils::pow(7.1 - pi, _I1[i]) * MathUtils::pow(tau - 1.222, _J1[i]);
1593 
1594  return sum;
1595 }
1596 
1597 template <typename T>
1598 T
1599 Water97FluidProperties::d2gamma1_dpi2(const T & pi, const T & tau) const
1600 {
1601  T sum = 0.0;
1602  for (std::size_t i = 0; i < _n1.size(); ++i)
1603  sum += _n1[i] * _I1[i] * (_I1[i] - 1) * MathUtils::pow(7.1 - pi, _I1[i] - 2) *
1604  MathUtils::pow(tau - 1.222, _J1[i]);
1605 
1606  return sum;
1607 }
1608 
1609 template <typename T>
1610 T
1611 Water97FluidProperties::dgamma1_dtau(const T & pi, const T & tau) const
1612 {
1613  T g = 0.0;
1614  for (std::size_t i = 0; i < _n1.size(); ++i)
1615  g += _n1[i] * _J1[i] * MathUtils::pow(7.1 - pi, _I1[i]) *
1616  MathUtils::pow(tau - 1.222, _J1[i] - 1);
1617 
1618  return g;
1619 }
1620 
1621 template <typename T>
1622 T
1623 Water97FluidProperties::d2gamma1_dtau2(const T & pi, const T & tau) const
1624 {
1625  T dg = 0.0;
1626  for (std::size_t i = 0; i < _n1.size(); ++i)
1627  dg += _n1[i] * _J1[i] * (_J1[i] - 1) * MathUtils::pow(7.1 - pi, _I1[i]) *
1628  MathUtils::pow(tau - 1.222, _J1[i] - 2);
1629 
1630  return dg;
1631 }
1632 
1633 template <typename T>
1634 T
1635 Water97FluidProperties::d2gamma1_dpitau(const T & pi, const T & tau) const
1636 {
1637  T dg = 0.0;
1638  for (std::size_t i = 0; i < _n1.size(); ++i)
1639  dg += -_n1[i] * _I1[i] * _J1[i] * MathUtils::pow(7.1 - pi, _I1[i] - 1) *
1640  MathUtils::pow(tau - 1.222, _J1[i] - 1);
1641 
1642  return dg;
1643 }
1644 
1645 template <typename T>
1646 T
1647 Water97FluidProperties::gamma2(const T & pi, const T & tau) const
1648 {
1649  using std::log;
1650 
1651  // Ideal gas part of the Gibbs free energy
1652  T sum0 = 0.0;
1653  for (std::size_t i = 0; i < _n02.size(); ++i)
1654  sum0 += _n02[i] * MathUtils::pow(tau, _J02[i]);
1655 
1656  T g0 = log(pi) + sum0;
1657 
1658  // Residual part of the Gibbs free energy
1659  T gr = 0.0;
1660  for (std::size_t i = 0; i < _n2.size(); ++i)
1661  gr += _n2[i] * MathUtils::pow(pi, _I2[i]) * MathUtils::pow(tau - 0.5, _J2[i]);
1662 
1663  return g0 + gr;
1664 }
1665 
1666 template <typename T>
1667 T
1668 Water97FluidProperties::d2gamma2_dpi2(const T & pi, const T & tau) const
1669 {
1670  // Ideal gas part of the Gibbs free energy
1671  T dg0 = -1.0 / pi / pi;
1672 
1673  // Residual part of the Gibbs free energy
1674  T dgr = 0.0;
1675  for (std::size_t i = 0; i < _n2.size(); ++i)
1676  dgr += _n2[i] * _I2[i] * (_I2[i] - 1) * MathUtils::pow(pi, _I2[i] - 2) *
1677  MathUtils::pow(tau - 0.5, _J2[i]);
1678 
1679  return dg0 + dgr;
1680 }
1681 
1682 template <typename T>
1683 T
1684 Water97FluidProperties::dgamma2_dtau(const T & pi, const T & tau) const
1685 {
1686  // Ideal gas part of the Gibbs free energy
1687  T dg0 = 0.0;
1688  for (std::size_t i = 0; i < _n02.size(); ++i)
1689  dg0 += _n02[i] * _J02[i] * MathUtils::pow(tau, _J02[i] - 1);
1690 
1691  // Residual part of the Gibbs free energy
1692  T dgr = 0.0;
1693  for (std::size_t i = 0; i < _n2.size(); ++i)
1694  dgr += _n2[i] * _J2[i] * MathUtils::pow(pi, _I2[i]) * MathUtils::pow(tau - 0.5, _J2[i] - 1);
1695 
1696  return dg0 + dgr;
1697 }
1698 
1699 template <typename T>
1700 T
1701 Water97FluidProperties::d2gamma2_dtau2(const T & pi, const T & tau) const
1702 {
1703  // Ideal gas part of the Gibbs free energy
1704  T dg0 = 0.0;
1705  for (std::size_t i = 0; i < _n02.size(); ++i)
1706  dg0 += _n02[i] * _J02[i] * (_J02[i] - 1) * MathUtils::pow(tau, _J02[i] - 2);
1707 
1708  // Residual part of the Gibbs free energy
1709  T dgr = 0.0;
1710  for (std::size_t i = 0; i < _n2.size(); ++i)
1711  dgr += _n2[i] * _J2[i] * (_J2[i] - 1) * MathUtils::pow(pi, _I2[i]) *
1712  MathUtils::pow(tau - 0.5, _J2[i] - 2);
1713 
1714  return dg0 + dgr;
1715 }
1716 
1717 template <typename T>
1718 T
1719 Water97FluidProperties::d2gamma2_dpitau(const T & pi, const T & tau) const
1720 {
1721  // Ideal gas part of the Gibbs free energy
1722  T dg0 = 0.0;
1723 
1724  // Residual part of the Gibbs free energy
1725  T dgr = 0.0;
1726  for (std::size_t i = 0; i < _n2.size(); ++i)
1727  dgr += _n2[i] * _I2[i] * _J2[i] * MathUtils::pow(pi, _I2[i] - 1) *
1728  MathUtils::pow(tau - 0.5, _J2[i] - 1);
1729 
1730  return dg0 + dgr;
1731 }
1732 
1733 template <typename T>
1734 T
1735 Water97FluidProperties::phi3(const T & delta, const T & tau) const
1736 {
1737  using std::log;
1738 
1739  T sum = 0.0;
1740  for (std::size_t i = 1; i < _n3.size(); ++i)
1741  sum += _n3[i] * MathUtils::pow(delta, _I3[i]) * MathUtils::pow(tau, _J3[i]);
1742 
1743  return _n3[0] * log(delta) + sum;
1744 }
1745 
1746 template <typename T>
1747 T
1748 Water97FluidProperties::dphi3_ddelta(const T & delta, const T & tau) const
1749 {
1750  T sum = 0.0;
1751  for (std::size_t i = 1; i < _n3.size(); ++i)
1752  sum += _n3[i] * _I3[i] * MathUtils::pow(delta, _I3[i] - 1) * MathUtils::pow(tau, _J3[i]);
1753 
1754  return _n3[0] / delta + sum;
1755 }
1756 
1757 template <typename T>
1758 T
1759 Water97FluidProperties::d2phi3_ddelta2(const T & delta, const T & tau) const
1760 {
1761  T sum = 0.0;
1762  for (std::size_t i = 1; i < _n3.size(); ++i)
1763  sum += _n3[i] * _I3[i] * (_I3[i] - 1) * MathUtils::pow(delta, _I3[i] - 2) *
1764  MathUtils::pow(tau, _J3[i]);
1765 
1766  return -_n3[0] / delta / delta + sum;
1767 }
1768 
1769 template <typename T>
1770 T
1771 Water97FluidProperties::dphi3_dtau(const T & delta, const T & tau) const
1772 {
1773  T sum = 0.0;
1774  for (std::size_t i = 1; i < _n3.size(); ++i)
1775  sum += _n3[i] * _J3[i] * MathUtils::pow(delta, _I3[i]) * MathUtils::pow(tau, _J3[i] - 1);
1776 
1777  return sum;
1778 }
1779 
1780 template <typename T>
1781 T
1782 Water97FluidProperties::d2phi3_dtau2(const T & delta, const T & tau) const
1783 {
1784  T sum = 0.0;
1785  for (std::size_t i = 1; i < _n3.size(); ++i)
1786  sum += _n3[i] * _J3[i] * (_J3[i] - 1) * MathUtils::pow(delta, _I3[i]) *
1787  MathUtils::pow(tau, _J3[i] - 2);
1788 
1789  return sum;
1790 }
1791 
1792 template <typename T>
1793 T
1794 Water97FluidProperties::d2phi3_ddeltatau(const T & delta, const T & tau) const
1795 {
1796  T sum = 0.0;
1797  for (std::size_t i = 1; i < _n3.size(); ++i)
1798  sum += _n3[i] * _I3[i] * _J3[i] * MathUtils::pow(delta, _I3[i] - 1) *
1799  MathUtils::pow(tau, _J3[i] - 1);
1800 
1801  return sum;
1802 }
1803 
1804 template <typename T>
1805 T
1806 Water97FluidProperties::gamma5(const T & pi, const T & tau) const
1807 {
1808  using std::log;
1809 
1810  // Ideal gas part of the Gibbs free energy
1811  T sum0 = 0.0;
1812  for (std::size_t i = 0; i < _n05.size(); ++i)
1813  sum0 += _n05[i] * MathUtils::pow(tau, _J05[i]);
1814 
1815  T g0 = log(pi) + sum0;
1816 
1817  // Residual part of the Gibbs free energy
1818  T gr = 0.0;
1819  for (std::size_t i = 0; i < _n5.size(); ++i)
1820  gr += _n5[i] * MathUtils::pow(pi, _I5[i]) * MathUtils::pow(tau, _J5[i]);
1821 
1822  return g0 + gr;
1823 }
1824 
1825 template <typename T>
1826 T
1827 Water97FluidProperties::dgamma5_dpi(const T & pi, const T & tau) const
1828 {
1829  // Ideal gas part of the Gibbs free energy
1830  T dg0 = 1.0 / pi;
1831 
1832  // Residual part of the Gibbs free energy
1833  T dgr = 0.0;
1834  for (std::size_t i = 0; i < _n5.size(); ++i)
1835  dgr += _n5[i] * _I5[i] * MathUtils::pow(pi, _I5[i] - 1) * MathUtils::pow(tau, _J5[i]);
1836 
1837  return dg0 + dgr;
1838 }
1839 
1840 template <typename T>
1841 T
1842 Water97FluidProperties::d2gamma5_dpi2(const T & pi, const T & tau) const
1843 {
1844  // Ideal gas part of the Gibbs free energy
1845  T dg0 = -1.0 / pi / pi;
1846 
1847  // Residual part of the Gibbs free energy
1848  T dgr = 0.0;
1849  for (std::size_t i = 0; i < _n5.size(); ++i)
1850  dgr += _n5[i] * _I5[i] * (_I5[i] - 1) * MathUtils::pow(pi, _I5[i] - 2) *
1851  MathUtils::pow(tau, _J5[i]);
1852 
1853  return dg0 + dgr;
1854 }
1855 
1856 template <typename T>
1857 T
1858 Water97FluidProperties::dgamma5_dtau(const T & pi, const T & tau) const
1859 {
1860  // Ideal gas part of the Gibbs free energy
1861  T dg0 = 0.0;
1862  for (std::size_t i = 0; i < _n05.size(); ++i)
1863  dg0 += _n05[i] * _J05[i] * MathUtils::pow(tau, _J05[i] - 1);
1864 
1865  // Residual part of the Gibbs free energy
1866  T dgr = 0.0;
1867  for (std::size_t i = 0; i < _n5.size(); ++i)
1868  dgr += _n5[i] * _J5[i] * MathUtils::pow(pi, _I5[i]) * MathUtils::pow(tau, _J5[i] - 1);
1869 
1870  return dg0 + dgr;
1871 }
1872 
1873 template <typename T>
1874 T
1875 Water97FluidProperties::d2gamma5_dtau2(const T & pi, const T & tau) const
1876 {
1877  // Ideal gas part of the Gibbs free energy
1878  T dg0 = 0.0;
1879  for (std::size_t i = 0; i < _n05.size(); ++i)
1880  dg0 += _n05[i] * _J05[i] * (_J05[i] - 1) * MathUtils::pow(tau, _J05[i] - 2);
1881 
1882  // Residual part of the Gibbs free energy
1883  T dgr = 0.0;
1884  for (std::size_t i = 0; i < _n5.size(); ++i)
1885  dgr += _n5[i] * _J5[i] * (_J5[i] - 1) * MathUtils::pow(pi, _I5[i]) *
1886  MathUtils::pow(tau, _J5[i] - 2);
1887 
1888  return dg0 + dgr;
1889 }
1890 
1891 template <typename T>
1892 T
1893 Water97FluidProperties::d2gamma5_dpitau(const T & pi, const T & tau) const
1894 {
1895  // Ideal gas part of the Gibbs free energy
1896  T dg0 = 0.0;
1897 
1898  // Residual part of the Gibbs free energy
1899  T dgr = 0.0;
1900  for (std::size_t i = 0; i < _n5.size(); ++i)
1901  dgr +=
1902  _n5[i] * _I5[i] * _J5[i] * MathUtils::pow(pi, _I5[i] - 1) * MathUtils::pow(tau, _J5[i] - 1);
1903 
1904  return dg0 + dgr;
1905 }
1906 
1907 template <typename T>
1908 T
1910 {
1911  T pi = pressure / 1.0e6;
1912 
1913  // Choose the constants based on the string xy
1914  unsigned int row;
1915 
1916  switch (xy)
1917  {
1918  case AB:
1919  row = 0;
1920  break;
1921  case CD:
1922  row = 1;
1923  break;
1924  case GH:
1925  row = 2;
1926  break;
1927  case IJ:
1928  row = 3;
1929  break;
1930  case JK:
1931  row = 4;
1932  break;
1933  case MN:
1934  row = 5;
1935  break;
1936  case OP:
1937  row = 6;
1938  break;
1939  case QU:
1940  row = 7;
1941  break;
1942  case RX:
1943  row = 8;
1944  break;
1945  case UV:
1946  row = 9;
1947  break;
1948  case WX:
1949  row = 10;
1950  break;
1951  default:
1952  row = 0;
1953  }
1954 
1955  T sum = 0.0;
1956 
1957  using std::log;
1958 
1959  if (xy == AB || xy == OP || xy == WX)
1960  for (std::size_t i = 0; i < _tempXY_n[row].size(); ++i)
1961  sum += _tempXY_n[row][i] * MathUtils::pow(log(pi), _tempXY_I[row][i]);
1962  else if (xy == EF)
1963  sum += 3.727888004 * (pi - _p_critical / 1.0e6) + _T_critical;
1964  else
1965  for (std::size_t i = 0; i < _tempXY_n[row].size(); ++i)
1966  sum += _tempXY_n[row][i] * MathUtils::pow(pi, _tempXY_I[row][i]);
1967 
1968  return sum;
1969 }
1970 
1971 template <typename T>
1972 void
1973 Water97FluidProperties::vaporPressureTemplate(const T & temperature, T & psat, T & dpsat_dT) const
1974 {
1975  using std::sqrt, std::pow;
1976 
1977  // Check whether the input temperature is within the region of validity of this equation.
1978  // Valid for 273.15 K <= t <= 647.096 K
1979  if (temperature < 273.15 || temperature > _T_critical)
1980  mooseException(name(),
1981  ": vaporPressure(): Temperature is outside range 273.15 K <= T <= 647.096 K");
1982 
1983  T theta, dtheta_dT, theta2, a, b, c, da_dtheta, db_dtheta, dc_dtheta;
1984  theta = temperature + _n4[8] / (temperature - _n4[9]);
1985  dtheta_dT = 1.0 - _n4[8] / (temperature - _n4[9]) / (temperature - _n4[9]);
1986  theta2 = theta * theta;
1987 
1988  a = theta2 + _n4[0] * theta + _n4[1];
1989  b = _n4[2] * theta2 + _n4[3] * theta + _n4[4];
1990  c = _n4[5] * theta2 + _n4[6] * theta + _n4[7];
1991 
1992  da_dtheta = 2.0 * theta + _n4[0];
1993  db_dtheta = 2.0 * _n4[2] * theta + _n4[3];
1994  dc_dtheta = 2.0 * _n4[5] * theta + _n4[6];
1995 
1996  T denominator = -b + sqrt(b * b - 4.0 * a * c);
1997 
1998  psat = Utility::pow<4>(2.0 * c / denominator) * 1.0e6;
1999 
2000  // The derivative wrt temperature is given by the chain rule
2001  T dpsat = 4.0 * Utility::pow<3>(2.0 * c / denominator);
2002  dpsat *= (2.0 * dc_dtheta / denominator -
2003  2.0 * c / denominator / denominator *
2004  (-db_dtheta + pow(b * b - 4.0 * a * c, -0.5) *
2005  (b * db_dtheta - 2.0 * da_dtheta * c - 2.0 * a * dc_dtheta)));
2006  dpsat_dT = dpsat * dtheta_dT * 1.0e6;
2007 }
2008 
2009 inline void
2010 Water97FluidProperties::vaporPressure(const Real temperature, Real & psat, Real & dpsat_dT) const
2011 {
2012  vaporPressureTemplate(temperature, psat, dpsat_dT);
2013 }
2014 
2015 template <typename T>
2016 void
2018  const T & pressure, const T & temperature, T & rho, T & drho_dp, T & drho_dT) const
2019 {
2020  auto functor = [this](const auto & pressure, const auto & temperature)
2021  { return this->rho_from_p_T_template(pressure, temperature); };
2022 
2023  xyDerivatives(pressure, temperature, rho, drho_dp, drho_dT, functor);
2024 }
2025 
2026 template <typename T>
2027 T
2029 {
2031 }
2032 
2033 template <typename T>
2034 void
2036  const T & pressure, const T & temperature, T & v, T & dv_dp, T & dv_dT) const
2037 {
2038  auto functor = [this](const auto & pressure, const auto & temperature)
2039  { return this->v_from_p_T_template(pressure, temperature); };
2040 
2041  xyDerivatives(pressure, temperature, v, dv_dp, dv_dT, functor);
2042 }
2043 
2044 template <typename T>
2045 T
2047 {
2048  T internal_energy, pi, tau;
2049 
2050  // Determine which region the point is in
2051  unsigned int region =
2053  switch (region)
2054  {
2055  case 1:
2056  pi = pressure / _p_star[0];
2057  tau = _T_star[0] / temperature;
2058  internal_energy =
2059  _Rw * temperature * (tau * dgamma1_dtau(pi, tau) - pi * dgamma1_dpi(pi, tau));
2060  break;
2061 
2062  case 2:
2063  pi = pressure / _p_star[1];
2064  tau = _T_star[1] / temperature;
2065  internal_energy =
2066  _Rw * temperature * (tau * dgamma2_dtau(pi, tau) - pi * dgamma2_dpi(pi, tau));
2067  break;
2068 
2069  case 3:
2070  {
2071  // Calculate density first, then use that in Helmholtz free energy
2072  T density3 = densityRegion3(pressure, temperature);
2073  T delta = density3 / _rho_critical;
2074  tau = _T_star[2] / temperature;
2075  internal_energy = _Rw * temperature * tau * dphi3_dtau(delta, tau);
2076  break;
2077  }
2078 
2079  case 5:
2080  pi = pressure / _p_star[4];
2081  tau = _T_star[4] / temperature;
2082  internal_energy =
2083  _Rw * temperature * (tau * dgamma5_dtau(pi, tau) - pi * dgamma5_dpi(pi, tau));
2084  break;
2085 
2086  default:
2087  mooseError("inRegion() has given an incorrect region");
2088  }
2089  // Output in J/kg
2090  return internal_energy;
2091 }
2092 
2093 template <typename T>
2094 void
2096  const T & pressure, const T & temperature, T & e, T & de_dp, T & de_dT) const
2097 {
2098  T pi, tau, dinternal_energy_dp, dinternal_energy_dT;
2099 
2100  // Determine which region the point is in
2101  unsigned int region =
2103  switch (region)
2104  {
2105  case 1:
2106  {
2107  pi = pressure / _p_star[0];
2108  tau = _T_star[0] / temperature;
2109  T dgdp = dgamma1_dpi(pi, tau);
2110  T d2gdpt = d2gamma1_dpitau(pi, tau);
2111  dinternal_energy_dp =
2112  _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma1_dpi2(pi, tau)) / _p_star[0];
2113  dinternal_energy_dT =
2114  _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma1_dtau2(pi, tau) - pi * dgdp);
2115  break;
2116  }
2117 
2118  case 2:
2119  {
2120  pi = pressure / _p_star[1];
2121  tau = _T_star[1] / temperature;
2122  T dgdp = dgamma2_dpi(pi, tau);
2123  T d2gdpt = d2gamma2_dpitau(pi, tau);
2124  dinternal_energy_dp =
2125  _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma2_dpi2(pi, tau)) / _p_star[1];
2126  dinternal_energy_dT =
2127  _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma2_dtau2(pi, tau) - pi * dgdp);
2128  break;
2129  }
2130 
2131  case 3:
2132  {
2133  // Calculate density first, then use that in Helmholtz free energy
2134  T density3 = densityRegion3(pressure, temperature);
2135  T delta = density3 / _rho_critical;
2136  tau = _T_star[2] / temperature;
2137  T dpdd = dphi3_ddelta(delta, tau);
2138  T d2pddt = d2phi3_ddeltatau(delta, tau);
2139  T d2pdd2 = d2phi3_ddelta2(delta, tau);
2140  dinternal_energy_dp =
2141  _T_star[2] * d2pddt / _rho_critical /
2142  (2.0 * temperature * delta * dpdd + temperature * delta * delta * d2pdd2);
2143  dinternal_energy_dT =
2144  -_Rw * (delta * tau * d2pddt * (dpdd - tau * d2pddt) / (2.0 * dpdd + delta * d2pdd2) +
2145  tau * tau * d2phi3_dtau2(delta, tau));
2146  break;
2147  }
2148 
2149  case 5:
2150  {
2151  pi = pressure / _p_star[4];
2152  tau = _T_star[4] / temperature;
2153  T dgdp = dgamma5_dpi(pi, tau);
2154  T d2gdpt = d2gamma5_dpitau(pi, tau);
2155  dinternal_energy_dp =
2156  _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma5_dpi2(pi, tau)) / _p_star[4];
2157  dinternal_energy_dT =
2158  _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma5_dtau2(pi, tau) - pi * dgdp);
2159  break;
2160  }
2161 
2162  default:
2163  mooseError("inRegion has given an incorrect region");
2164  }
2165 
2166  e = this->e_from_p_T(pressure, temperature);
2167  de_dp = dinternal_energy_dp;
2168  de_dT = dinternal_energy_dT;
2169 }
2170 
2171 template <typename T>
2172 T
2174  const T & pi, const T & theta, Real a, Real b, Real c, Real d, Real e, unsigned int sid) const
2175 {
2176  using std::pow;
2177 
2178  T sum = 0.0;
2179 
2180  for (std::size_t i = 0; i < _n3s[sid].size(); ++i)
2181  sum += _n3s[sid][i] * MathUtils::pow(pow(pi - a, c), _I3s[sid][i]) *
2182  MathUtils::pow(pow(theta - b, d), _J3s[sid][i]);
2183 
2184  return pow(sum, e);
2185 }
2186 
2187 template <typename T>
2188 std::pair<T, T>
2189 Water97FluidProperties::T_drhodT_from_p_rho(const T & p, const T & rho) const
2190 {
2191  // For NewtonSolve
2192  // y = rho
2193  // x = p
2194  // z = T
2195  auto lambda =
2196  [this](const T & p, const T & temperature, T & rho, T & drho_dp, T & drho_dtemperature)
2197  { rho_from_p_T_template(p, temperature, rho, drho_dp, drho_dtemperature); };
2199  p, rho, _T_initial_guess, _tolerance, lambda, name() + "::T_drhodT_from_p_rho");
2200 }
2201 
2202 template <typename T>
2203 std::pair<T, T>
2204 Water97FluidProperties::p_T_from_v_e(const T & v, const T & e) const
2205 {
2206  const auto rho = 1 / v;
2207  const auto p = p_from_v_e(v, e);
2208  return {p, T_drhodT_from_p_rho(p, rho).first};
2209 }
2210 
2211 template <typename T>
2212 std::pair<T, T>
2213 Water97FluidProperties::rho_T_from_v_e(const T & v, const T & e) const
2214 {
2215  const auto rho = 1 / v;
2216  const auto p = p_from_v_e(v, e);
2217  return {rho, T_drhodT_from_p_rho(p, rho).first};
2218 }
2219 
2220 template <typename T>
2221 T
2223 {
2224  T speed2, pi, tau, delta;
2225 
2226  // Determine which region the point is in
2227  unsigned int region =
2229  switch (region)
2230  {
2231  case 1:
2232  pi = pressure / _p_star[0];
2233  tau = _T_star[0] / temperature;
2234  speed2 = _Rw * temperature * Utility::pow<2>(dgamma1_dpi(pi, tau)) /
2235  (Utility::pow<2>(dgamma1_dpi(pi, tau) - tau * d2gamma1_dpitau(pi, tau)) /
2236  (tau * tau * d2gamma1_dtau2(pi, tau)) -
2237  d2gamma1_dpi2(pi, tau));
2238  break;
2239 
2240  case 2:
2241  pi = pressure / _p_star[1];
2242  tau = _T_star[1] / temperature;
2243  speed2 = _Rw * temperature * Utility::pow<2>(pi * dgamma2_dpi(pi, tau)) /
2244  ((-pi * pi * d2gamma2_dpi2(pi, tau)) +
2245  Utility::pow<2>(pi * dgamma2_dpi(pi, tau) - tau * pi * d2gamma2_dpitau(pi, tau)) /
2246  (tau * tau * d2gamma2_dtau2(pi, tau)));
2247  break;
2248 
2249  case 3:
2250  {
2251  // Calculate density first, then use that in Helmholtz free energy
2252  T density3 = densityRegion3(pressure, temperature);
2253  delta = density3 / _rho_critical;
2254  tau = _T_star[2] / temperature;
2255  speed2 =
2256  _Rw * temperature *
2257  (2.0 * delta * dphi3_ddelta(delta, tau) + delta * delta * d2phi3_ddelta2(delta, tau) -
2258  Utility::pow<2>(delta * dphi3_ddelta(delta, tau) -
2259  delta * tau * d2phi3_ddeltatau(delta, tau)) /
2260  (tau * tau * d2phi3_dtau2(delta, tau)));
2261  break;
2262  }
2263 
2264  case 5:
2265  pi = pressure / _p_star[4];
2266  tau = _T_star[4] / temperature;
2267  speed2 = _Rw * temperature * Utility::pow<2>(pi * dgamma5_dpi(pi, tau)) /
2268  ((-pi * pi * d2gamma5_dpi2(pi, tau)) +
2269  Utility::pow<2>(pi * dgamma5_dpi(pi, tau) - tau * pi * d2gamma5_dpitau(pi, tau)) /
2270  (tau * tau * d2gamma5_dtau2(pi, tau)));
2271  break;
2272 
2273  default:
2274  mooseError("inRegion() has given an incorrect region");
2275  }
2276 
2277  using std::sqrt;
2278  return sqrt(speed2);
2279 }
2280 
2281 template <typename T>
2282 T
2284 {
2285  T specific_heat, pi, tau, delta;
2286 
2287  // Determine which region the point is in
2288  unsigned int region =
2290  switch (region)
2291  {
2292  case 1:
2293  pi = pressure / _p_star[0];
2294  tau = _T_star[0] / temperature;
2295  specific_heat = -_Rw * tau * tau * d2gamma1_dtau2(pi, tau);
2296  break;
2297 
2298  case 2:
2299  pi = pressure / _p_star[1];
2300  tau = _T_star[1] / temperature;
2301  specific_heat = -_Rw * tau * tau * d2gamma2_dtau2(pi, tau);
2302  break;
2303 
2304  case 3:
2305  {
2306  // Calculate density first, then use that in Helmholtz free energy
2307  T density3 = densityRegion3(pressure, temperature);
2308  delta = density3 / _rho_critical;
2309  tau = _T_star[2] / temperature;
2310  specific_heat =
2311  _Rw *
2312  (-tau * tau * d2phi3_dtau2(delta, tau) +
2313  (delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau)) *
2314  (delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau)) /
2315  (2.0 * delta * dphi3_ddelta(delta, tau) +
2316  delta * delta * d2phi3_ddelta2(delta, tau)));
2317  break;
2318  }
2319 
2320  case 5:
2321  pi = pressure / _p_star[4];
2322  tau = _T_star[4] / temperature;
2323  specific_heat = -_Rw * tau * tau * d2gamma5_dtau2(pi, tau);
2324  break;
2325 
2326  default:
2327  mooseError("inRegion() has given an incorrect region");
2328  }
2329  return specific_heat;
2330 }
2331 
2332 template <typename T>
2333 T
2335 {
2336  T specific_heat, pi, tau, delta;
2337 
2338  // Determine which region the point is in
2339  unsigned int region =
2341  switch (region)
2342  {
2343  case 1:
2344  pi = pressure / _p_star[0];
2345  tau = _T_star[0] / temperature;
2346  specific_heat =
2347  _Rw * (-tau * tau * d2gamma1_dtau2(pi, tau) +
2348  Utility::pow<2>(dgamma1_dpi(pi, tau) - tau * d2gamma1_dpitau(pi, tau)) /
2349  d2gamma1_dpi2(pi, tau));
2350  break;
2351 
2352  case 2:
2353  pi = pressure / _p_star[1];
2354  tau = _T_star[1] / temperature;
2355  specific_heat =
2356  _Rw * (-tau * tau * d2gamma2_dtau2(pi, tau) +
2357  Utility::pow<2>(dgamma2_dpi(pi, tau) - tau * d2gamma2_dpitau(pi, tau)) /
2358  d2gamma2_dpi2(pi, tau));
2359  break;
2360 
2361  case 3:
2362  {
2363  // Calculate density first, then use that in Helmholtz free energy
2364  T density3 = densityRegion3(pressure, temperature);
2365  delta = density3 / _rho_critical;
2366  tau = _T_star[2] / temperature;
2367  specific_heat = _Rw * (-tau * tau * d2phi3_dtau2(delta, tau));
2368  break;
2369  }
2370 
2371  case 5:
2372  pi = pressure / _p_star[4];
2373  tau = _T_star[4] / temperature;
2374  specific_heat =
2375  _Rw * (-tau * tau * d2gamma5_dtau2(pi, tau) +
2376  Utility::pow<2>(dgamma5_dpi(pi, tau) - tau * d2gamma5_dpitau(pi, tau)) /
2377  d2gamma5_dpi2(pi, tau));
2378  break;
2379 
2380  default:
2381  mooseError("inRegion() has given an incorrect region");
2382  }
2383  return specific_heat;
2384 }
2385 
2386 template <typename T>
2387 T
2389 {
2390  using std::sqrt, std::exp, std::abs, std::pow;
2391  // Scale the density and temperature. Note that the scales are slightly
2392  // different to the critical values used in IAPWS-IF97
2393  T Tbar = temperature / 647.26;
2394  T rhobar = density / 317.7;
2395 
2396  // Ideal gas component
2397  T sum0 = 0.0;
2398 
2399  for (std::size_t i = 0; i < _k_a.size(); ++i)
2400  sum0 += _k_a[i] * MathUtils::pow(Tbar, i);
2401 
2402  T lambda0 = sqrt(Tbar) * sum0;
2403 
2404  // The contribution due to finite density
2405  T lambda1 =
2406  -0.39707 + 0.400302 * rhobar + 1.06 * exp(-0.171587 * Utility::pow<2>(rhobar + 2.392190));
2407 
2408  // Critical enhancement
2409  T DeltaT = abs(Tbar - 1.0) + 0.00308976;
2410  T Q = 2.0 + 0.0822994 / pow(DeltaT, 0.6);
2411  T S = (Tbar >= 1.0 ? 1.0 / DeltaT : 10.0932 / pow(DeltaT, 0.6));
2412 
2413  T lambda2 =
2414  (0.0701309 / Utility::pow<10>(Tbar) + 0.011852) * pow(rhobar, 1.8) *
2415  exp(0.642857 * (1.0 - pow(rhobar, 2.8))) +
2416  0.00169937 * S * pow(rhobar, Q) * exp((Q / (1.0 + Q)) * (1.0 - pow(rhobar, 1.0 + Q))) -
2417  1.02 * exp(-4.11717 * pow(Tbar, 1.5) - 6.17937 / Utility::pow<5>(rhobar));
2418 
2419  return lambda0 + lambda1 + lambda2;
2420 }
2421 
2422 template <typename T>
2423 T
2425 {
2426  T rho = this->rho_from_p_T(pressure, temperature);
2427  return this->k_from_rho_T_template(rho, temperature);
2428 }
2429 
2430 template <typename T>
2431 T
2432 Water97FluidProperties::k_from_v_e_template(const T & v, const T & e) const
2433 {
2434  const auto [p, temperature] = p_T_from_v_e(v, e);
2436 }
2437 
2438 template <typename T>
2439 T
2440 Water97FluidProperties::p_from_v_e_template(const T & v, const T & e) const
2441 {
2442  const auto rho = 1 / v;
2443  // For NewtonSolve
2444  // y = e
2445  // x = rho
2446  // z = p
2447  auto lambda = [this](const T & rho, const T & p, T & e, T & de_drho, T & de_dp)
2448  { e_from_p_rho(p, rho, e, de_dp, de_drho); };
2450  rho, e, _p_initial_guess, _tolerance, lambda, name() + "::p_from_v_e")
2451  .first;
2452 }
2453 
2454 template <typename T>
2455 T
2456 Water97FluidProperties::e_from_p_rho_template(const T & p, const T & rho) const
2457 {
2458  const auto temperature = T_drhodT_from_p_rho(p, rho).first;
2460 }
2461 
2462 template <typename T>
2463 void
2465  const T & p, const T & rho, T & e, T & de_dp, T & de_drho) const
2466 {
2467  auto functor = [this](const auto & p, const auto & rho)
2468  { return this->e_from_p_rho_template(p, rho); };
2469 
2470  xyDerivatives(p, rho, e, de_dp, de_drho, functor);
2471 }
2472 
2473 template <typename T>
2474 T
2476 {
2477  using std::sqrt, std::exp;
2478 
2479  const T mu_star = 1.e-6;
2480  const T rhobar = density / _rho_critical;
2481  const T Tbar = temperature / _T_critical;
2482 
2483  // Viscosity in limit of zero density
2484  T sum0 = 0.0;
2485  for (std::size_t i = 0; i < _mu_H0.size(); ++i)
2486  sum0 += _mu_H0[i] / MathUtils::pow(Tbar, i);
2487 
2488  const T mu0 = 100.0 * sqrt(Tbar) / sum0;
2489 
2490  // Residual component due to finite density
2491  T sum1 = 0.0;
2492  for (unsigned int i = 0; i < 6; ++i)
2493  {
2494  const T fact = MathUtils::pow(1.0 / Tbar - 1.0, i);
2495  for (unsigned int j = 0; j < 7; ++j)
2496  sum1 += fact * _mu_Hij[i][j] * MathUtils::pow(rhobar - 1.0, j);
2497  }
2498 
2499  const T mu1 = exp(rhobar * sum1);
2500 
2501  // The water viscosity (in Pa.s) is then given by
2502  return mu_star * mu0 * mu1;
2503 }
2504 
2505 template <typename T>
2506 std::pair<T, T>
2507 Water97FluidProperties::p_T_from_v_h(const T & v, const T & h) const
2508 {
2510  bool conversion_succeeded;
2511  p_T_from_v_h(
2512  v, h, _p_initial_guess, _T_initial_guess, pressure, temperature, conversion_succeeded);
2513  return {std::move(pressure), std::move(temperature)};
2514 }
2515 
2516 template <typename T>
2517 T
2519 {
2520  T enthalpy, pi, tau, delta;
2521 
2522  // Determine which region the point is in
2523  unsigned int region =
2525  switch (region)
2526  {
2527  case 1:
2528  pi = pressure / _p_star[0];
2529  tau = _T_star[0] / temperature;
2530  enthalpy = _Rw * _T_star[0] * dgamma1_dtau(pi, tau);
2531  break;
2532 
2533  case 2:
2534  pi = pressure / _p_star[1];
2535  tau = _T_star[1] / temperature;
2536  enthalpy = _Rw * _T_star[1] * dgamma2_dtau(pi, tau);
2537  break;
2538 
2539  case 3:
2540  {
2541  // Calculate density first, then use that in Helmholtz free energy
2542  T density3 = densityRegion3(pressure, temperature);
2543  delta = density3 / _rho_critical;
2544  tau = _T_star[2] / temperature;
2545  enthalpy =
2546  _Rw * temperature * (tau * dphi3_dtau(delta, tau) + delta * dphi3_ddelta(delta, tau));
2547  break;
2548  }
2549 
2550  case 5:
2551  pi = pressure / _p_star[4];
2552  tau = _T_star[4] / temperature;
2553  enthalpy = _Rw * _T_star[4] * dgamma5_dtau(pi, tau);
2554  break;
2555 
2556  default:
2557  mooseError("Water97FluidProperties::inRegion has given an incorrect region");
2558  }
2559  return enthalpy;
2560 }
2561 
2562 template <typename T>
2563 void
2565  const T & pressure, const T & temperature, T & h, T & dh_dp, T & dh_dT) const
2566 {
2567  auto functor = [this](const auto & pressure, const auto & temperature)
2568  { return this->h_from_p_T_template(pressure, temperature); };
2569 
2570  xyDerivatives(pressure, temperature, h, dh_dp, dh_dT, functor);
2571 }
2572 
2573 template <typename T>
2574 T
2576 {
2577  T density, pi, tau;
2578 
2579  // Determine which region the point is in
2580  unsigned int region =
2582 
2583  switch (region)
2584  {
2585  case 1:
2586  pi = pressure / _p_star[0];
2587  tau = _T_star[0] / temperature;
2588  density = pressure / (pi * _Rw * temperature * dgamma1_dpi(pi, tau));
2589  break;
2590 
2591  case 2:
2592  pi = pressure / _p_star[1];
2593  tau = _T_star[1] / temperature;
2594  density = pressure / (pi * _Rw * temperature * dgamma2_dpi(pi, tau));
2595  break;
2596 
2597  case 3:
2599  break;
2600 
2601  case 5:
2602  pi = pressure / _p_star[4];
2603  tau = _T_star[4] / temperature;
2604  density = pressure / (pi * _Rw * temperature * dgamma5_dpi(pi, tau));
2605  break;
2606 
2607  default:
2608  mooseError("inRegion() has given an incorrect region");
2609  }
2610  return density;
2611 }
2612 
2613 template <typename T>
2614 T
2616 {
2617  T rho = this->rho_from_p_T_template(pressure, temperature);
2618  return this->mu_from_rho_T_template(rho, temperature);
2619 }
2620 
2621 template <typename T>
2622 T
2623 Water97FluidProperties::s_from_p_T_template(const T & pressure, const T & temperature) const
2624 {
2625  T entropy, pi, tau, density3, delta;
2626 
2627  // Determine which region the point is in
2628  unsigned int region =
2630  switch (region)
2631  {
2632  case 1:
2633  pi = pressure / _p_star[0];
2634  tau = _T_star[0] / temperature;
2635  entropy = _Rw * (tau * dgamma1_dtau(pi, tau) - gamma1(pi, tau));
2636  break;
2637 
2638  case 2:
2639  pi = pressure / _p_star[1];
2640  tau = _T_star[1] / temperature;
2641  entropy = _Rw * (tau * dgamma2_dtau(pi, tau) - gamma2(pi, tau));
2642  break;
2643 
2644  case 3:
2645  // Calculate density first, then use that in Helmholtz free energy
2646  density3 = densityRegion3(pressure, temperature);
2647  delta = density3 / _rho_critical;
2648  tau = _T_star[2] / temperature;
2649  entropy = _Rw * (tau * dphi3_dtau(delta, tau) - phi3(delta, tau));
2650  break;
2651 
2652  case 5:
2653  pi = pressure / _p_star[4];
2654  tau = _T_star[4] / temperature;
2655  entropy = _Rw * (tau * dgamma5_dtau(pi, tau) - gamma5(pi, tau));
2656  break;
2657 
2658  default:
2659  mooseError("inRegion() has given an incorrect region");
2660  }
2661  return entropy;
2662 }
2663 
2664 template <typename T>
2665 void
2666 Water97FluidProperties::s_from_p_T_template(
2667  const T & pressure, const T & temperature, T & s, T & ds_dp, T & ds_dT) const
2668 {
2669  auto functor = [this](const auto & pressure, const auto & temperature)
2670  { return this->s_from_p_T_template(pressure, temperature); };
2671 
2672  xyDerivatives(pressure, temperature, s, ds_dp, ds_dT, functor);
2673 }
const std::array< Real, 43 > _n2
T v_from_p_T_template(const T &pressure, const T &temperature) const
const Real _T_critical
Critical temperature (K)
T dgamma2_dtau(const T &pi, const T &tau) const
Derivative of Gibbs free energy in Region 2 wrt tau.
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
virtual Real mu_from_p_T(Real pressure, Real temperature) const override
T gamma2(const T &pi, const T &tau) const
Gibbs free energy in Region 2 - superheated steam.
T dphi3_dtau(const T &delta, const T &tau) const
Derivative of Helmholtz free energy in Region 3 wrt tau.
static InputParameters validParams()
virtual void rho_mu_from_p_T(Real pressure, Real temperature, Real &rho, Real &mu) const override
Combined methods.
virtual Real triplePointTemperature() const override
Triple point temperature.
void vaporPressureTemplate(const T &temperature, T &psat, T &dpsat_dT) const
ADReal temperature_from_ph3b(const ADReal &pressure, const ADReal &enthalpy) const
Backwards equation T(p, h) in Region 3b Eq.
const std::array< int, 6 > _I5
ADReal temperature_from_ph3a(const ADReal &pressure, const ADReal &enthalpy) const
Backwards equation T(p, h) in Region 3a Eq.
unsigned int inRegion(Real pressure, Real temperature) const
Determines the phase region that the given pressure and temperature values lie in.
const std::array< Real, 5 > _n23
Constants for the boundary between regions 2 and 3.
T tempXY(const T &pressure, subregionEnum xy) const
Boundaries between subregions in region 3.
T d2phi3_ddelta2(const T &delta, const T &tau) const
Second derivative of Helmholtz free energy in Region 3 wrt delta.
std::pair< T, T > rho_T_from_v_e(const T &v, const T &e) const
Computes the density (first member of the pair) and temperature (second member of the pair) as functi...
Real vaporTemperature(Real pressure) const override
Saturation temperature as a function of pressure.
void p_T_from_v_h(const T &v, const T &h, Real p0, Real T0, T &pressure, T &temperature, bool &conversion_succeeded) const
Determines (p,T) from (v,h) using Newton Solve in 2D Useful for conversion between different sets of ...
Real p_from_v_e(Real v, Real e) const override
ADReal temperature_from_ph2b(const ADReal &pressure, const ADReal &enthalpy) const
Backwards equation T(p, h) in Region 2b Eq.
auto exp(const T &)
const Real _p_critical
Critical pressure (Pa)
const std::array< int, 33 > _Jph3b
const std::array< std::array< Real, 8 >, 26 > _par3
Real T_from_v_e(Real v, Real e) const override
const std::array< int, 20 > _Iph1
virtual Real e_from_p_T(Real pressure, Real temperature) const override
Real b2bc(Real pressure) const
Boundary between subregions b and c in region 2.
T subregionVolume(const T &pi, const T &theta, Real a, Real b, Real c, Real d, Real e, unsigned int sid) const
Specific volume in all subregions of region 3 EXCEPT subregion n (13).
virtual std::string fluidName() const override
Fluid name.
const std::array< int, 23 > _Jph2c
const std::array< Real, 36 > _nph2a
T k_from_v_e_template(const T &v, const T &e) const
T e_from_p_T_template(const T &pressure, const T &temperature) const
propfuncWithDefinitionOverride(s, p, T)
virtual Real criticalDensity() const override
Critical density.
T k_from_p_T_template(const T &pressure, const T &temperature) const
Real b23T(Real pressure) const
Auxillary equation for the boundary between regions 2 and 3.
const InputParameters & parameters() const
T d2gamma5_dpitau(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 5 wrt pi and tau.
Water97FluidProperties(const InputParameters &parameters)
virtual ADReal cv_from_v_e(const ADReal &v, const ADReal &e) const override
T dgamma5_dtau(const T &pi, const T &tau) const
Derivative of Gibbs free energy in Region 5 wrt tau.
T c_from_p_T_template(const T &pressure, const T &temperature) const
T d2gamma2_dpitau(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 2 wrt pi and tau.
T rho_from_p_T_template(const T &pressure, const T &temperature) const
unsigned int subregion3(Real pressure, Real temperature) const
Provides the correct subregion index for a (P,T) point in region 3.
const std::array< int, 36 > _Jph2a
static const std::string density
Definition: NS.h:33
T d2gamma5_dtau2(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 5 wrt tau.
auto raw_value(const Eigen::Map< T > &in)
int delta(unsigned int i, unsigned int j)
Delta function, which returns zero if $i j$ and unity if $i=j$.
const std::array< std::array< Real, 7 >, 6 > _mu_Hij
ADReal vaporTemperature_ad(const ADReal &pressure) const
AD version of saturation temperature as a function of pressure (used internally)
const std::vector< std::vector< Real > > _n3s
Constants for all 26 subregions in region 3.
T h_from_p_T_template(const T &pressure, const T &temperature) const
const std::array< Real, 33 > _nph3b
static const std::string temperature
Definition: NS.h:59
const std::array< Real, 38 > _nph2b
std::pair< T, T > p_T_from_v_h(const T &v, const T &h) const
Computes the pressure (first member of the pair) and temperature (second member of the pair) as funct...
virtual Real h_from_p_T(Real pressure, Real temperature) const override
virtual Real criticalTemperature() const override
Critical temperature.
const std::array< int, 33 > _Iph3b
virtual Real k_from_p_T(Real pressure, Real temperature) const override
ADReal mu_from_v_e(const ADReal &v, const ADReal &e) const override
DualNumber< Real, DNDerivativeType, true > ADReal
const std::array< int, 34 > _J1
unsigned int subregion2ph(Real pressure, Real enthalpy) const
Provides the correct subregion index for a (P,h) point in region 2.
const std::array< int, 9 > _J02
const Real _T_triple
Triple point temperature (K)
const std::array< int, 43 > _J2
const std::array< int, 6 > _J5
std::pair< T, T > p_T_from_v_e(const T &v, const T &e) const
Computes the pressure (first member of the pair) and temperature (second member of the pair) as funct...
const std::array< int, 23 > _Iph2c
std::pair< T, T > NewtonSolve(const T &x, const T &y, const Real z_initial_guess, const Real tolerance, const Functor &func, const std::string &caller_name, const unsigned int max_its=100)
NewtonSolve does a 1D Newton Solve to solve the equation y = f(x, z) for variable z...
ADReal temperature_from_ph2a(const ADReal &pressure, const ADReal &enthalpy) const
Backwards equation T(p, h) in Region 2a Eq.
const std::array< Real, 5 > _p_star
Pressure scale for each region.
static void xyDerivatives(const T x, const T &y, T &z, T &dz_dx, T &dz_dy, const Functor &z_from_x_y)
Computes the dependent variable z and its derivatives with respect to the independent variables x and...
subregionEnum
Enum of subregion ids for region 3.
T cv_from_p_T_template(const T &pressure, const T &temperature) const
const std::array< int, 20 > _Jph1
const std::array< Real, 34 > _n1
Reference constants used in to calculate thermophysical properties of water.
const std::array< Real, 5 > _T_star
Temperature scale for each region.
const Real _tolerance
Newton&#39;s method may be used to convert between variable sets.
const std::vector< std::vector< int > > _I3s
e e e e s T T T T T rho v v T e h
const std::string & name() const
const std::array< int, 36 > _Iph2a
T d2gamma1_dpitau(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 1 wrt pi and tau.
Real k_from_v_e(Real v, Real e) const override
static const std::string S
Definition: NS.h:163
const std::array< int, 38 > _Jph2b
const std::array< int, 31 > _Jph3a
virtual Real s_from_h_p(Real enthalpy, Real pressure) const override
T dphi3_ddelta(const T &delta, const T &tau) const
Derivative of Helmholtz free energy in Region 3 wrt delta.
T d2gamma2_dpi2(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 2 wrt pi.
static const std::string mu
Definition: NS.h:123
T gamma5(const T &pi, const T &tau) const
Gibbs free energy in Region 5.
T d2gamma1_dtau2(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 1 wrt tau.
ADReal T_from_p_h_ad(const ADReal &pressure, const ADReal &enthalpy) const
AD version of backwards equation T(p, h) (used internally) From Revised Release on the IAPWS Industri...
const std::array< int, 31 > _Iph3a
virtual Real v_from_p_T(Real pressure, Real temperature) const override
T dgamma2_dpi(const T &pi, const T &tau) const
Derivative of Gibbs free energy in Region 2 wrt pi.
virtual Real triplePointPressure() const override
Triple point pressure.
T dgamma1_dtau(const T &pi, const T &tau) const
Derivative of Gibbs free energy in Region 1 wrt tau.
Common class for single phase fluid properties.
const std::array< Real, 20 > _nph1
ADReal temperature_from_ph1(const ADReal &pressure, const ADReal &enthalpy) const
Backwards equation T(p, h) in Region 1 Eq.
const std::array< int, 40 > _J3
T d2gamma2_dtau2(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 2 wrt tau.
ADReal e_from_v_h(const ADReal &v, const ADReal &h) const override
const std::array< int, 34 > _I1
const std::array< int, 6 > _J05
Constants for region 5.
auto log(const T &)
unsigned int subregion3ph(Real pressure, Real enthalpy) const
Provides the correct subregion index for a (P,h) point in region 3.
Real volume(const MeshBase &mesh, unsigned int dim=libMesh::invalid_uint)
T dgamma1_dpi(const T &pi, const T &tau) const
Derivative of Gibbs free energy in Region 1 wrt pi.
T dgamma5_dpi(const T &pi, const T &tau) const
Derivative of Gibbs free energy in Region 5 wrt pi.
virtual Real T_from_p_h(Real pressure, Real enthalpy) const override
Backwards equation T(p, h) From Revised Release on the IAPWS Industrial Formulation 1997 for the Ther...
T mu_from_rho_T_template(const T &density, const T &temperature) const
T d2phi3_dtau2(const T &delta, const T &tau) const
Second derivative of Helmholtz free energy in Region 3 wrt tau.
virtual Real c_from_p_T(Real pressure, Real temperature) const override
const std::vector< std::vector< Real > > _tempXY_n
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
T p_from_v_e_template(const T &v, const T &e) const
const std::array< Real, 31 > _nph3a
T k_from_rho_T_template(const T &density, const T &temperature) const
Real b23p(Real temperature) const
Auxillary equation for the boundary between regions 2 and 3.
const std::array< Real, 6 > _n5
const Real _rho_critical
Critical density (kg/m^3)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string v
Definition: NS.h:84
T phi3(const T &delta, const T &tau) const
Helmholtz free energy in Region 3.
const std::array< int, 38 > _Iph2b
int N
T densityRegion3(const T &pressure, const T &temperature) const
Density function for Region 3 - supercritical water and steam.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sqrt(_arg)) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tanh
virtual Real e_from_p_rho(Real p, Real rho) const override
const std::array< Real, 6 > _n05
virtual Real criticalPressure() const override
Critical pressure.
Real b3ab(Real pressure) const
Boundary between subregions a and b in region 3.
Water (H2O) fluid properties as a function of pressure (Pa) and temperature (K) from IAPWS-IF97: Revi...
static const std::string pressure
Definition: NS.h:56
virtual Real vaporPressure(Real temperature) const override
Vapor pressure.
const Real _p_initial_guess
Initial guess for pressure (or pressure used to compute the initial guess)
T d2gamma5_dpi2(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 5 wrt pi.
void mooseError(Args &&... args) const
T gamma1(const T &pi, const T &tau) const
Gibbs free energy in Region 1 - single phase liquid region.
Real henryConstant(Real temperature, const std::vector< Real > &coeffs) const
IAPWS formulation of Henry&#39;s law constant for dissolution in water From Guidelines on the Henry&#39;s con...
unsigned int inRegionPH(Real pressure, Real enthalpy) const
Determines the phase region that the given pressure and enthaply values lie in.
const std::vector< std::vector< int > > _J3s
static const std::string internal_energy
Definition: NS.h:61
virtual Real k_from_rho_T(Real density, Real temperature) const override
T cp_from_p_T_template(const T &pressure, const T &temperature) const
const Real _Rw
Specific gas constant for H2O (universal gas constant / molar mass of water - J/kg/K) ...
virtual Real rho_from_p_T(Real pressure, Real temperature) const override
static const std::string specific_heat
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
const std::array< int, 40 > _I3
const Real _p_triple
Triple point pressure (Pa)
T mu_from_p_T_template(const T &pressure, const T &temperature) const
const std::array< Real, 9 > _n02
Constants for region 2.
T e_from_p_rho_template(const T &p, const T &rho) const
T pow(T x, int e)
virtual ADReal c_from_v_e(const ADReal &v, const ADReal &e) const override
virtual Real cv_from_p_T(Real pressure, Real temperature) const override
ADReal temperature_from_ph2c(const ADReal &pressure, const ADReal &enthalpy) const
Backwards equation T(p, h) in Region 2c Eq.
T d2gamma1_dpi2(const T &pi, const T &tau) const
Second derivative of Gibbs free energy in Region 1 wrt pi.
virtual Real mu_from_rho_T(Real density, Real temperature) const override
std::pair< T, T > T_drhodT_from_p_rho(const T &p, const T &rho) const
Computes the temperature (first member of the pair) and the derivative of density (second member of t...
const std::array< Real, 23 > _nph2c
const std::array< int, 43 > _I2
const Real _T_initial_guess
Initial guess for temperature (or temperature used to compute the initial guess)
virtual Real cp_from_p_T(Real pressure, Real temperature) const override
const Real _Mh2o
Water molar mass (kg/mol)
MooseUnits pow(const MooseUnits &, int)
const std::array< Real, 40 > _n3
Constants for region 3.
static const std::string k
Definition: NS.h:130
std::array< Real, 4 > _k_a
Constants for thermal conductivity.
const std::array< unsigned int, 26 > _par3N
virtual ADReal cp_from_v_e(const ADReal &v, const ADReal &e) const override
T d2phi3_ddeltatau(const T &delta, const T &tau) const
Second derivative of Helmholtz free energy in Region 3 wrt delta and tau.
virtual Real molarMass() const override
Molar mass [kg/mol].
const Real pi
const std::vector< std::vector< int > > _tempXY_I
Constnats for the tempXY() method.
const std::array< Real, 10 > _n4
Constants for region 4 (the saturation curve up to the critical point)
const std::array< Real, 4 > _mu_H0
Constants from Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance...