LCOV - code coverage report
Current view: top level - include/fluidproperties - Water97FluidProperties.h (source / functions) Hit Total Coverage
Test: idaholab/moose fluid_properties: #31405 (292dce) with base fef103 Lines: 542 589 92.0 %
Date: 2025-09-04 07:53:14 Functions: 115 155 74.2 %
Legend: Lines: hit not hit

          Line data    Source code
       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             : 
      12             : #include "SinglePhaseFluidProperties.h"
      13             : #include "NewtonInversion.h"
      14             : #include <array>
      15             : 
      16             : #pragma GCC diagnostic push
      17             : #pragma GCC diagnostic ignored "-Woverloaded-virtual"
      18             : 
      19             : /**
      20             :  * Water (H2O) fluid properties as a function of pressure (Pa)
      21             :  * and temperature (K) from IAPWS-IF97:
      22             :  * Revised Release on the IAPWS Industrial Formulation 1997 for
      23             :  * the Thermodynamic Properties of Water and Steam.
      24             :  *
      25             :  * To avoid iteration in Region 3, the backwards equations from:
      26             :  * Revised Supplementary Release on Backward Equations for Specific Volume
      27             :  * as a Function of Pressure and Temperature v(p,T) for Region 3 of the
      28             :  * IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water
      29             :  * and Steam are implemented.
      30             :  *
      31             :  * Water viscosity as a function of density (kg/m^3) and temperature (K) from:
      32             :  * Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance.
      33             :  *
      34             :  * Thermal conductivity as a function of density (kg/m^3) and temperature (K) from:
      35             :  * Revised Release on the IAPS Formulation 1985 for the Thermal Conductivity of
      36             :  * Ordinary Water Substance
      37             :  */
      38             : class Water97FluidProperties : public SinglePhaseFluidProperties
      39             : {
      40             : public:
      41             :   static InputParameters validParams();
      42             : 
      43             :   Water97FluidProperties(const InputParameters & parameters);
      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>
      81             :   void rho_from_p_T_template(
      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,
     174             :                                Real temperature,
     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             : 
     199             :   propfuncWithDefinitionOverride(s, p, T);
     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             : 
     229             :   /**
     230             :    * Saturation temperature as a function of pressure
     231             :    *
     232             :    * Eq. (31) from Revised Release on the IAPWS Industrial
     233             :    * Formulation 1997 for the Thermodynamic Properties of Water
     234             :    * and Steam
     235             :    *
     236             :    * Valid for 611.213 Pa <= p <= 22.064 MPa
     237             :    *
     238             :    * @param pressure water pressure (Pa)
     239             :    * @return saturation temperature (K)
     240             :    */
     241             :   Real vaporTemperature(Real pressure) const override;
     242             :   virtual void vaporTemperature(Real pressure, Real & Tsat, Real & dTsat_dp) const override;
     243             : 
     244             :   /**
     245             :    * Auxillary equation for the boundary between regions 2 and 3
     246             :    *
     247             :    * Eq. (5) from Revised Release on the IAPWS Industrial
     248             :    * Formulation 1997 for the Thermodynamic Properties of Water
     249             :    * and Steam
     250             :    *
     251             :    * @param temperature water temperature (K)
     252             :    * @return pressure (Pa) on the boundary between region 2 and 3
     253             :    */
     254             :   Real b23p(Real temperature) const;
     255             : 
     256             :   /**
     257             :    * Auxillary equation for the boundary between regions 2 and 3
     258             :    *
     259             :    * Eq. (6) from Revised Release on the IAPWS Industrial
     260             :    * Formulation 1997 for the Thermodynamic Properties of Water
     261             :    * and Steam
     262             :    *
     263             :    * @param pressure water pressure (Pa)
     264             :    * @return temperature (K) on the boundary between regions 2 and 3
     265             :    */
     266             :   Real b23T(Real pressure) const;
     267             : 
     268             :   /**
     269             :    * Determines the phase region that the given pressure and temperature values
     270             :    * lie in
     271             :    *
     272             :    * @param pressure water pressure (Pa)
     273             :    * @param temperature water temperature (K)
     274             :    * @return region phase region index
     275             :    */
     276             :   unsigned int inRegion(Real pressure, Real temperature) const;
     277             : 
     278             :   /**
     279             :    * Provides the correct subregion index for a (P,T) point in
     280             :    * region 3. From Revised Supplementary Release on Backward Equations for
     281             :    * Specific Volume as a Function of Pressure and Temperature v(p,T)
     282             :    * for Region 3 of the IAPWS Industrial Formulation 1997 for the
     283             :    * Thermodynamic Properties of Water and Steam
     284             :    *
     285             :    * @param pressure water pressure (Pa)
     286             :    * @param temperature water temperature (K)
     287             :    * @return subregion index
     288             :    */
     289             :   unsigned int subregion3(Real pressure, Real temperature) const;
     290             : 
     291             :   /**
     292             :    * Specific volume in all subregions of region 3 EXCEPT subregion n (13).
     293             :    *
     294             :    * @param pi scaled water pressure
     295             :    * @param theta scaled water temperature
     296             :    * @param a to e constants
     297             :    * @param sid subregion ID of the subregion
     298             :    * @return volume water specific volume (m^3/kg)
     299             :    */
     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             : 
     310             :   /**
     311             :    * Density function for Region 3 - supercritical water and steam
     312             :    *
     313             :    * To avoid iteration, use the backwards equations for region 3
     314             :    * from Revised Supplementary Release on Backward Equations for
     315             :    * Specific Volume as a Function of Pressure and Temperature v(p,T)
     316             :    * for Region 3 of the IAPWS Industrial Formulation 1997 for the
     317             :    * Thermodynamic Properties of Water and Steam.
     318             :    *
     319             :    * @param pressure water pressure (Pa)
     320             :    * @param temperature water temperature (K)
     321             :    * @return density (kg/m^3) in region 3
     322             :    */
     323             :   template <typename T>
     324             :   T densityRegion3(const T & pressure, const T & temperature) const;
     325             : 
     326             :   /**
     327             :    * Backwards equation T(p, h)
     328             :    * From Revised Release on the IAPWS Industrial Formulation 1997 for the
     329             :    * Thermodynamic Properties of Water and Steam
     330             :    *
     331             :    * @param pressure water pressure (Pa)
     332             :    * @param enthalpy water enthalpy (J/kg)
     333             :    * @return temperature water temperature (K)
     334             :    */
     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             : 
     339             :   /**
     340             :    * Boundary between subregions b and c in region 2.
     341             :    * Equation 21 from Revised Release on the IAPWS Industrial
     342             :    * Formulation 1997 for the Thermodynamic Properties of Water
     343             :    * and Steam.
     344             :    *
     345             :    * @param pressure water pressure (Pa)
     346             :    * @return enthalpy at boundary (J/kg)
     347             :    */
     348             :   Real b2bc(Real pressure) const;
     349             : 
     350             :   /**
     351             :    * Boundary between subregions a and b in region 3.
     352             :    * Equation 1 from Revised Supplementary Release on Backward Equations for
     353             :    * the Functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS
     354             :    * Industrial Formulation 1997 for the Thermodynamic Properties of Water and
     355             :    * Steam
     356             :    *
     357             :    * @param pressure water pressure (Pa)
     358             :    * @return enthalpy at boundary (J/kg)
     359             :    */
     360             :   Real b3ab(Real pressure) const;
     361             : 
     362             :   /**
     363             :    * IAPWS formulation of Henry's law constant for dissolution in water
     364             :    * From Guidelines on the Henry's constant and vapour
     365             :    * liquid distribution constant for gases in H20 and D20 at high
     366             :    * temperatures, IAPWS (2004)
     367             :    * @param T fluid temperature (K)
     368             :    * @param coeffs Henry's constant coefficients of gas
     369             :    * @param[out] Kh Henry's constant
     370             :    * @param[out] dKh_dT derivative of Kh wrt temperature
     371             :    */
     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             : 
     377             :   /**
     378             :    * Computes the pressure (first member of the pair) and temperature (second member of the pair) as
     379             :    * functions of specific volume and specific internal energy
     380             :    */
     381             :   template <typename T>
     382             :   std::pair<T, T> p_T_from_v_e(const T & v, const T & e) const;
     383             : 
     384             :   /**
     385             :    * Computes the density (first member of the pair) and temperature (second member of the pair) as
     386             :    * functions of specific volume and specific internal energy
     387             :    */
     388             :   template <typename T>
     389             :   std::pair<T, T> rho_T_from_v_e(const T & v, const T & e) const;
     390             : 
     391             :   /**
     392             :    * Computes the pressure (first member of the pair) and temperature (second member of the pair) as
     393             :    * functions of specific volume and specific enthalpy
     394             :    */
     395             :   template <typename T>
     396             :   std::pair<T, T> p_T_from_v_h(const T & v, const T & h) const;
     397             :   using SinglePhaseFluidProperties::p_T_from_v_h;
     398             : 
     399             : protected:
     400             :   /**
     401             :    * Gibbs free energy in Region 1 - single phase liquid region
     402             :    *
     403             :    * From Eq. (7) From Revised Release on the IAPWS Industrial
     404             :    * Formulation 1997 for the Thermodynamic Properties of Water
     405             :    * and Steam, IAPWS 2007.
     406             :    *
     407             :    * @param pi reduced pressure (-)
     408             :    * @param tau reduced temperature (-)
     409             :    * @return Gibbs free energy (-)
     410             :    */
     411             :   template <typename T>
     412             :   T gamma1(const T & pi, const T & tau) const;
     413             : 
     414             :   /**
     415             :    * Derivative of Gibbs free energy in Region 1 wrt pi
     416             :    *
     417             :    * @param pi reduced pressure (-)
     418             :    * @param tau reduced temperature (-)
     419             :    * @return derivative of Gibbs free energy wrt pi (-)
     420             :    */
     421             :   template <typename T>
     422             :   T dgamma1_dpi(const T & pi, const T & tau) const;
     423             : 
     424             :   /**
     425             :    * Second derivative of Gibbs free energy in Region 1 wrt pi
     426             :    *
     427             :    * @param pi reduced pressure (-)
     428             :    * @param tau reduced temperature (-)
     429             :    * @return second derivative of Gibbs free energy wrt pi (-)
     430             :    */
     431             :   template <typename T>
     432             :   T d2gamma1_dpi2(const T & pi, const T & tau) const;
     433             : 
     434             :   /**
     435             :    * Derivative of Gibbs free energy in Region 1 wrt tau
     436             :    *
     437             :    * @param pi reduced pressure (-)
     438             :    * @param tau reduced temperature (-)
     439             :    * @return derivative of Gibbs free energy wrt tau (-)
     440             :    */
     441             :   template <typename T>
     442             :   T dgamma1_dtau(const T & pi, const T & tau) const;
     443             : 
     444             :   /**
     445             :    * Second derivative of Gibbs free energy in Region 1 wrt tau
     446             :    *
     447             :    * @param pi reduced pressure (-)
     448             :    * @param tau reduced temperature (-)
     449             :    * @return second derivative of Gibbs free energy wrt tau (-)
     450             :    */
     451             :   template <typename T>
     452             :   T d2gamma1_dtau2(const T & pi, const T & tau) const;
     453             : 
     454             :   /**
     455             :    * Second derivative of Gibbs free energy in Region 1 wrt pi and tau
     456             :    *
     457             :    * @param pi reduced pressure (-)
     458             :    * @param tau reduced temperature (-)
     459             :    * @return second derivative of Gibbs free energy wrt pi and tau (-)
     460             :    */
     461             :   template <typename T>
     462             :   T d2gamma1_dpitau(const T & pi, const T & tau) const;
     463             : 
     464             :   /**
     465             :    * Gibbs free energy in Region 2 - superheated steam
     466             :    *
     467             :    * From Eq. (15) From Revised Release on the IAPWS Industrial
     468             :    * Formulation 1997 for the Thermodynamic Properties of Water
     469             :    * and Steam, IAPWS 2007.
     470             :    *
     471             :    * @param pi reduced pressure (-)
     472             :    * @param tau reduced temperature (-)
     473             :    * @return Gibbs free energy (-)
     474             :    */
     475             :   template <typename T>
     476             :   T gamma2(const T & pi, const T & tau) const;
     477             : 
     478             :   /**
     479             :    * Derivative of Gibbs free energy in Region 2 wrt pi
     480             :    *
     481             :    * @param pi reduced pressure (-)
     482             :    * @param tau reduced temperature (-)
     483             :    * @return derivative of Gibbs free energy wrt pi (-)
     484             :    */
     485             :   template <typename T>
     486             :   T dgamma2_dpi(const T & pi, const T & tau) const;
     487             : 
     488             :   /**
     489             :    * Second derivative of Gibbs free energy in Region 2 wrt pi
     490             :    *
     491             :    * @param pi reduced pressure (-)
     492             :    * @param tau reduced temperature (-)
     493             :    * @return second derivative of Gibbs free energy wrt pi (-)
     494             :    */
     495             :   template <typename T>
     496             :   T d2gamma2_dpi2(const T & pi, const T & tau) const;
     497             : 
     498             :   /**
     499             :    * Derivative of Gibbs free energy in Region 2 wrt tau
     500             :    *
     501             :    * @param pi reduced pressure (-)
     502             :    * @param tau reduced temperature (-)
     503             :    * @return derivative of Gibbs free energy wrt tau (-)
     504             :    */
     505             :   template <typename T>
     506             :   T dgamma2_dtau(const T & pi, const T & tau) const;
     507             : 
     508             :   /**
     509             :    * Second derivative of Gibbs free energy in Region 2 wrt tau
     510             :    *
     511             :    * @param pi reduced pressure (-)
     512             :    * @param tau reduced temperature (-)
     513             :    * @return second derivative of Gibbs free energy wrt tau (-)
     514             :    */
     515             :   template <typename T>
     516             :   T d2gamma2_dtau2(const T & pi, const T & tau) const;
     517             : 
     518             :   /**
     519             :    * Second derivative of Gibbs free energy in Region 2 wrt pi and tau
     520             :    *
     521             :    * @param pi reduced pressure (-)
     522             :    * @param tau reduced temperature (-)
     523             :    * @return second derivative of Gibbs free energy wrt pi and tau (-)
     524             :    */
     525             :   template <typename T>
     526             :   T d2gamma2_dpitau(const T & pi, const T & tau) const;
     527             : 
     528             :   /**
     529             :    * Helmholtz free energy in Region 3
     530             :    *
     531             :    * From Eq. (28) From Revised Release on the IAPWS Industrial
     532             :    * Formulation 1997 for the Thermodynamic Properties of Water
     533             :    * and Steam, IAPWS 2007.
     534             :    *
     535             :    * @param delta reduced density (-)
     536             :    * @param tau reduced temperature (-)
     537             :    * @return Helmholtz free energy (-)
     538             :    */
     539             :   template <typename T>
     540             :   T phi3(const T & delta, const T & tau) const;
     541             : 
     542             :   /**
     543             :    * Derivative of Helmholtz free energy in Region 3 wrt delta
     544             :    *
     545             :    * @param delta reduced density (-)
     546             :    * @param tau reduced temperature (-)
     547             :    * @return derivative of Helmholtz free energy wrt delta (-)
     548             :    */
     549             :   template <typename T>
     550             :   T dphi3_ddelta(const T & delta, const T & tau) const;
     551             : 
     552             :   /**
     553             :    * Second derivative of Helmholtz free energy in Region 3 wrt delta
     554             :    *
     555             :    * @param delta reduced density (-)
     556             :    * @param tau reduced temperature (-)
     557             :    * @return second derivative of Helmholtz free energy wrt delta (-)
     558             :    */
     559             :   template <typename T>
     560             :   T d2phi3_ddelta2(const T & delta, const T & tau) const;
     561             : 
     562             :   /**
     563             :    * Derivative of Helmholtz free energy in Region 3 wrt tau
     564             :    *
     565             :    * @param delta reduced density (-)
     566             :    * @param tau reduced temperature (-)
     567             :    * @return derivative of Helmholtz free energy wrt tau (-)
     568             :    */
     569             :   template <typename T>
     570             :   T dphi3_dtau(const T & delta, const T & tau) const;
     571             : 
     572             :   /**
     573             :    * Second derivative of Helmholtz free energy in Region 3 wrt tau
     574             :    *
     575             :    * @param delta reduced density (-)
     576             :    * @param tau reduced temperature (-)
     577             :    * @return second derivative of Helmholtz free energy wrt tau (-)
     578             :    */
     579             :   template <typename T>
     580             :   T d2phi3_dtau2(const T & delta, const T & tau) const;
     581             : 
     582             :   /**
     583             :    * Second derivative of Helmholtz free energy in Region 3 wrt delta and tau
     584             :    *
     585             :    * @param delta reduced density (-)
     586             :    * @param tau reduced temperature (-)
     587             :    * @return second derivative of Helmholtz free energy wrt delta and tau (-)
     588             :    */
     589             :   template <typename T>
     590             :   T d2phi3_ddeltatau(const T & delta, const T & tau) const;
     591             : 
     592             :   /**
     593             :    * Gibbs free energy in Region 5
     594             :    *
     595             :    * From Eq. (32) From Revised Release on the IAPWS Industrial
     596             :    * Formulation 1997 for the Thermodynamic Properties of Water
     597             :    * and Steam, IAPWS 2007.
     598             :    *
     599             :    * @param pi reduced pressure (-)
     600             :    * @param tau reduced temperature (-)
     601             :    * @return Gibbs free energy (-)
     602             :    */
     603             :   template <typename T>
     604             :   T gamma5(const T & pi, const T & tau) const;
     605             : 
     606             :   /**
     607             :    * Derivative of Gibbs free energy in Region 5 wrt pi
     608             :    *
     609             :    * @param pi reduced pressure (-)
     610             :    * @param tau reduced temperature (-)
     611             :    * @return derivative of Gibbs free energy wrt pi (-)
     612             :    */
     613             :   template <typename T>
     614             :   T dgamma5_dpi(const T & pi, const T & tau) const;
     615             : 
     616             :   /**
     617             :    * Second derivative of Gibbs free energy in Region 5 wrt pi
     618             :    *
     619             :    * @param pi reduced pressure (-)
     620             :    * @param tau reduced temperature (-)
     621             :    * @return second derivative of Gibbs free energy wrt pi (-)
     622             :    */
     623             :   template <typename T>
     624             :   T d2gamma5_dpi2(const T & pi, const T & tau) const;
     625             : 
     626             :   /**
     627             :    * Derivative of Gibbs free energy in Region 5 wrt tau
     628             :    *
     629             :    * @param pi reduced pressure (-)
     630             :    * @param tau reduced temperature (-)
     631             :    * @return derivative of Gibbs free energy wrt tau (-)
     632             :    */
     633             :   template <typename T>
     634             :   T dgamma5_dtau(const T & pi, const T & tau) const;
     635             : 
     636             :   /**
     637             :    * Second derivative of Gibbs free energy in Region 5 wrt tau
     638             :    *
     639             :    * @param pi reduced pressure (-)
     640             :    * @param tau reduced temperature (-)
     641             :    * @return second derivative of Gibbs free energy wrt tau (-)
     642             :    */
     643             :   template <typename T>
     644             :   T d2gamma5_dtau2(const T & pi, const T & tau) const;
     645             : 
     646             :   /**
     647             :    * Second derivative of Gibbs free energy in Region 5 wrt pi and tau
     648             :    *
     649             :    * @param pi reduced pressure (-)
     650             :    * @param tau reduced temperature (-)
     651             :    * @return second derivative of Gibbs free energy wrt pi and tau (-)
     652             :    */
     653             :   template <typename T>
     654             :   T d2gamma5_dpitau(const T & pi, const T & tau) const;
     655             : 
     656             :   /// Enum of subregion ids for region 3
     657             :   enum subregionEnum
     658             :   {
     659             :     AB,
     660             :     CD,
     661             :     GH,
     662             :     IJ,
     663             :     JK,
     664             :     MN,
     665             :     OP,
     666             :     QU,
     667             :     RX,
     668             :     UV,
     669             :     WX,
     670             :     EF
     671             :   };
     672             : 
     673             :   /**
     674             :    * Boundaries between subregions in region 3.
     675             :    * From Revised Supplementary Release on Backward Equations for
     676             :    * Specific Volume as a Function of Pressure and Temperature v(p,T)
     677             :    * for Region 3 of the IAPWS Industrial Formulation 1997 for the
     678             :    * Thermodynamic Properties of Water and Steam
     679             :    *
     680             :    * @param pressure water pressure (Pa)
     681             :    * @param xy string to select the boundary between two subregions
     682             :    * @return temperature (K) along the boundary
     683             :    */
     684             :   template <typename T>
     685             :   T tempXY(const T & pressure, subregionEnum xy) const;
     686             : 
     687             :   /**
     688             :    * Determines the phase region that the given pressure and enthaply values
     689             :    * lie in.
     690             :    *
     691             :    * @param pressure water pressure (Pa)
     692             :    * @param enthalpy water enthalpy (J/kg)
     693             :    * @return region phase region index
     694             :    */
     695             :   unsigned int inRegionPH(Real pressure, Real enthalpy) const;
     696             : 
     697             :   /**
     698             :    * Provides the correct subregion index for a (P,h) point in
     699             :    * region 2. From Revised Release on the IAPWS Industrial
     700             :    * Formulation 1997 for the Thermodynamic Properties of Water
     701             :    * and Steam
     702             :    *
     703             :    * @param pressure water pressure (Pa)
     704             :    * @param enthalpy water enthalpy (J/kg)
     705             :    * @return subregion index
     706             :    */
     707             :   unsigned int subregion2ph(Real pressure, Real enthalpy) const;
     708             : 
     709             :   /**
     710             :    * Provides the correct subregion index for a (P,h) point in
     711             :    * region 3. From Revised Supplementary Release on Backward Equations for
     712             :    * the Functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS
     713             :    * Industrial Formulation 1997 for the Thermodynamic Properties of Water and
     714             :    * Steam
     715             :    *
     716             :    * @param pressure water pressure (Pa)
     717             :    * @param enthalpy water enthalpy (J/kg)
     718             :    * @return subregion index
     719             :    */
     720             :   unsigned int subregion3ph(Real pressure, Real enthalpy) const;
     721             : 
     722             :   /**
     723             :    * AD version of backwards equation T(p, h) (used internally)
     724             :    * From Revised Release on the IAPWS Industrial Formulation 1997 for the
     725             :    * Thermodynamic Properties of Water and Steam
     726             :    *
     727             :    * @param pressure water pressure (Pa)
     728             :    * @param enthalpy water enthalpy (J/kg)
     729             :    * @return temperature water temperature (K)
     730             :    */
     731             :   ADReal T_from_p_h_ad(const ADReal & pressure, const ADReal & enthalpy) const;
     732             : 
     733             :   /**
     734             :    * Backwards equation T(p, h) in Region 1
     735             :    * Eq. (11) from Revised Release on the IAPWS Industrial
     736             :    * Formulation 1997 for the Thermodynamic Properties of Water
     737             :    * and Steam
     738             :    *
     739             :    * @param pressure water pressure (Pa)
     740             :    * @param enthalpy water enthalpy (J/kg)
     741             :    * @return temperature water temperature (K)
     742             :    */
     743             :   ADReal temperature_from_ph1(const ADReal & pressure, const ADReal & enthalpy) const;
     744             : 
     745             :   /**
     746             :    * Backwards equation T(p, h) in Region 2a
     747             :    * Eq. (22) from Revised Release on the IAPWS Industrial
     748             :    * Formulation 1997 for the Thermodynamic Properties of Water
     749             :    * and Steam
     750             :    *
     751             :    * @param pressure water pressure (Pa)
     752             :    * @param enthalpy water enthalpy (J/kg)
     753             :    * @return temperature water temperature (K)
     754             :    */
     755             :   ADReal temperature_from_ph2a(const ADReal & pressure, const ADReal & enthalpy) const;
     756             : 
     757             :   /**
     758             :    * Backwards equation T(p, h) in Region 2b
     759             :    * Eq. (23) from Revised Release on the IAPWS Industrial
     760             :    * Formulation 1997 for the Thermodynamic Properties of Water
     761             :    * and Steam
     762             :    *
     763             :    * @param pressure water pressure (Pa)
     764             :    * @param enthalpy water enthalpy (J/kg)
     765             :    * @return temperature water temperature (K)
     766             :    */
     767             :   ADReal temperature_from_ph2b(const ADReal & pressure, const ADReal & enthalpy) const;
     768             : 
     769             :   /**
     770             :    * Backwards equation T(p, h) in Region 2c
     771             :    * Eq. (24) from Revised Release on the IAPWS Industrial
     772             :    * Formulation 1997 for the Thermodynamic Properties of Water
     773             :    * and Steam
     774             :    *
     775             :    * @param pressure water pressure (Pa)
     776             :    * @param enthalpy water enthalpy (J/kg)
     777             :    * @return temperature water temperature (K)
     778             :    */
     779             :   ADReal temperature_from_ph2c(const ADReal & pressure, const ADReal & enthalpy) const;
     780             : 
     781             :   /**
     782             :    * Backwards equation T(p, h) in Region 3a
     783             :    * Eq. (2) from Revised Supplementary Release on Backward Equations for
     784             :    * the Functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS
     785             :    * Industrial Formulation 1997 for the Thermodynamic Properties of Water and
     786             :    * Steam
     787             :    *
     788             :    * @param pressure water pressure (Pa)
     789             :    * @param enthalpy water enthalpy (J/kg)
     790             :    * @return temperature water temperature (K)
     791             :    */
     792             :   ADReal temperature_from_ph3a(const ADReal & pressure, const ADReal & enthalpy) const;
     793             : 
     794             :   /**
     795             :    * Backwards equation T(p, h) in Region 3b
     796             :    * Eq. (3) from Revised Supplementary Release on Backward Equations for
     797             :    * the Functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS
     798             :    * Industrial Formulation 1997 for the Thermodynamic Properties of Water and
     799             :    * Steam
     800             :    *
     801             :    * @param pressure water pressure (Pa)
     802             :    * @param enthalpy water enthalpy (J/kg)
     803             :    * @return temperature water temperature (K)
     804             :    */
     805             :   ADReal temperature_from_ph3b(const ADReal & pressure, const ADReal & enthalpy) const;
     806             : 
     807             :   /**
     808             :    * AD version of saturation temperature as a function of pressure (used internally)
     809             :    *
     810             :    * Eq. (31) from Revised Release on the IAPWS Industrial
     811             :    * Formulation 1997 for the Thermodynamic Properties of Water
     812             :    * and Steam
     813             :    *
     814             :    * Valid for 611.213 Pa <= p <= 22.064 MPa
     815             :    *
     816             :    * @param pressure water pressure (Pa)
     817             :    * @return saturation temperature (K)
     818             :    */
     819             :   ADReal vaporTemperature_ad(const ADReal & pressure) const;
     820             : 
     821             :   /// Water molar mass (kg/mol)
     822             :   const Real _Mh2o;
     823             :   /// Specific gas constant for H2O (universal gas constant / molar mass of water - J/kg/K)
     824             :   const Real _Rw;
     825             :   /// Critical pressure (Pa)
     826             :   const Real _p_critical;
     827             :   /// Critical temperature (K)
     828             :   const Real _T_critical;
     829             :   /// Critical density (kg/m^3)
     830             :   const Real _rho_critical;
     831             :   /// Triple point pressure (Pa)
     832             :   const Real _p_triple;
     833             :   /// Triple point temperature (K)
     834             :   const Real _T_triple;
     835             : 
     836             :   /**
     837             :    * Reference constants used in to calculate thermophysical properties of water.
     838             :    * Taken from Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic
     839             :    * Properties of Water and Steam and from Revised Supplementary Release on Backward Equations for
     840             :    * Specific Volume as a Function of Pressure and Temperature v(p,T) for Region 3 of the IAPWS
     841             :    * Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam
     842             :    */
     843             : 
     844             :   /// Constants for region 1
     845             :   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             : 
     875             :   /// Constants for region 2
     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             : 
     960             :   /// Constants for the boundary between regions 2 and 3
     961             :   const std::array<Real, 5> _n23{{0.34805185628969e3,
     962             :                                   -0.11671859879975e1,
     963             :                                   0.10192970039326e-2,
     964             :                                   0.57254459862746e3,
     965             :                                   0.13918839778870e2}};
     966             : 
     967             :   /// Constants for region 3
     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             : 
    1018             :   /// Constants for all 26 subregions in region 3
    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             : 
    1404             :   /// Constants for region 4 (the saturation curve up to the critical point)
    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             : 
    1416             :   /// Constants for region 5
    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             : 
    1434             :   /// Constnats for the tempXY() method
    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             : 
    1484             :   /// Constants from Release on the IAPWS Formulation 2008 for the Viscosity of
    1485             :   /// Ordinary Water Substance
    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             : 
    1495             :   /// Constants for thermal conductivity
    1496             :   std::array<Real, 4> _k_a{{0.0102811, 0.0299621, 0.0156146, -0.00422464}};
    1497             : 
    1498             :   /// Temperature scale for each region
    1499             :   const std::array<Real, 5> _T_star{{1386.0, 540.0, _T_critical, 1.0, 1000.0}};
    1500             :   /// Pressure scale for each region
    1501             :   const std::array<Real, 5> _p_star{{16.53e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6}};
    1502             : 
    1503             : private:
    1504             :   /**
    1505             :    * Computes the temperature (first member of the pair) and the derivative of density (second
    1506             :    * member of the pair) with respect to temperature as a function of pressure and density
    1507             :    */
    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      653074 : Water97FluidProperties::dgamma1_dpi(const T & pi, const T & tau) const
    1517             : {
    1518         336 :   T sum = 0.0;
    1519    22857590 :   for (std::size_t i = 0; i < _n1.size(); ++i)
    1520    22227364 :     sum += -_n1[i] * _I1[i] * MathUtils::pow(7.1 - pi, _I1[i] - 1) *
    1521    22190508 :            MathUtils::pow(tau - 1.222, _J1[i]);
    1522             : 
    1523      653074 :   return sum;
    1524             : }
    1525             : 
    1526             : template <typename T>
    1527             : T
    1528        9201 : Water97FluidProperties::dgamma2_dpi(const T & pi, const T & tau) const
    1529             : {
    1530             :   // Ideal gas part of the Gibbs free energy
    1531        9201 :   T dg0 = 1.0 / pi;
    1532             : 
    1533             :   // Residual part of the Gibbs free energy
    1534           0 :   T dgr = 0.0;
    1535      404844 :   for (std::size_t i = 0; i < _n2.size(); ++i)
    1536      395643 :     dgr += _n2[i] * _I2[i] * MathUtils::pow(pi, _I2[i] - 1) * MathUtils::pow(tau - 0.5, _J2[i]);
    1537             : 
    1538        9201 :   return dg0 + dgr;
    1539             : }
    1540             : 
    1541             : template <typename T>
    1542             : T
    1543         102 : Water97FluidProperties::densityRegion3(const T & pressure, const T & temperature) const
    1544             : {
    1545             :   // Region 3 is subdivided into 26 subregions, each with a given backwards equation
    1546             :   // to directly calculate density from pressure and temperature without the need for
    1547             :   // expensive iterations. Find the subregion id that the point is in:
    1548             :   unsigned int sid =
    1549         102 :       subregion3(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    1550             : 
    1551           0 :   T vstar, pi, theta;
    1552             :   Real a, b, c, d, e;
    1553             :   unsigned int N;
    1554             : 
    1555         102 :   vstar = _par3[sid][0];
    1556         105 :   pi = pressure / _par3[sid][1] / 1.0e6;
    1557         102 :   theta = temperature / _par3[sid][2];
    1558         102 :   a = _par3[sid][3];
    1559         102 :   b = _par3[sid][4];
    1560         102 :   c = _par3[sid][5];
    1561         102 :   d = _par3[sid][6];
    1562         102 :   e = _par3[sid][7];
    1563         102 :   N = _par3N[sid];
    1564             : 
    1565           0 :   T sum = 0.0;
    1566           0 :   T volume = 0.0;
    1567             : 
    1568             :   // Note that subregion 13 is the only different formulation
    1569         102 :   if (sid == 13)
    1570             :   {
    1571          80 :     for (std::size_t i = 0; i < N; ++i)
    1572          78 :       sum += _n3s[sid][i] * MathUtils::pow(pi - a, _I3s[sid][i]) *
    1573          78 :              MathUtils::pow(theta - b, _J3s[sid][i]);
    1574             : 
    1575           2 :     volume = vstar * std::exp(sum);
    1576             :   }
    1577             :   else
    1578         103 :     volume = vstar * subregionVolume(pi, theta, a, b, c, d, e, sid);
    1579             : 
    1580             :   // Density is the inverse of volume
    1581         102 :   return 1.0 / volume;
    1582             : }
    1583             : 
    1584             : template <typename T>
    1585             : T
    1586         242 : Water97FluidProperties::gamma1(const T & pi, const T & tau) const
    1587             : {
    1588           5 :   T sum = 0.0;
    1589        8470 :   for (std::size_t i = 0; i < _n1.size(); ++i)
    1590        8568 :     sum += _n1[i] * MathUtils::pow(7.1 - pi, _I1[i]) * MathUtils::pow(tau - 1.222, _J1[i]);
    1591             : 
    1592         242 :   return sum;
    1593             : }
    1594             : 
    1595             : template <typename T>
    1596             : T
    1597         457 : Water97FluidProperties::d2gamma1_dpi2(const T & pi, const T & tau) const
    1598             : {
    1599          11 :   T sum = 0.0;
    1600       15995 :   for (std::size_t i = 0; i < _n1.size(); ++i)
    1601       16286 :     sum += _n1[i] * _I1[i] * (_I1[i] - 1) * MathUtils::pow(7.1 - pi, _I1[i] - 2) *
    1602       15164 :            MathUtils::pow(tau - 1.222, _J1[i]);
    1603             : 
    1604         457 :   return sum;
    1605             : }
    1606             : 
    1607             : template <typename T>
    1608             : T
    1609      435224 : Water97FluidProperties::dgamma1_dtau(const T & pi, const T & tau) const
    1610             : {
    1611          61 :   T g = 0.0;
    1612    15232840 :   for (std::size_t i = 0; i < _n1.size(); ++i)
    1613    14801764 :     g += _n1[i] * _J1[i] * MathUtils::pow(7.1 - pi, _I1[i]) *
    1614    14794488 :          MathUtils::pow(tau - 1.222, _J1[i] - 1);
    1615             : 
    1616      435224 :   return g;
    1617             : }
    1618             : 
    1619             : template <typename T>
    1620             : T
    1621        1119 : Water97FluidProperties::d2gamma1_dtau2(const T & pi, const T & tau) const
    1622             : {
    1623          16 :   T dg = 0.0;
    1624       39165 :   for (std::size_t i = 0; i < _n1.size(); ++i)
    1625       39134 :     dg += _n1[i] * _J1[i] * (_J1[i] - 1) * MathUtils::pow(7.1 - pi, _I1[i]) *
    1626       37502 :           MathUtils::pow(tau - 1.222, _J1[i] - 2);
    1627             : 
    1628        1119 :   return dg;
    1629             : }
    1630             : 
    1631             : template <typename T>
    1632             : T
    1633         457 : Water97FluidProperties::d2gamma1_dpitau(const T & pi, const T & tau) const
    1634             : {
    1635          11 :   T dg = 0.0;
    1636       15995 :   for (std::size_t i = 0; i < _n1.size(); ++i)
    1637       16286 :     dg += -_n1[i] * _I1[i] * _J1[i] * MathUtils::pow(7.1 - pi, _I1[i] - 1) *
    1638       15164 :           MathUtils::pow(tau - 1.222, _J1[i] - 1);
    1639             : 
    1640         457 :   return dg;
    1641             : }
    1642             : 
    1643             : template <typename T>
    1644             : T
    1645          12 : Water97FluidProperties::gamma2(const T & pi, const T & tau) const
    1646             : {
    1647             :   // Ideal gas part of the Gibbs free energy
    1648           0 :   T sum0 = 0.0;
    1649         120 :   for (std::size_t i = 0; i < _n02.size(); ++i)
    1650         108 :     sum0 += _n02[i] * MathUtils::pow(tau, _J02[i]);
    1651             : 
    1652          13 :   T g0 = std::log(pi) + sum0;
    1653             : 
    1654             :   // Residual part of the Gibbs free energy
    1655           0 :   T gr = 0.0;
    1656         528 :   for (std::size_t i = 0; i < _n2.size(); ++i)
    1657         516 :     gr += _n2[i] * MathUtils::pow(pi, _I2[i]) * MathUtils::pow(tau - 0.5, _J2[i]);
    1658             : 
    1659          12 :   return g0 + gr;
    1660             : }
    1661             : 
    1662             : template <typename T>
    1663             : T
    1664           5 : Water97FluidProperties::d2gamma2_dpi2(const T & pi, const T & tau) const
    1665             : {
    1666             :   // Ideal gas part of the Gibbs free energy
    1667           5 :   T dg0 = -1.0 / pi / pi;
    1668             : 
    1669             :   // Residual part of the Gibbs free energy
    1670           0 :   T dgr = 0.0;
    1671         220 :   for (std::size_t i = 0; i < _n2.size(); ++i)
    1672         215 :     dgr += _n2[i] * _I2[i] * (_I2[i] - 1) * MathUtils::pow(pi, _I2[i] - 2) *
    1673         215 :            MathUtils::pow(tau - 0.5, _J2[i]);
    1674             : 
    1675           5 :   return dg0 + dgr;
    1676             : }
    1677             : 
    1678             : template <typename T>
    1679             : T
    1680        6190 : Water97FluidProperties::dgamma2_dtau(const T & pi, const T & tau) const
    1681             : {
    1682             :   // Ideal gas part of the Gibbs free energy
    1683           0 :   T dg0 = 0.0;
    1684       61900 :   for (std::size_t i = 0; i < _n02.size(); ++i)
    1685       55710 :     dg0 += _n02[i] * _J02[i] * MathUtils::pow(tau, _J02[i] - 1);
    1686             : 
    1687             :   // Residual part of the Gibbs free energy
    1688           0 :   T dgr = 0.0;
    1689      272360 :   for (std::size_t i = 0; i < _n2.size(); ++i)
    1690      266170 :     dgr += _n2[i] * _J2[i] * MathUtils::pow(pi, _I2[i]) * MathUtils::pow(tau - 0.5, _J2[i] - 1);
    1691             : 
    1692        6190 :   return dg0 + dgr;
    1693             : }
    1694             : 
    1695             : template <typename T>
    1696             : T
    1697           8 : Water97FluidProperties::d2gamma2_dtau2(const T & pi, const T & tau) const
    1698             : {
    1699             :   // Ideal gas part of the Gibbs free energy
    1700           0 :   T dg0 = 0.0;
    1701          80 :   for (std::size_t i = 0; i < _n02.size(); ++i)
    1702          72 :     dg0 += _n02[i] * _J02[i] * (_J02[i] - 1) * MathUtils::pow(tau, _J02[i] - 2);
    1703             : 
    1704             :   // Residual part of the Gibbs free energy
    1705           0 :   T dgr = 0.0;
    1706         352 :   for (std::size_t i = 0; i < _n2.size(); ++i)
    1707         344 :     dgr += _n2[i] * _J2[i] * (_J2[i] - 1) * MathUtils::pow(pi, _I2[i]) *
    1708         344 :            MathUtils::pow(tau - 0.5, _J2[i] - 2);
    1709             : 
    1710           8 :   return dg0 + dgr;
    1711             : }
    1712             : 
    1713             : template <typename T>
    1714             : T
    1715           5 : Water97FluidProperties::d2gamma2_dpitau(const T & pi, const T & tau) const
    1716             : {
    1717             :   // Ideal gas part of the Gibbs free energy
    1718           0 :   T dg0 = 0.0;
    1719             : 
    1720             :   // Residual part of the Gibbs free energy
    1721           0 :   T dgr = 0.0;
    1722         220 :   for (std::size_t i = 0; i < _n2.size(); ++i)
    1723         215 :     dgr += _n2[i] * _I2[i] * _J2[i] * MathUtils::pow(pi, _I2[i] - 1) *
    1724         215 :            MathUtils::pow(tau - 0.5, _J2[i] - 1);
    1725             : 
    1726           5 :   return dg0 + dgr;
    1727             : }
    1728             : 
    1729             : template <typename T>
    1730             : T
    1731          12 : Water97FluidProperties::phi3(const T & delta, const T & tau) const
    1732             : {
    1733           0 :   T sum = 0.0;
    1734         480 :   for (std::size_t i = 1; i < _n3.size(); ++i)
    1735         468 :     sum += _n3[i] * MathUtils::pow(delta, _I3[i]) * MathUtils::pow(tau, _J3[i]);
    1736             : 
    1737          12 :   return _n3[0] * std::log(delta) + sum;
    1738             : }
    1739             : 
    1740             : template <typename T>
    1741             : T
    1742          19 : Water97FluidProperties::dphi3_ddelta(const T & delta, const T & tau) const
    1743             : {
    1744           0 :   T sum = 0.0;
    1745         760 :   for (std::size_t i = 1; i < _n3.size(); ++i)
    1746         741 :     sum += _n3[i] * _I3[i] * MathUtils::pow(delta, _I3[i] - 1) * MathUtils::pow(tau, _J3[i]);
    1747             : 
    1748          19 :   return _n3[0] / delta + sum;
    1749             : }
    1750             : 
    1751             : template <typename T>
    1752             : T
    1753           7 : Water97FluidProperties::d2phi3_ddelta2(const T & delta, const T & tau) const
    1754             : {
    1755           0 :   T sum = 0.0;
    1756         280 :   for (std::size_t i = 1; i < _n3.size(); ++i)
    1757         273 :     sum += _n3[i] * _I3[i] * (_I3[i] - 1) * MathUtils::pow(delta, _I3[i] - 2) *
    1758         273 :            MathUtils::pow(tau, _J3[i]);
    1759             : 
    1760           7 :   return -_n3[0] / delta / delta + sum;
    1761             : }
    1762             : 
    1763             : template <typename T>
    1764             : T
    1765          33 : Water97FluidProperties::dphi3_dtau(const T & delta, const T & tau) const
    1766             : {
    1767           0 :   T sum = 0.0;
    1768        1320 :   for (std::size_t i = 1; i < _n3.size(); ++i)
    1769        1287 :     sum += _n3[i] * _J3[i] * MathUtils::pow(delta, _I3[i]) * MathUtils::pow(tau, _J3[i] - 1);
    1770             : 
    1771          33 :   return sum;
    1772             : }
    1773             : 
    1774             : template <typename T>
    1775             : T
    1776           8 : Water97FluidProperties::d2phi3_dtau2(const T & delta, const T & tau) const
    1777             : {
    1778           0 :   T sum = 0.0;
    1779         320 :   for (std::size_t i = 1; i < _n3.size(); ++i)
    1780         312 :     sum += _n3[i] * _J3[i] * (_J3[i] - 1) * MathUtils::pow(delta, _I3[i]) *
    1781         312 :            MathUtils::pow(tau, _J3[i] - 2);
    1782             : 
    1783           8 :   return sum;
    1784             : }
    1785             : 
    1786             : template <typename T>
    1787             : T
    1788           7 : Water97FluidProperties::d2phi3_ddeltatau(const T & delta, const T & tau) const
    1789             : {
    1790           0 :   T sum = 0.0;
    1791         280 :   for (std::size_t i = 1; i < _n3.size(); ++i)
    1792         273 :     sum += _n3[i] * _I3[i] * _J3[i] * MathUtils::pow(delta, _I3[i] - 1) *
    1793         273 :            MathUtils::pow(tau, _J3[i] - 1);
    1794             : 
    1795           7 :   return sum;
    1796             : }
    1797             : 
    1798             : template <typename T>
    1799             : T
    1800           9 : Water97FluidProperties::gamma5(const T & pi, const T & tau) const
    1801             : {
    1802             :   // Ideal gas part of the Gibbs free energy
    1803           0 :   T sum0 = 0.0;
    1804          63 :   for (std::size_t i = 0; i < _n05.size(); ++i)
    1805          54 :     sum0 += _n05[i] * MathUtils::pow(tau, _J05[i]);
    1806             : 
    1807          10 :   T g0 = std::log(pi) + sum0;
    1808             : 
    1809             :   // Residual part of the Gibbs free energy
    1810           0 :   T gr = 0.0;
    1811          63 :   for (std::size_t i = 0; i < _n5.size(); ++i)
    1812          54 :     gr += _n5[i] * MathUtils::pow(pi, _I5[i]) * MathUtils::pow(tau, _J5[i]);
    1813             : 
    1814           9 :   return g0 + gr;
    1815             : }
    1816             : 
    1817             : template <typename T>
    1818             : T
    1819          25 : Water97FluidProperties::dgamma5_dpi(const T & pi, const T & tau) const
    1820             : {
    1821             :   // Ideal gas part of the Gibbs free energy
    1822          25 :   T dg0 = 1.0 / pi;
    1823             : 
    1824             :   // Residual part of the Gibbs free energy
    1825           0 :   T dgr = 0.0;
    1826         175 :   for (std::size_t i = 0; i < _n5.size(); ++i)
    1827         150 :     dgr += _n5[i] * _I5[i] * MathUtils::pow(pi, _I5[i] - 1) * MathUtils::pow(tau, _J5[i]);
    1828             : 
    1829          25 :   return dg0 + dgr;
    1830             : }
    1831             : 
    1832             : template <typename T>
    1833             : T
    1834           5 : Water97FluidProperties::d2gamma5_dpi2(const T & pi, const T & tau) const
    1835             : {
    1836             :   // Ideal gas part of the Gibbs free energy
    1837           5 :   T dg0 = -1.0 / pi / pi;
    1838             : 
    1839             :   // Residual part of the Gibbs free energy
    1840           0 :   T dgr = 0.0;
    1841          35 :   for (std::size_t i = 0; i < _n5.size(); ++i)
    1842          30 :     dgr += _n5[i] * _I5[i] * (_I5[i] - 1) * MathUtils::pow(pi, _I5[i] - 2) *
    1843          30 :            MathUtils::pow(tau, _J5[i]);
    1844             : 
    1845           5 :   return dg0 + dgr;
    1846             : }
    1847             : 
    1848             : template <typename T>
    1849             : T
    1850          27 : Water97FluidProperties::dgamma5_dtau(const T & pi, const T & tau) const
    1851             : {
    1852             :   // Ideal gas part of the Gibbs free energy
    1853           0 :   T dg0 = 0.0;
    1854         189 :   for (std::size_t i = 0; i < _n05.size(); ++i)
    1855         162 :     dg0 += _n05[i] * _J05[i] * MathUtils::pow(tau, _J05[i] - 1);
    1856             : 
    1857             :   // Residual part of the Gibbs free energy
    1858           0 :   T dgr = 0.0;
    1859         189 :   for (std::size_t i = 0; i < _n5.size(); ++i)
    1860         162 :     dgr += _n5[i] * _J5[i] * MathUtils::pow(pi, _I5[i]) * MathUtils::pow(tau, _J5[i] - 1);
    1861             : 
    1862          27 :   return dg0 + dgr;
    1863             : }
    1864             : 
    1865             : template <typename T>
    1866             : T
    1867           8 : Water97FluidProperties::d2gamma5_dtau2(const T & pi, const T & tau) const
    1868             : {
    1869             :   // Ideal gas part of the Gibbs free energy
    1870           0 :   T dg0 = 0.0;
    1871          56 :   for (std::size_t i = 0; i < _n05.size(); ++i)
    1872          48 :     dg0 += _n05[i] * _J05[i] * (_J05[i] - 1) * MathUtils::pow(tau, _J05[i] - 2);
    1873             : 
    1874             :   // Residual part of the Gibbs free energy
    1875           0 :   T dgr = 0.0;
    1876          56 :   for (std::size_t i = 0; i < _n5.size(); ++i)
    1877          48 :     dgr += _n5[i] * _J5[i] * (_J5[i] - 1) * MathUtils::pow(pi, _I5[i]) *
    1878          48 :            MathUtils::pow(tau, _J5[i] - 2);
    1879             : 
    1880           8 :   return dg0 + dgr;
    1881             : }
    1882             : 
    1883             : template <typename T>
    1884             : T
    1885           5 : Water97FluidProperties::d2gamma5_dpitau(const T & pi, const T & tau) const
    1886             : {
    1887             :   // Ideal gas part of the Gibbs free energy
    1888           0 :   T dg0 = 0.0;
    1889             : 
    1890             :   // Residual part of the Gibbs free energy
    1891           0 :   T dgr = 0.0;
    1892          35 :   for (std::size_t i = 0; i < _n5.size(); ++i)
    1893          30 :     dgr +=
    1894          30 :         _n5[i] * _I5[i] * _J5[i] * MathUtils::pow(pi, _I5[i] - 1) * MathUtils::pow(tau, _J5[i] - 1);
    1895             : 
    1896           5 :   return dg0 + dgr;
    1897             : }
    1898             : 
    1899             : template <typename T>
    1900             : T
    1901         682 : Water97FluidProperties::tempXY(const T & pressure, const subregionEnum xy) const
    1902             : {
    1903         682 :   T pi = pressure / 1.0e6;
    1904             : 
    1905             :   // Choose the constants based on the string xy
    1906             :   unsigned int row;
    1907             : 
    1908             :   switch (xy)
    1909             :   {
    1910             :     case AB:
    1911             :       row = 0;
    1912             :       break;
    1913             :     case CD:
    1914             :       row = 1;
    1915             :       break;
    1916             :     case GH:
    1917             :       row = 2;
    1918             :       break;
    1919             :     case IJ:
    1920             :       row = 3;
    1921             :       break;
    1922             :     case JK:
    1923             :       row = 4;
    1924             :       break;
    1925             :     case MN:
    1926             :       row = 5;
    1927             :       break;
    1928             :     case OP:
    1929             :       row = 6;
    1930             :       break;
    1931             :     case QU:
    1932             :       row = 7;
    1933             :       break;
    1934             :     case RX:
    1935             :       row = 8;
    1936             :       break;
    1937             :     case UV:
    1938             :       row = 9;
    1939             :       break;
    1940             :     case WX:
    1941             :       row = 10;
    1942             :       break;
    1943             :     default:
    1944             :       row = 0;
    1945             :   }
    1946             : 
    1947             :   T sum = 0.0;
    1948             : 
    1949         682 :   if (xy == AB || xy == OP || xy == WX)
    1950         384 :     for (std::size_t i = 0; i < _tempXY_n[row].size(); ++i)
    1951         320 :       sum += _tempXY_n[row][i] * MathUtils::pow(std::log(pi), _tempXY_I[row][i]);
    1952         618 :   else if (xy == EF)
    1953          76 :     sum += 3.727888004 * (pi - _p_critical / 1.0e6) + _T_critical;
    1954             :   else
    1955        2826 :     for (std::size_t i = 0; i < _tempXY_n[row].size(); ++i)
    1956        2284 :       sum += _tempXY_n[row][i] * MathUtils::pow(pi, _tempXY_I[row][i]);
    1957             : 
    1958         682 :   return sum;
    1959             : }
    1960             : 
    1961             : template <typename T>
    1962             : void
    1963           1 : Water97FluidProperties::vaporPressureTemplate(const T & temperature, T & psat, T & dpsat_dT) const
    1964             : {
    1965             :   // Check whether the input temperature is within the region of validity of this equation.
    1966             :   // Valid for 273.15 K <= t <= 647.096 K
    1967           1 :   if (temperature < 273.15 || temperature > _T_critical)
    1968           0 :     mooseException(name(),
    1969             :                    ": vaporPressure(): Temperature is outside range 273.15 K <= T <= 647.096 K");
    1970             : 
    1971             :   T theta, dtheta_dT, theta2, a, b, c, da_dtheta, db_dtheta, dc_dtheta;
    1972           1 :   theta = temperature + _n4[8] / (temperature - _n4[9]);
    1973           1 :   dtheta_dT = 1.0 - _n4[8] / (temperature - _n4[9]) / (temperature - _n4[9]);
    1974           1 :   theta2 = theta * theta;
    1975             : 
    1976           1 :   a = theta2 + _n4[0] * theta + _n4[1];
    1977           1 :   b = _n4[2] * theta2 + _n4[3] * theta + _n4[4];
    1978           1 :   c = _n4[5] * theta2 + _n4[6] * theta + _n4[7];
    1979             : 
    1980           1 :   da_dtheta = 2.0 * theta + _n4[0];
    1981           1 :   db_dtheta = 2.0 * _n4[2] * theta + _n4[3];
    1982           1 :   dc_dtheta = 2.0 * _n4[5] * theta + _n4[6];
    1983             : 
    1984           1 :   T denominator = -b + std::sqrt(b * b - 4.0 * a * c);
    1985             : 
    1986           1 :   psat = Utility::pow<4>(2.0 * c / denominator) * 1.0e6;
    1987             : 
    1988             :   // The derivative wrt temperature is given by the chain rule
    1989           1 :   T dpsat = 4.0 * Utility::pow<3>(2.0 * c / denominator);
    1990           1 :   dpsat *= (2.0 * dc_dtheta / denominator -
    1991           1 :             2.0 * c / denominator / denominator *
    1992           1 :                 (-db_dtheta + std::pow(b * b - 4.0 * a * c, -0.5) *
    1993           1 :                                   (b * db_dtheta - 2.0 * da_dtheta * c - 2.0 * a * dc_dtheta)));
    1994           1 :   dpsat_dT = dpsat * dtheta_dT * 1.0e6;
    1995           1 : }
    1996             : 
    1997             : inline void
    1998           1 : Water97FluidProperties::vaporPressure(const Real temperature, Real & psat, Real & dpsat_dT) const
    1999             : {
    2000           1 :   vaporPressureTemplate(temperature, psat, dpsat_dT);
    2001           1 : }
    2002             : 
    2003             : template <typename T>
    2004             : void
    2005         256 : Water97FluidProperties::rho_from_p_T_template(
    2006             :     const T & pressure, const T & temperature, T & rho, T & drho_dp, T & drho_dT) const
    2007             : {
    2008         316 :   auto functor = [this](const auto & pressure, const auto & temperature)
    2009         316 :   { return this->rho_from_p_T_template(pressure, temperature); };
    2010             : 
    2011         276 :   xyDerivatives(pressure, temperature, rho, drho_dp, drho_dT, functor);
    2012         256 : }
    2013             : 
    2014             : template <typename T>
    2015             : T
    2016          16 : Water97FluidProperties::v_from_p_T_template(const T & pressure, const T & temperature) const
    2017             : {
    2018          16 :   return 1 / rho_from_p_T_template(pressure, temperature);
    2019             : }
    2020             : 
    2021             : template <typename T>
    2022             : void
    2023          12 : Water97FluidProperties::v_from_p_T_template(
    2024             :     const T & pressure, const T & temperature, T & v, T & dv_dp, T & dv_dT) const
    2025             : {
    2026          16 :   auto functor = [this](const auto & pressure, const auto & temperature)
    2027          16 :   { return this->v_from_p_T_template(pressure, temperature); };
    2028             : 
    2029          16 :   xyDerivatives(pressure, temperature, v, dv_dp, dv_dT, functor);
    2030          12 : }
    2031             : 
    2032             : template <typename T>
    2033             : T
    2034      220316 : Water97FluidProperties::e_from_p_T_template(const T & pressure, const T & temperature) const
    2035             : {
    2036             :   T internal_energy, pi, tau;
    2037             : 
    2038             :   // Determine which region the point is in
    2039             :   unsigned int region =
    2040      220316 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2041      220316 :   switch (region)
    2042             :   {
    2043      217231 :     case 1:
    2044      217231 :       pi = pressure / _p_star[0];
    2045      217231 :       tau = _T_star[0] / temperature;
    2046      217311 :       internal_energy =
    2047      217231 :           _Rw * temperature * (tau * dgamma1_dtau(pi, tau) - pi * dgamma1_dpi(pi, tau));
    2048      217231 :       break;
    2049             : 
    2050        3067 :     case 2:
    2051        3067 :       pi = pressure / _p_star[1];
    2052        3067 :       tau = _T_star[1] / temperature;
    2053        3067 :       internal_energy =
    2054        3067 :           _Rw * temperature * (tau * dgamma2_dtau(pi, tau) - pi * dgamma2_dpi(pi, tau));
    2055        3067 :       break;
    2056             : 
    2057           9 :     case 3:
    2058             :     {
    2059             :       // Calculate density first, then use that in Helmholtz free energy
    2060           9 :       T density3 = densityRegion3(pressure, temperature);
    2061           9 :       T delta = density3 / _rho_critical;
    2062           9 :       tau = _T_star[2] / temperature;
    2063           9 :       internal_energy = _Rw * temperature * tau * dphi3_dtau(delta, tau);
    2064             :       break;
    2065             :     }
    2066             : 
    2067           9 :     case 5:
    2068           9 :       pi = pressure / _p_star[4];
    2069           9 :       tau = _T_star[4] / temperature;
    2070           9 :       internal_energy =
    2071           9 :           _Rw * temperature * (tau * dgamma5_dtau(pi, tau) - pi * dgamma5_dpi(pi, tau));
    2072           9 :       break;
    2073             : 
    2074           0 :     default:
    2075           0 :       mooseError("inRegion() has given an incorrect region");
    2076             :   }
    2077             :   // Output in J/kg
    2078      220316 :   return internal_energy;
    2079             : }
    2080             : 
    2081             : template <typename T>
    2082             : void
    2083           8 : Water97FluidProperties::e_from_p_T_template(
    2084             :     const T & pressure, const T & temperature, T & e, T & de_dp, T & de_dT) const
    2085             : {
    2086             :   T pi, tau, dinternal_energy_dp, dinternal_energy_dT;
    2087             : 
    2088             :   // Determine which region the point is in
    2089             :   unsigned int region =
    2090           8 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2091           8 :   switch (region)
    2092             :   {
    2093           5 :     case 1:
    2094             :     {
    2095           5 :       pi = pressure / _p_star[0];
    2096           5 :       tau = _T_star[0] / temperature;
    2097           5 :       T dgdp = dgamma1_dpi(pi, tau);
    2098           5 :       T d2gdpt = d2gamma1_dpitau(pi, tau);
    2099           6 :       dinternal_energy_dp =
    2100           5 :           _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma1_dpi2(pi, tau)) / _p_star[0];
    2101           6 :       dinternal_energy_dT =
    2102           4 :           _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma1_dtau2(pi, tau) - pi * dgdp);
    2103           4 :       break;
    2104             :     }
    2105             : 
    2106           1 :     case 2:
    2107             :     {
    2108           1 :       pi = pressure / _p_star[1];
    2109           1 :       tau = _T_star[1] / temperature;
    2110           1 :       T dgdp = dgamma2_dpi(pi, tau);
    2111           1 :       T d2gdpt = d2gamma2_dpitau(pi, tau);
    2112           1 :       dinternal_energy_dp =
    2113           1 :           _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma2_dpi2(pi, tau)) / _p_star[1];
    2114           1 :       dinternal_energy_dT =
    2115           1 :           _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma2_dtau2(pi, tau) - pi * dgdp);
    2116           1 :       break;
    2117             :     }
    2118             : 
    2119           1 :     case 3:
    2120             :     {
    2121             :       // Calculate density first, then use that in Helmholtz free energy
    2122           1 :       T density3 = densityRegion3(pressure, temperature);
    2123           1 :       T delta = density3 / _rho_critical;
    2124           1 :       tau = _T_star[2] / temperature;
    2125           1 :       T dpdd = dphi3_ddelta(delta, tau);
    2126           1 :       T d2pddt = d2phi3_ddeltatau(delta, tau);
    2127           1 :       T d2pdd2 = d2phi3_ddelta2(delta, tau);
    2128           1 :       dinternal_energy_dp =
    2129           1 :           _T_star[2] * d2pddt / _rho_critical /
    2130           1 :           (2.0 * temperature * delta * dpdd + temperature * delta * delta * d2pdd2);
    2131           1 :       dinternal_energy_dT =
    2132           1 :           -_Rw * (delta * tau * d2pddt * (dpdd - tau * d2pddt) / (2.0 * dpdd + delta * d2pdd2) +
    2133           1 :                   tau * tau * d2phi3_dtau2(delta, tau));
    2134             :       break;
    2135             :     }
    2136             : 
    2137           1 :     case 5:
    2138             :     {
    2139           1 :       pi = pressure / _p_star[4];
    2140           1 :       tau = _T_star[4] / temperature;
    2141           1 :       T dgdp = dgamma5_dpi(pi, tau);
    2142           1 :       T d2gdpt = d2gamma5_dpitau(pi, tau);
    2143           1 :       dinternal_energy_dp =
    2144           1 :           _Rw * temperature * (tau * d2gdpt - dgdp - pi * d2gamma5_dpi2(pi, tau)) / _p_star[4];
    2145           1 :       dinternal_energy_dT =
    2146           1 :           _Rw * (pi * tau * d2gdpt - tau * tau * d2gamma5_dtau2(pi, tau) - pi * dgdp);
    2147           1 :       break;
    2148             :     }
    2149             : 
    2150           0 :     default:
    2151           0 :       mooseError("inRegion has given an incorrect region");
    2152             :   }
    2153             : 
    2154           8 :   e = this->e_from_p_T(pressure, temperature);
    2155           8 :   de_dp = dinternal_energy_dp;
    2156           8 :   de_dT = dinternal_energy_dT;
    2157           8 : }
    2158             : 
    2159             : template <typename T>
    2160             : T
    2161         100 : Water97FluidProperties::subregionVolume(
    2162             :     const T & pi, const T & theta, Real a, Real b, Real c, Real d, Real e, unsigned int sid) const
    2163             : {
    2164           0 :   T sum = 0.0;
    2165             : 
    2166        3394 :   for (std::size_t i = 0; i < _n3s[sid].size(); ++i)
    2167        3294 :     sum += _n3s[sid][i] * MathUtils::pow(std::pow(pi - a, c), _I3s[sid][i]) *
    2168        3189 :            MathUtils::pow(std::pow(theta - b, d), _J3s[sid][i]);
    2169             : 
    2170         100 :   return std::pow(sum, e);
    2171             : }
    2172             : 
    2173             : template <typename T>
    2174             : std::pair<T, T>
    2175          74 : Water97FluidProperties::T_drhodT_from_p_rho(const T & p, const T & rho) const
    2176             : {
    2177             :   // For NewtonSolve
    2178             :   // y = rho
    2179             :   // x = p
    2180             :   // z = T
    2181          74 :   auto lambda =
    2182         256 :       [this](const T & p, const T & temperature, T & rho, T & drho_dp, T & drho_dtemperature)
    2183         296 :   { rho_from_p_T_template(p, temperature, rho, drho_dp, drho_dtemperature); };
    2184          10 :   return FluidPropertiesUtils::NewtonSolve(
    2185         148 :       p, rho, _T_initial_guess, _tolerance, lambda, name() + "::T_drhodT_from_p_rho");
    2186             : }
    2187             : 
    2188             : template <typename T>
    2189             : std::pair<T, T>
    2190          18 : Water97FluidProperties::p_T_from_v_e(const T & v, const T & e) const
    2191             : {
    2192          18 :   const auto rho = 1 / v;
    2193          18 :   const auto p = p_from_v_e(v, e);
    2194          18 :   return {p, T_drhodT_from_p_rho(p, rho).first};
    2195             : }
    2196             : 
    2197             : template <typename T>
    2198             : std::pair<T, T>
    2199           4 : Water97FluidProperties::rho_T_from_v_e(const T & v, const T & e) const
    2200             : {
    2201           4 :   const auto rho = 1 / v;
    2202           4 :   const auto p = p_from_v_e(v, e);
    2203           4 :   return {rho, T_drhodT_from_p_rho(p, rho).first};
    2204             : }
    2205             : 
    2206             : template <typename T>
    2207             : T
    2208         236 : Water97FluidProperties::c_from_p_T_template(const T & pressure, const T & temperature) const
    2209             : {
    2210             :   T speed2, pi, tau, delta;
    2211             : 
    2212             :   // Determine which region the point is in
    2213             :   unsigned int region =
    2214         236 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2215         236 :   switch (region)
    2216             :   {
    2217         227 :     case 1:
    2218         227 :       pi = pressure / _p_star[0];
    2219         227 :       tau = _T_star[0] / temperature;
    2220         252 :       speed2 = _Rw * temperature * Utility::pow<2>(dgamma1_dpi(pi, tau)) /
    2221         222 :                (Utility::pow<2>(dgamma1_dpi(pi, tau) - tau * d2gamma1_dpitau(pi, tau)) /
    2222         222 :                     (tau * tau * d2gamma1_dtau2(pi, tau)) -
    2223         222 :                 d2gamma1_dpi2(pi, tau));
    2224         227 :       break;
    2225             : 
    2226           3 :     case 2:
    2227           3 :       pi = pressure / _p_star[1];
    2228           3 :       tau = _T_star[1] / temperature;
    2229           3 :       speed2 = _Rw * temperature * Utility::pow<2>(pi * dgamma2_dpi(pi, tau)) /
    2230           3 :                ((-pi * pi * d2gamma2_dpi2(pi, tau)) +
    2231           3 :                 Utility::pow<2>(pi * dgamma2_dpi(pi, tau) - tau * pi * d2gamma2_dpitau(pi, tau)) /
    2232           3 :                     (tau * tau * d2gamma2_dtau2(pi, tau)));
    2233           3 :       break;
    2234             : 
    2235           3 :     case 3:
    2236             :     {
    2237             :       // Calculate density first, then use that in Helmholtz free energy
    2238           3 :       T density3 = densityRegion3(pressure, temperature);
    2239           3 :       delta = density3 / _rho_critical;
    2240           3 :       tau = _T_star[2] / temperature;
    2241           3 :       speed2 =
    2242           3 :           _Rw * temperature *
    2243           3 :           (2.0 * delta * dphi3_ddelta(delta, tau) + delta * delta * d2phi3_ddelta2(delta, tau) -
    2244           3 :            Utility::pow<2>(delta * dphi3_ddelta(delta, tau) -
    2245           3 :                            delta * tau * d2phi3_ddeltatau(delta, tau)) /
    2246           3 :                (tau * tau * d2phi3_dtau2(delta, tau)));
    2247           3 :       break;
    2248             :     }
    2249             : 
    2250           3 :     case 5:
    2251           3 :       pi = pressure / _p_star[4];
    2252           3 :       tau = _T_star[4] / temperature;
    2253           3 :       speed2 = _Rw * temperature * Utility::pow<2>(pi * dgamma5_dpi(pi, tau)) /
    2254           3 :                ((-pi * pi * d2gamma5_dpi2(pi, tau)) +
    2255           3 :                 Utility::pow<2>(pi * dgamma5_dpi(pi, tau) - tau * pi * d2gamma5_dpitau(pi, tau)) /
    2256           3 :                     (tau * tau * d2gamma5_dtau2(pi, tau)));
    2257           3 :       break;
    2258             : 
    2259           0 :     default:
    2260           0 :       mooseError("inRegion() has given an incorrect region");
    2261             :   }
    2262             : 
    2263         236 :   return std::sqrt(speed2);
    2264             : }
    2265             : 
    2266             : template <typename T>
    2267             : T
    2268         671 : Water97FluidProperties::cp_from_p_T_template(const T & pressure, const T & temperature) const
    2269             : {
    2270             :   T specific_heat, pi, tau, delta;
    2271             : 
    2272             :   // Determine which region the point is in
    2273             :   unsigned int region =
    2274         671 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2275         671 :   switch (region)
    2276             :   {
    2277         662 :     case 1:
    2278         662 :       pi = pressure / _p_star[0];
    2279         662 :       tau = _T_star[0] / temperature;
    2280         667 :       specific_heat = -_Rw * tau * tau * d2gamma1_dtau2(pi, tau);
    2281         662 :       break;
    2282             : 
    2283           3 :     case 2:
    2284           3 :       pi = pressure / _p_star[1];
    2285           3 :       tau = _T_star[1] / temperature;
    2286           3 :       specific_heat = -_Rw * tau * tau * d2gamma2_dtau2(pi, tau);
    2287           3 :       break;
    2288             : 
    2289           3 :     case 3:
    2290             :     {
    2291             :       // Calculate density first, then use that in Helmholtz free energy
    2292           3 :       T density3 = densityRegion3(pressure, temperature);
    2293           3 :       delta = density3 / _rho_critical;
    2294           3 :       tau = _T_star[2] / temperature;
    2295           3 :       specific_heat =
    2296           3 :           _Rw *
    2297           3 :           (-tau * tau * d2phi3_dtau2(delta, tau) +
    2298           3 :            (delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau)) *
    2299           3 :                (delta * dphi3_ddelta(delta, tau) - delta * tau * d2phi3_ddeltatau(delta, tau)) /
    2300           3 :                (2.0 * delta * dphi3_ddelta(delta, tau) +
    2301           3 :                 delta * delta * d2phi3_ddelta2(delta, tau)));
    2302           3 :       break;
    2303             :     }
    2304             : 
    2305           3 :     case 5:
    2306           3 :       pi = pressure / _p_star[4];
    2307           3 :       tau = _T_star[4] / temperature;
    2308           3 :       specific_heat = -_Rw * tau * tau * d2gamma5_dtau2(pi, tau);
    2309           3 :       break;
    2310             : 
    2311           0 :     default:
    2312           0 :       mooseError("inRegion() has given an incorrect region");
    2313             :   }
    2314         671 :   return specific_heat;
    2315             : }
    2316             : 
    2317             : template <typename T>
    2318             : T
    2319         228 : Water97FluidProperties::cv_from_p_T_template(const T & pressure, const T & temperature) const
    2320             : {
    2321             :   T specific_heat, pi, tau, delta;
    2322             : 
    2323             :   // Determine which region the point is in
    2324             :   unsigned int region =
    2325         228 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2326         228 :   switch (region)
    2327             :   {
    2328         225 :     case 1:
    2329         225 :       pi = pressure / _p_star[0];
    2330         225 :       tau = _T_star[0] / temperature;
    2331         240 :       specific_heat =
    2332         225 :           _Rw * (-tau * tau * d2gamma1_dtau2(pi, tau) +
    2333         220 :                  Utility::pow<2>(dgamma1_dpi(pi, tau) - tau * d2gamma1_dpitau(pi, tau)) /
    2334         220 :                      d2gamma1_dpi2(pi, tau));
    2335         225 :       break;
    2336             : 
    2337           1 :     case 2:
    2338           1 :       pi = pressure / _p_star[1];
    2339           1 :       tau = _T_star[1] / temperature;
    2340           1 :       specific_heat =
    2341           1 :           _Rw * (-tau * tau * d2gamma2_dtau2(pi, tau) +
    2342           1 :                  Utility::pow<2>(dgamma2_dpi(pi, tau) - tau * d2gamma2_dpitau(pi, tau)) /
    2343           1 :                      d2gamma2_dpi2(pi, tau));
    2344           1 :       break;
    2345             : 
    2346           1 :     case 3:
    2347             :     {
    2348             :       // Calculate density first, then use that in Helmholtz free energy
    2349           1 :       T density3 = densityRegion3(pressure, temperature);
    2350           1 :       delta = density3 / _rho_critical;
    2351           1 :       tau = _T_star[2] / temperature;
    2352           1 :       specific_heat = _Rw * (-tau * tau * d2phi3_dtau2(delta, tau));
    2353           1 :       break;
    2354             :     }
    2355             : 
    2356           1 :     case 5:
    2357           1 :       pi = pressure / _p_star[4];
    2358           1 :       tau = _T_star[4] / temperature;
    2359           1 :       specific_heat =
    2360           1 :           _Rw * (-tau * tau * d2gamma5_dtau2(pi, tau) +
    2361           1 :                  Utility::pow<2>(dgamma5_dpi(pi, tau) - tau * d2gamma5_dpitau(pi, tau)) /
    2362           1 :                      d2gamma5_dpi2(pi, tau));
    2363           1 :       break;
    2364             : 
    2365           0 :     default:
    2366           0 :       mooseError("inRegion() has given an incorrect region");
    2367             :   }
    2368         228 :   return specific_heat;
    2369             : }
    2370             : 
    2371             : template <typename T>
    2372             : T
    2373         241 : Water97FluidProperties::k_from_rho_T_template(const T & density, const T & temperature) const
    2374             : {
    2375             :   // Scale the density and temperature. Note that the scales are slightly
    2376             :   // different to the critical values used in IAPWS-IF97
    2377         241 :   T Tbar = temperature / 647.26;
    2378         241 :   T rhobar = density / 317.7;
    2379             : 
    2380             :   // Ideal gas component
    2381           8 :   T sum0 = 0.0;
    2382             : 
    2383        1205 :   for (std::size_t i = 0; i < _k_a.size(); ++i)
    2384         996 :     sum0 += _k_a[i] * MathUtils::pow(Tbar, i);
    2385             : 
    2386         249 :   T lambda0 = std::sqrt(Tbar) * sum0;
    2387             : 
    2388             :   // The contribution due to finite density
    2389         257 :   T lambda1 = -0.39707 + 0.400302 * rhobar +
    2390         257 :               1.06 * std::exp(-0.171587 * Utility::pow<2>(rhobar + 2.392190));
    2391             : 
    2392             :   // Critical enhancement
    2393         241 :   T DeltaT = std::abs(Tbar - 1.0) + 0.00308976;
    2394         257 :   T Q = 2.0 + 0.0822994 / std::pow(DeltaT, 0.6);
    2395         249 :   T S = (Tbar >= 1.0 ? 1.0 / DeltaT : 10.0932 / std::pow(DeltaT, 0.6));
    2396             : 
    2397         273 :   T lambda2 = (0.0701309 / Utility::pow<10>(Tbar) + 0.011852) * std::pow(rhobar, 1.8) *
    2398         249 :                   std::exp(0.642857 * (1.0 - std::pow(rhobar, 2.8))) +
    2399         241 :               0.00169937 * S * std::pow(rhobar, Q) *
    2400         257 :                   std::exp((Q / (1.0 + Q)) * (1.0 - std::pow(rhobar, 1.0 + Q))) -
    2401         265 :               1.02 * std::exp(-4.11717 * std::pow(Tbar, 1.5) - 6.17937 / Utility::pow<5>(rhobar));
    2402             : 
    2403         241 :   return lambda0 + lambda1 + lambda2;
    2404             : }
    2405             : 
    2406             : template <typename T>
    2407             : T
    2408         241 : Water97FluidProperties::k_from_p_T_template(const T & pressure, const T & temperature) const
    2409             : {
    2410         241 :   T rho = this->rho_from_p_T(pressure, temperature);
    2411         241 :   return this->k_from_rho_T_template(rho, temperature);
    2412             : }
    2413             : 
    2414             : template <typename T>
    2415             : T
    2416           8 : Water97FluidProperties::k_from_v_e_template(const T & v, const T & e) const
    2417             : {
    2418           8 :   const auto [p, temperature] = p_T_from_v_e(v, e);
    2419           8 :   return k_from_p_T_template(p, temperature);
    2420             : }
    2421             : 
    2422             : template <typename T>
    2423             : T
    2424          22 : Water97FluidProperties::p_from_v_e_template(const T & v, const T & e) const
    2425             : {
    2426          22 :   const auto rho = 1 / v;
    2427             :   // For NewtonSolve
    2428             :   // y = e
    2429             :   // x = rho
    2430             :   // z = p
    2431          66 :   auto lambda = [this](const T & rho, const T & p, T & e, T & de_drho, T & de_dp)
    2432          44 :   { e_from_p_rho(p, rho, e, de_dp, de_drho); };
    2433           7 :   return FluidPropertiesUtils::NewtonSolve(
    2434          22 :              rho, e, _p_initial_guess, _tolerance, lambda, name() + "::p_from_v_e")
    2435          37 :       .first;
    2436             : }
    2437             : 
    2438             : template <typename T>
    2439             : T
    2440          52 : Water97FluidProperties::e_from_p_rho_template(const T & p, const T & rho) const
    2441             : {
    2442          52 :   const auto temperature = T_drhodT_from_p_rho(p, rho).first;
    2443          52 :   return e_from_p_T_template(p, temperature);
    2444             : }
    2445             : 
    2446             : template <typename T>
    2447             : void
    2448          31 : Water97FluidProperties::e_from_p_rho_template(
    2449             :     const T & p, const T & rho, T & e, T & de_dp, T & de_drho) const
    2450             : {
    2451          46 :   auto functor = [this](const auto & p, const auto & rho)
    2452          46 :   { return this->e_from_p_rho_template(p, rho); };
    2453             : 
    2454          46 :   xyDerivatives(p, rho, e, de_dp, de_drho, functor);
    2455          31 : }
    2456             : 
    2457             : template <typename T>
    2458             : T
    2459      220266 : Water97FluidProperties::mu_from_rho_T_template(const T & density, const T & temperature) const
    2460             : {
    2461           5 :   const T mu_star = 1.e-6;
    2462      220266 :   const T rhobar = density / _rho_critical;
    2463      220266 :   const T Tbar = temperature / _T_critical;
    2464             : 
    2465             :   // Viscosity in limit of zero density
    2466           5 :   T sum0 = 0.0;
    2467     1101330 :   for (std::size_t i = 0; i < _mu_H0.size(); ++i)
    2468      881084 :     sum0 += _mu_H0[i] / MathUtils::pow(Tbar, i);
    2469             : 
    2470      220266 :   const T mu0 = 100.0 * std::sqrt(Tbar) / sum0;
    2471             : 
    2472             :   // Residual component due to finite density
    2473           5 :   T sum1 = 0.0;
    2474     1541862 :   for (unsigned int i = 0; i < 6; ++i)
    2475             :   {
    2476     1321626 :     const T fact = MathUtils::pow(1.0 / Tbar - 1.0, i);
    2477    10572768 :     for (unsigned int j = 0; j < 7; ++j)
    2478     9251382 :       sum1 += fact * _mu_Hij[i][j] * MathUtils::pow(rhobar - 1.0, j);
    2479             :   }
    2480             : 
    2481      220266 :   const T mu1 = std::exp(rhobar * sum1);
    2482             : 
    2483             :   // The water viscosity (in Pa.s) is then given by
    2484      220266 :   return mu_star * mu0 * mu1;
    2485             : }
    2486             : 
    2487             : template <typename T>
    2488             : std::pair<T, T>
    2489           4 : Water97FluidProperties::p_T_from_v_h(const T & v, const T & h) const
    2490             : {
    2491             :   T pressure, temperature;
    2492             :   bool conversion_succeeded;
    2493           4 :   p_T_from_v_h(
    2494           4 :       v, h, _p_initial_guess, _T_initial_guess, pressure, temperature, conversion_succeeded);
    2495           4 :   return {std::move(pressure), std::move(temperature)};
    2496             : }
    2497             : 
    2498             : template <typename T>
    2499             : T
    2500      220883 : Water97FluidProperties::h_from_p_T_template(const T & pressure, const T & temperature) const
    2501             : {
    2502             :   T enthalpy, pi, tau, delta;
    2503             : 
    2504             :   // Determine which region the point is in
    2505             :   unsigned int region =
    2506      220883 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2507      220883 :   switch (region)
    2508             :   {
    2509      217751 :     case 1:
    2510      217751 :       pi = pressure / _p_star[0];
    2511      217751 :       tau = _T_star[0] / temperature;
    2512      217767 :       enthalpy = _Rw * _T_star[0] * dgamma1_dtau(pi, tau);
    2513      217751 :       break;
    2514             : 
    2515        3111 :     case 2:
    2516        3111 :       pi = pressure / _p_star[1];
    2517        3111 :       tau = _T_star[1] / temperature;
    2518        3111 :       enthalpy = _Rw * _T_star[1] * dgamma2_dtau(pi, tau);
    2519        3111 :       break;
    2520             : 
    2521          12 :     case 3:
    2522             :     {
    2523             :       // Calculate density first, then use that in Helmholtz free energy
    2524          12 :       T density3 = densityRegion3(pressure, temperature);
    2525          12 :       delta = density3 / _rho_critical;
    2526          12 :       tau = _T_star[2] / temperature;
    2527          13 :       enthalpy =
    2528          12 :           _Rw * temperature * (tau * dphi3_dtau(delta, tau) + delta * dphi3_ddelta(delta, tau));
    2529          11 :       break;
    2530             :     }
    2531             : 
    2532           9 :     case 5:
    2533           9 :       pi = pressure / _p_star[4];
    2534           9 :       tau = _T_star[4] / temperature;
    2535           9 :       enthalpy = _Rw * _T_star[4] * dgamma5_dtau(pi, tau);
    2536           9 :       break;
    2537             : 
    2538           0 :     default:
    2539           0 :       mooseError("Water97FluidProperties::inRegion has given an incorrect region");
    2540             :   }
    2541      220883 :   return enthalpy;
    2542             : }
    2543             : 
    2544             : template <typename T>
    2545             : void
    2546          15 : Water97FluidProperties::h_from_p_T_template(
    2547             :     const T & pressure, const T & temperature, T & h, T & dh_dp, T & dh_dT) const
    2548             : {
    2549          30 :   auto functor = [this](const auto & pressure, const auto & temperature)
    2550          30 :   { return this->h_from_p_T_template(pressure, temperature); };
    2551             : 
    2552          30 :   xyDerivatives(pressure, temperature, h, dh_dp, dh_dT, functor);
    2553          15 : }
    2554             : 
    2555             : template <typename T>
    2556             : T
    2557      441530 : Water97FluidProperties::rho_from_p_T_template(const T & pressure, const T & temperature) const
    2558             : {
    2559         124 :   T density, pi, tau;
    2560             : 
    2561             :   // Determine which region the point is in
    2562             :   unsigned int region =
    2563      441530 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2564             : 
    2565      441530 :   switch (region)
    2566             :   {
    2567      435381 :     case 1:
    2568      435381 :       pi = pressure / _p_star[0];
    2569      435381 :       tau = _T_star[0] / temperature;
    2570      435661 :       density = pressure / (pi * _Rw * temperature * dgamma1_dpi(pi, tau));
    2571      435381 :       break;
    2572             : 
    2573        6129 :     case 2:
    2574        6129 :       pi = pressure / _p_star[1];
    2575        6129 :       tau = _T_star[1] / temperature;
    2576        6129 :       density = pressure / (pi * _Rw * temperature * dgamma2_dpi(pi, tau));
    2577        6129 :       break;
    2578             : 
    2579           9 :     case 3:
    2580           9 :       density = densityRegion3(pressure, temperature);
    2581           1 :       break;
    2582             : 
    2583          11 :     case 5:
    2584          11 :       pi = pressure / _p_star[4];
    2585          11 :       tau = _T_star[4] / temperature;
    2586          11 :       density = pressure / (pi * _Rw * temperature * dgamma5_dpi(pi, tau));
    2587          11 :       break;
    2588             : 
    2589           0 :     default:
    2590           0 :       mooseError("inRegion() has given an incorrect region");
    2591             :   }
    2592      441530 :   return density;
    2593             : }
    2594             : 
    2595             : template <typename T>
    2596             : T
    2597      220247 : Water97FluidProperties::mu_from_p_T_template(const T & pressure, const T & temperature) const
    2598             : {
    2599      220247 :   T rho = this->rho_from_p_T_template(pressure, temperature);
    2600      220247 :   return this->mu_from_rho_T_template(rho, temperature);
    2601             : }
    2602             : 
    2603             : template <typename T>
    2604             : T
    2605         275 : Water97FluidProperties::s_from_p_T_template(const T & pressure, const T & temperature) const
    2606             : {
    2607             :   T entropy, pi, tau, density3, delta;
    2608             : 
    2609             :   // Determine which region the point is in
    2610             :   unsigned int region =
    2611         275 :       inRegion(MetaPhysicL::raw_value(pressure), MetaPhysicL::raw_value(temperature));
    2612         275 :   switch (region)
    2613             :   {
    2614         242 :     case 1:
    2615         242 :       pi = pressure / _p_star[0];
    2616         242 :       tau = _T_star[0] / temperature;
    2617         256 :       entropy = _Rw * (tau * dgamma1_dtau(pi, tau) - gamma1(pi, tau));
    2618         242 :       break;
    2619             : 
    2620          12 :     case 2:
    2621          12 :       pi = pressure / _p_star[1];
    2622          12 :       tau = _T_star[1] / temperature;
    2623          13 :       entropy = _Rw * (tau * dgamma2_dtau(pi, tau) - gamma2(pi, tau));
    2624          12 :       break;
    2625             : 
    2626          12 :     case 3:
    2627             :       // Calculate density first, then use that in Helmholtz free energy
    2628          12 :       density3 = densityRegion3(pressure, temperature);
    2629          12 :       delta = density3 / _rho_critical;
    2630          12 :       tau = _T_star[2] / temperature;
    2631          13 :       entropy = _Rw * (tau * dphi3_dtau(delta, tau) - phi3(delta, tau));
    2632          12 :       break;
    2633             : 
    2634           9 :     case 5:
    2635           9 :       pi = pressure / _p_star[4];
    2636           9 :       tau = _T_star[4] / temperature;
    2637          10 :       entropy = _Rw * (tau * dgamma5_dtau(pi, tau) - gamma5(pi, tau));
    2638           9 :       break;
    2639             : 
    2640           0 :     default:
    2641           0 :       mooseError("inRegion() has given an incorrect region");
    2642             :   }
    2643         275 :   return entropy;
    2644             : }
    2645             : 
    2646             : template <typename T>
    2647             : void
    2648           3 : Water97FluidProperties::s_from_p_T_template(
    2649             :     const T & pressure, const T & temperature, T & s, T & ds_dp, T & ds_dT) const
    2650             : {
    2651          10 :   auto functor = [this](const auto & pressure, const auto & temperature)
    2652          10 :   { return this->s_from_p_T_template(pressure, temperature); };
    2653             : 
    2654          10 :   xyDerivatives(pressure, temperature, s, ds_dp, ds_dT, functor);
    2655           3 : }

Generated by: LCOV version 1.14