/** This class contains equations about the properties of water.

    It is primarily based on The International Association for the Properties
    of Water and Steam (IAPWS, http://iapws.org/) document IAPWS R6-95(2018),
    a PDF of which is usually available from:

    http://www.iapws.org/relguide/IAPWS-95.html

    Specifically, the file
    http://www.iapws.org/relguide/IAPWS95-2018.pdf

    All references to equation and table numbers are the equation numbers in
    this document.

    The development of this document is discussed in
    https://aip.scitation.org/doi/10.1063/1.1461829

    There is an online calculator that can be used to validate these equations
    at:
    http://twt.mpei.ac.ru/mcs/worksheets/iapws/IAPWS95.xmcd

    IAPWS also produces a simpler formulation "for industry", called IAPWS-IF97
    which is described at:
    http://www.iapws.org/relguide/IF97-Rev.html

    That page contains a PDF of IAPWS-IF97 and some "backward equations" that
    make it easier to calculate without iterating inverse solutions.  These
    may be used to find good "first guesses" for inverse solutions.

    Simpler equations for vapor pressure and densities at the saturation point
    are found here:
    http://www.iapws.org/relguide/Supp-sat.html

    especially in the PDF file:
    http://www.iapws.org/relguide/Supp-sat.pdf

    These can be used to obtain good initial guesses for finding
    higher-accuracy results in the full IAPWS95 model.

    A scientific paper describing the derivation of the saturation properties
    can be found at:
    https://aip.scitation.org/doi/10.1063/1.555926

    All of the methods in this class are class-level methods so you don't
    need an instance of the Water class to call them.  Just call them like:

    Water.boilingPoint[1 atmosphere]
*/
class Water
{
   /** Constants as defined in IAPWS RS-95(2018).  These may not exactly match
       the current best-known values of these constants in the SI, but they are
       used as defined there for purposes of matching its output exactly.
   */

   /** The temperature of the critical point of water. (eq.1) */
   class var Tc = 647.096 K

   /** The critical density of water. (eq.2) */
   class var rhoc = 322 kg m^-3

   /** The specific gas constant.  Due to the use of the specific gas constant,
       Eq. (4) corresponds to a mass-based formulation.  In order to convert
       values of specific properties to molar properties, an appropriate value
       for the molar mass should be used  (eq. 3) */
   class var R = 0.46151805 kJ kg^-1 K^-1

   /** The critical pressure of water from IAPWS SR1-86(1992) */
   class var pc = 22.064 MPa

   /** The triple point of water. */
   class var TriplePoint = 273.16 K    

   /** IAPWS doesn't directly define the molar mass of water directly but
       indicates that the best-known values of several of these constants have
       changed since the model was published.  It sort of waves its hand at
       George S. Kell, Effects of isotopic composition, temperature, pressure,
       and dissolved gases on the density of liquid water Journal of Physical
       and Chemical Reference Data 6, 1109 (1977);
       https://doi.org/10.1063/1.555561

       to hint that it contains a suitable definition of the molar mass.  Its
       masses and isotope ratios are taken from "Atomic Weights of the
       Elements 1973", Pure Appl. Chem 37, 589-603 (1974) Of course, none of
       this matches the current definition of Avogadro's constant, the gas
       constant, etc.  Anyway, this is the version of the molar mass of water
       defined in that paper, and the one that should probably be used with
       this model and its definition of R.  See the notes for R above.  */
   class var molarMassH2O = 18.015242 g/mol
       
       
   /** Calculates the pressure at the specified density and temperature.
       
       p = ρ^2 (df/dρ)_T
       
       See Table 3 for this equation.
   */
   class pressure[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      return ρ R T (1 + δ dφrδDimensioned[ρ, T])
   }


   /** Calculates the compressibility factor Z = p / (ρ R T) */
   class compressibilityFactor[ρ is mass_density, T is temperature] :=
   {
      return pressure[ρ, T] / (ρ R T)
   }

   
   /** Calculates the specific internal energy of the system at the specified
       density and temperature.  Multiply this by mass to get the energy.

       u = f - T (df/dT)_rho

       returns energy / mass
   */
   class specificInternalEnergy[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return τ R T dφdτ[δ, τ]
   }


   /** Calculates the specific entropy

       s = -(df/dT)_rho

       Multiply by mass to get the total entropy.
   
       This has units of specific heat capacity (e.g. cal/g/degC)
   */
   class specificEntropy[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return R (τ dφdτ[δ, τ] - φ[δ, τ])
   }
   
   
   /** Calculates the specific enthalpy

       h = f - T(df/dT)_rho + ρ(df/dρ)_T
       Multiply by mass to get the total enthalpy.
   
       This has units of energy/mass
   */
   class specificEnthalpy[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return R T (1 + τ dφdτ[δ, τ] + δ dφrδ[δ, τ])
   }

   /** Calculates the specific isochoric heat capacity

       c_v = (du/dT)_rho

       This has units of specific heat capacity (e.g. cal/g/degC)
   */
   class specificIsochoricHeatCapacity[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return R (-τ^2 d2φτ[δ, τ])
   }
   
   
   /** Calculates the specific isobaric heat capacity

       c_p = (dh/dT)_p
   
       This has units of specific heat capacity (e.g. cal/g/degC)
   */
   class specificIsobaricHeatCapacity[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return R (-τ^2 d2φτ[δ, τ] + (1 + δ dφrδ[δ, τ] - δ τ d2φrδτ[δ, τ])^2 / (1 + 2 δ dφrδ[δ, τ] + δ^2 d2φrδ[δ, τ]))
   }

   
   /** Calculates the speed of sound.

       w = (dp/dρ)^(1/2)
   
       This has units of velocity
   */
   class speedOfSound[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return sqrt[R T (1 + 2 δ dφrδ[δ, τ] + δ^2 d2φrδ[δ,τ] - (1 + δ dφrδ[δ, τ] - δ τ d2φrδτ[δ, τ])^2/(τ^2 (d2φτ[δ, τ])))]
   }

   
   /** The Joule-Thompson coefficient

       mu = (dT/dρ)_h

       Results have units of temperature / pressure
   */
   class JouleThompsonCoefficient[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return (-(δ dφrδ[δ, τ] + δ^2 d2φrδ[δ,τ] + δ τ d2φrδτ[δ,τ]) / ((1 + δ dφrδ[δ, τ] - δ τ d2φrδτ[δ,τ])^2 - τ^2 d2φτ[δ,τ] (1 + 2 δ dφrδ[δ, τ] + δ^2 d2φrδ[δ, τ]))) / (R ρ)
   }

   /** The isentropic temperature-pressure coefficent

       βₛ = (dT/dp)_s

       Results have units of temperature / pressure
   */
   class isentropicTemperaturePressureCoefficient[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return  (1 + δ dφrδ[δ, τ] - δ τ d2φrδτ[δ, τ]) / ((1 + δ dφrδ[δ, τ] - δ τ d2φrδτ[δ,τ])^2 - τ^2 d2φτ[δ,τ] (1 + 2 δ dφrδ[δ, τ] + δ^2 d2φrδ[δ, τ])) / (R ρ)
   }

   
   /** The isothermal throttling coefficient

       delta_T = (dh/dp)_T

       This has units of inverse density (e.g. m^3 / kg)
   */
   class isothermalThrottlingCoefficient[ρ is mass_density, T is temperature] :=
   {
      δ = ρ / rhoc
      τ = Tc / T
      return (1-(1 + δ dφrδ[δ, τ] - δ τ d2φrδτ[δ, τ])/ (1 + 2 δ dφrδ[δ, τ] + δ^2 d2φrδ[δ,τ])) / ρ
   }


   /** The second virial coefficient B(T).
   
       "Even at low density, real gases don't quite obey the ideal gas law.
        A systematic way to account for deviations from ideal behavior is the
        virial expansion,

        P V = n R T (1 + B(T)/(V/n) + C(T)/(V/n)^2 + ...)

        where the functions B(T), C(T), and so on are called the "virial
        coefficents".  When the density of the gas is fairly low, so that the
        volume per mole is large, each term in the series is much smaller than
        the one before.  In many situations it's sufficient to omit the third
        term and concentrate on the second, whose coefficent B(T) is called the
        second virial coefficient (the first coefficient being 1.)"
        --Daniel V. Schroeder, An Introduction to Thermal Physics, p.9

        B(τ) rhoc = lim[φrδ(δ, τ), δ->0]

        Units are m^3 kg^-1

        This actually returns the *specific* version of B.
        To get the dimensionless units that will work in the equation above,
        you need to multiply by Water.molarMassH2O
   */
   class B[T is temperature] :=
   {
      τ = Tc / T
      sum = 0
      
      for i = [1,2,3,8,9,10,23]
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ τ^tᵢ
      }

      for i = 55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         θᵢ = Aᵢ + 1 - τ       // See footnote b, p.12
         sum = sum + nᵢ (θᵢ^2 + Bᵢ)^bᵢ exp[-Cᵢ - Dᵢ(τ-1)^2]
      }
      
      return sum/rhoc
   }
   
   /** The third virial coefficient C(T).  See the discussion for the second
       virial coefficient B[T] above.

       Units are m^6 kg^-2

       This actually returns the *specific* version of C.
       To get the dimensionless units for the equation documented in B[T],
       you need to multiply by Water.molarMassH2O
   */
   class C[T is temperature] :=
   {
      τ = Tc / T
      sum = 0
      
      for i = [4,5,11,12,24,25,26]
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + 2 nᵢ τ^tᵢ
      }

      for i = 8 to 10
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum - 2 nᵢ τ^tᵢ
      }
      
      for i = 55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         θᵢ = Aᵢ + 1 - τ       // See footnote b, p.12
         sum = sum + 4 nᵢ (Cᵢ (θᵢ^2 + Bᵢ) - bᵢ (Aᵢ θᵢ / βᵢ + Bᵢ aᵢ))* (θᵢ^2 + Bᵢ)^(bᵢ-1) exp[-Cᵢ - Dᵢ(τ-1)^2]
      }
      
      return sum/rhoc^2
   }
   
   /** Calculates the specific internal energy of the system at the specified
       atmospheric pressure equals the saturated vapor pressure.   You can use
       StandardAtmosphereTest.frink which calls StandardAtmosphere.frink
       to find the approximate atmospheric pressure for any altitude.
   */
   class boilingPoint[ps is pressure] :=
   {
      return saturatedTemperature[ps]
   }

   
   /** Returns the approximate saturated vapor pressure of water from
       IAPWS SR1-86(1992).  This does not use the full IAPWS95 model, but it
       can be used to get a very good first guess.   This is Section 3 */
   class saturatedVaporPressure[T is temperature] :=
   {
      θ = T/Tc
      τ = 1 - θ
      p = pc exp[Tc/T (-7.85951783 τ +
                        1.84408259 τ^(3/2) +
                       -11.7866497 τ^3 +
                        22.6807411 τ^(7/2) +
                       -15.9618719 τ^4 +
                        1.80122502 τ^(15/2))]
      return p
   }

   
   /** Returns the approximate saturated vapor pressure of water from
       IAPWS IF-97.  This does not use the full IAPWS95 model, but it
       can be used to get a very good first guess.   This is Section 8 and
       8.1 of IAPWS IF-97 */
   class saturatedVaporPressure97[Ts is temperature] :=
   {
      Tstar = 1 K
      Θ = Ts / Tstar + Table34@9 / ((Ts/Tstar) - Table34@10)  // eq. 29b
      A = Θ^2 + Table34@1 Θ + Table34@2
      B = Table34@3 Θ^2 + Table34@4 Θ + Table34@5
      C = Table34@6 Θ^2 + Table34@7 Θ + Table34@8
      p = 1 MPa ((2 C) / (-B + (B^2 - 4 A C)^(1/2)))^4   // eq. 30
      return p
   }

   
   /** Returns the approximate saturated vapor temperature of water from
       IAPWS97.  This does not use the full IAPWS95 model, but it
       can be used to get a very good first guess.   This is Section 8.2 of
       IAPWS97. */
   class saturatedTemperature[ps is pressure] :=
   {
      // Equation 31
      pstar = 1 MPa
      β = (ps / pstar)^(1/4)  // Eq. 29a
      
      E = β^2 + Table34@3 β + Table34@6
      F = Table34@1 β^2 + Table34@4 β + Table34@7
      G = Table34@2 β^2 + Table34@5 β + Table34@8
      D = (2 G)/(-F - (F^2 - 4 E G)^(1/2))
      
      Ts = 1 K ((Table34@10 + D - ((Table34@10 + D)^2 - 4(Table34@9 + Table34@10 D))^(1/2)) / 2)
      return Ts
   }

   
   /** Returns the approximate saturated liquid density of water from
       IAPWS SR1-86(1992).  This does not use the full IAPWS95 model, but it
       can be used to get a very good first guess.   This is Section 4.1 */
   class saturatedLiquidDensity[T is temperature] :=
   {
      θ = T/Tc
      τ = 1-θ
      ρ = rhoc (1 + 1.99274064 τ^(1/3)    + // b1
                      1.09965342 τ^(2/3)    + // b2
                    -0.510839303 τ^(5/3)    + // b3
                     -1.75493479 τ^(16/3)   + // b4 
                     -45.5170352 τ^(43/3)   + // b5
                   -6.74694450e5 τ^(110/3))   // b6
      return ρ
   }

   
   /** Returns the approximate saturated vapor density of water from
       IAPWS SR1-86(1992).  This does not use the full IAPWS95 model, but it
       can be used to get a very good first guess.   This is Section 4.2 */
   class saturatedVaporDensity[T is temperature] :=
   {
      θ = T/Tc
      τ = 1-θ
      ρ = rhoc exp[1 + -2.03150240 τ^(2/6)    + // c1
                         -2.68302940 τ^(4/6)    + // c2
                         -5.38626492 τ^(8/6)    + // c3
                         -17.2991605 τ^(18/6)   + // c4 
                         -44.7586581 τ^(37/6)   + // c5
                         -63.9201063 τ^(71/6)]    // c6
      return ρ
   }

   /** The specific Helmholtz free energy f (also sometimes called A).
       The IAPWS formulation is implemented as a fundamental equation for the
       specific Helmholtz free energy f.
       
       This equation is expressed in dimensionless form, φ = f/(R T),
       or, conversely f = φ R T

       (see eq. 4 and following)
   */
   class f[ρ is mass_density, T is temperature, debug=false] :=
   {
      return R T phiDimensioned[ρ, T, debug]
   }

   
   /** The formulation of the φ equation is implemented as a dimensionless
       version φ[δ, τ] where
       δ = ρ / rhoc
       τ = Tc / T
   
       This turns a dimensionally-correct call with density and temperature
       into the internal dimensionless form.  (see eq.4 and following)
   */
   class phiDimensioned[ρ is mass_density, T is temperature, debug=false] :=
   {
      return φ[ρ / rhoc, Tc / T, debug]
   }

   
   /** The dimensionless φ function φ[δ, τ] is separated into two
       parts, an ideal-gas part φ0[δ, τ] and a residual part
       φr[δ, τ]   (eq. 4)
   */
   class φ[δ is dimensionless, τ is dimensionless, debug=false] :=
   {
      φ0 = φ0[δ, τ]
      φr = φr[δ, τ]

      if debug
         println["φ0=$φ0\nφr=$φr\n"]
      
      return φ0 + φr
   }


   /** Calculates the first derivative of φ with respect to τ.  This
       includes both the ideal gas and residual parts. */
   class dφdτ[δ is dimensionless, τ is dimensionless] :=
   {
      return dφ0dτ[δ, τ] + dφrτ[δ, τ]
   }

   
   /** Calculates the second derivative of φ with respect to τ.  This
       includes both the ideal gas and residual parts. */
   class d2φτ[δ is dimensionless, τ is dimensionless] :=
   {
      return d2φ0dτ[δ, τ] + d2φrτ[δ, τ]
   }

   /** Calculates the second derivative of φ with respect to δ.  This
       includes both the ideal gas and residual parts. */
   class d2φδ[δ is dimensionless, τ is dimensionless] :=
   {
      return d2φ0dδ[δ, τ] + d2φrδ[δ, τ]
   }

   /** The ideal gas part φ0 of the dimensionless Helmholtz free energy.
       This is obtained from an equation for the specific isobaric heat
       capacity in the ideal-gas state developed by Cooper [6]

       This implements eq. 5
   */
   class φ0[δ is dimensionless, τ is dimensionless] :=
   {
      sum = ln[δ] + Table1@1  +       /* n1 0 */
                        Table1@2 τ +    /* n2 0 */
                        Table1@3 ln[τ]  /* n3 0 */

      for i = 4 to 8
      {
         [nᵢ, θᵢ] = Table1@i
         sum = sum + nᵢ ln[1 - exp[-θᵢ τ]]
      }
      
      return sum
   }

   
   /** The first partial derivative of the ideal gas part of the dimensionless
       Helmholtz free energy with respect to δ. */
   class dφ0dδ[δ is dimensionless, τ is dimensionless] :=
   {
      return 1/δ
   }

   
   /** The second partial derivative of the ideal gas part of the dimensionless
       Helmholtz free energy with respect to δ. */
   class d2φ0dδ[δ is dimensionless, τ is dimensionless] :=
   {
      return -1 / δ^2
   }

   
   /** The first partial derivative of the ideal gas part of the dimensionless
       Helmholtz free energy with respect to τ. */
   class dφ0dτ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = Table1@2 + Table1@3 / τ
      for i=4 to 8
      {
         [nᵢ, θᵢ] = Table1@i
         sum = sum + nᵢ θᵢ ((1 - exp[-θᵢ τ])^-1 - 1)
      }

      return sum
   }

   
   /** The second partial derivative of the ideal gas part of the dimensionless
       Helmholtz free energy with respect to τ. */
   class d2φ0dτ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = -Table1@3 / τ^2
      for i=4 to 8
      {
         [nᵢ, θᵢ] = Table1@i
         sum = sum - nᵢ θᵢ^2 exp[-θᵢ τ] (1 - exp[-θᵢ τ])^-2
      }

      return sum
   }

   
   /** The second partial derivative of the ideal gas part of the dimensionless
       Helmholtz free energy with respect to δ then with respect to τ. */
   class d2φ0dδdτ[δ is dimensionless, τ is dimensionless] :=
   {
      return 0
   }

   
   /** The residual part φr[δ, τ] of the dimensionless Helmholtz free
       energy.

       This implements eq. 6
   */
   class φr[δ is dimensionless, τ is dimensionless] :=
   {
      sum = 0
      for i=1 to 7
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ
      }

      for i=8 to 51
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ e^(-δ^cᵢ)
      }

      for i=52 to 54
      {
         [cᵢ, dᵢ, tᵢ, nᵢ, αᵢ, βᵢ, θᵢ, εᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ * exp[-αᵢ (δ - εᵢ)^2 - βᵢ(τ - θᵢ)^2]
      }

      for i=55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         Ψ = exp[-Cᵢ (δ-1)^2 - Dᵢ (τ - 1)^2]
         θ = (1 - τ) + Aᵢ ((δ-1)^2)^(1/(2 βᵢ))
         Δ = θ^2 + Bᵢ ((δ-1)^2)^aᵢ
         sum = sum + nᵢ Δ^bᵢ δ Ψ
      }

      return sum
   }

   
   /** The partial derivative of φr[δ, τ] with respect to δ, with
       correct dimensions. */
   class dφrδDimensioned[ρ is mass_density, T is temperature] :=
   {
      return dφrδ[ρ / rhoc, Tc / T]
   }

   
   /** The partial derivative of φr[δ, τ] with respect to δ.

       This implements eq. 6
   */
   class dφrδ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = 0
      for i=1 to 7
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ dᵢ δ^(dᵢ-1) τ^tᵢ
      }

      for i=8 to 51
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ exp[-δ^cᵢ] (δ^(dᵢ-1) τ^tᵢ(dᵢ - cᵢ δ^cᵢ))
      }

      for i=52 to 54
      {
         [cᵢ, dᵢ, tᵢ, nᵢ, αᵢ, βᵢ, θᵢ, εᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ exp[-αᵢ (δ-εᵢ)^2 - βᵢ (τ-θᵢ)^2] (dᵢ/δ - 2 αᵢ (δ - εᵢ))
      }
      
      for i=55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         Ψ = exp[-Cᵢ (δ-1)^2 - Dᵢ (τ - 1)^2]
         θ = (1 - τ) + Aᵢ ((δ-1)^2)^(1/(2 βᵢ))
         Δ = θ^2 + Bᵢ ((δ-1)^2)^aᵢ
         dΨdδ = -2 Cᵢ (δ-1) Ψ
         dΔdδ = (δ-1) (Aᵢ θ 2 / βᵢ ((δ-1)^2)^(1/(2 βᵢ)-1) + 2 Bᵢ aᵢ((δ-1)^2)^(aᵢ-1))
         dΔbidδ = bᵢ Δ^(bᵢ-1) dΔdδ
         sum = sum + nᵢ (Δ^bᵢ (Ψ + δ dΨdδ) + dΔbidδ δ Ψ)
      }

      return sum
   }

   
   /** The second partial derivative of φr[δ, τ] with respect to δ,
       with correct dimensions. */
   class d2φrδDimensioned[ρ is mass_density, T is temperature] :=
   {
      return d2φrδ[ρ / rhoc, Tc / T]
   }

   
   /** The second partial derivative of φr[δ, τ] with respect to δ.
   */
   class d2φrδ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = 0
      for i=1 to 7
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ dᵢ (dᵢ-1) δ^(dᵢ-2) τ^tᵢ
      }

      for i=8 to 51
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ exp[-δ^cᵢ] (δ^(dᵢ-2) τ^tᵢ ((dᵢ - cᵢ δ^cᵢ)(dᵢ-1-cᵢ δ^cᵢ) - cᵢ^2 δ^cᵢ))
      }

      for i=52 to 54
      {
         [cᵢ, dᵢ, tᵢ, nᵢ, αᵢ, βᵢ, θᵢ, εᵢ] = Table2@i
         sum = sum + nᵢ τ^tᵢ exp[-αᵢ (δ-εᵢ)^2 - βᵢ(τ-θᵢ)^2] * (-2 αᵢ δ^dᵢ + 4 αᵢ^2 δ^dᵢ (δ-εᵢ)^2 - 4 dᵢ αᵢ δ^(dᵢ-1) (δ-εᵢ) + dᵢ (dᵢ-1) δ^(dᵢ-2))
      }
      
      for i=55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         Ψ = exp[-Cᵢ (δ-1)^2 - Dᵢ (τ - 1)^2]
         θ = (1 - τ) + Aᵢ ((δ-1)^2)^(1/(2 βᵢ))
         Δ = θ^2 + Bᵢ ((δ-1)^2)^aᵢ
         dΨdδ = -2 Cᵢ (δ-1) Ψ
         dΔdδ = (δ-1) (Aᵢ θ 2 / βᵢ ((δ-1)^2)^(1/(2 βᵢ)-1) + 2 Bᵢ aᵢ((δ-1)^2)^(aᵢ-1))
         d2Ψdδ = (2 Cᵢ (δ-1)^2 -1) 2 Cᵢ Ψ
         dΔbidδ = bᵢ Δ^(bᵢ-1) dΔdδ
         d2Δdδ = 1 / (δ-1) dΔdδ + (δ-1)^2 (4 Bᵢ aᵢ (aᵢ-1)((δ-1)^2)^(aᵢ-2) + 2 Aᵢ^2 (1/βᵢ)^2 (((δ-1)^2)^(1/(2 βᵢ) - 1))^2 + Aᵢ θ 4/βᵢ (1 / (2 βᵢ) - 1)((δ-1)^2)^(1/(2 βᵢ) - 2))
         d2Δbidδ = bᵢ (Δ^(bᵢ-1) d2Δdδ + (bᵢ - 1) Δ^(bᵢ-2) dΔdδ^2)
         sum = sum + nᵢ (δ^bᵢ (2 dΨdδ + δ d2Ψdδ) + 2 dΔbidδ (Ψ + δ dΨdδ) + d2Δbidδ δ Ψ)
      }

      return sum
   }

   
   /** The partial derivative of φr[δ, τ] with respect to τ, with
       correct dimensions. */
   class dφrτDimensioned[ρ is mass_density, T is temperature] :=
   {
      return dφrτ[ρ / rhoc, Tc / T]
   }

   
   /** The partial derivative of φr[δ, τ] with respect to τ.
   */
   class dφrτ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = 0
      for i=1 to 7
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ tᵢ δ^dᵢ τ^(tᵢ-1)
      }

      for i=8 to 51
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ tᵢ δ^dᵢ τ^(tᵢ-1) exp[-δ^cᵢ]
      }

      for i=52 to 54
      {
         [cᵢ, dᵢ, tᵢ, nᵢ, αᵢ, βᵢ, θᵢ, εᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ exp[-αᵢ (δ-εᵢ)^2 - βᵢ (τ-θᵢ)^2] (tᵢ/τ - 2 βᵢ(τ - θᵢ))
      }
      
      for i=55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         Ψ = exp[-Cᵢ (δ-1)^2 - Dᵢ (τ - 1)^2]
         θ = (1 - τ) + Aᵢ ((δ-1)^2)^(1/(2 βᵢ))
         Δ = θ^2 + Bᵢ ((δ-1)^2)^aᵢ
         dΔbidτ = -2 θ bᵢ Δ^(bᵢ-1)
         dΨdτ = -2 Dᵢ (τ-1) Ψ
         sum = sum + nᵢ Δ (dΔbidτ Ψ + Δ^bᵢ dΨdτ)
      }

      return sum
   }

   
   /** The partial second derivative of φr[δ, τ] with respect to τ,
       with correct dimensions. */
   class d2φrτDimensioned[ρ is mass_density, T is temperature] :=
   {
      return d2φrτ[ρ / rhoc, Tc / T]
   }

   
   /** The partial second derivative of φr[δ, τ] with respect to τ.
   */
   class d2φrτ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = 0
      for i=1 to 7
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ tᵢ (tᵢ-1) δ^dᵢ τ^(tᵢ-2)
      }

      for i=8 to 51
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ tᵢ (tᵢ-1) δ^dᵢ τ^(tᵢ-2) exp[-δ^cᵢ]
      }

      for i=52 to 54
      {
         [cᵢ, dᵢ, tᵢ, nᵢ, αᵢ, βᵢ, θᵢ, εᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ exp[-αᵢ (δ-εᵢ)^2 - βᵢ (τ-θᵢ)^2] ((tᵢ/τ - 2 βᵢ(τ - θᵢ))^2 - tᵢ/τ^2 - 2 βᵢ)
      }
      
      for i=55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         Ψ = exp[-Cᵢ (δ-1)^2 - Dᵢ (τ - 1)^2]
         θ = (1 - τ) + Aᵢ ((δ-1)^2)^(1/(2 βᵢ))
         Δ = θ^2 + Bᵢ ((δ-1)^2)^aᵢ
         dΔbidτ = -2 θ bᵢ Δ^(bᵢ-1)
         dΨdτ = -2 Dᵢ (τ-1) Ψ
         d2Ψdτ = (2 Dᵢ (τ-1)^2 - 1) 2 Dᵢ Ψ
         d2Δbidτ = 2 bᵢ Δ^(bᵢ-1) + 4 θ^2 bᵢ (bᵢ-1) Δ^(bᵢ-2)
         sum = sum + nᵢ δ (d2Δbidτ Ψ + 2 dΔbidτ dΨdτ + Δ^bᵢ d2Ψdτ)
      }

      return sum
   }

   
   /** The partial second derivative of φr[δ, τ] with respect to τ,
       with correct dimensions. */
   class d2φrδτDimensioned[ρ is mass_density, T is temperature] :=
   {
      return d2φrδτ[ρ / rhoc, Tc / T]
   }

   
   /** The partial second derivative of φr[δ, τ] with respect to 
       δ then to τ.
   */
   class d2φrδτ[δ is dimensionless, τ is dimensionless] :=
   {
      sum = 0
      for i=1 to 7
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ dᵢ tᵢ δ^(dᵢ-1) τ^(tᵢ-1)
      }

      for i=8 to 51
      {
         [cᵢ, dᵢ, tᵢ, nᵢ] = Table2@i
         sum = sum + nᵢ tᵢ δ^(dᵢ-1) τ^(tᵢ-1) (dᵢ - cᵢ δ^cᵢ) exp[-δ^cᵢ]
      }

      for i=52 to 54
      {
         [cᵢ, dᵢ, tᵢ, nᵢ, αᵢ, βᵢ, θᵢ, εᵢ] = Table2@i
         sum = sum + nᵢ δ^dᵢ τ^tᵢ exp[-αᵢ (δ-εᵢ)^2 - βᵢ (τ-θᵢ)^2] (dᵢ/δ - 2 αᵢ (δ - εᵢ))(tᵢ/τ - 2 βᵢ(τ - θᵢ))
      }
      
      for i=55 to 56
      {
         [aᵢ, bᵢ, Bᵢ, nᵢ, Cᵢ, Dᵢ, Aᵢ, βᵢ] = Table2@i
         Ψ = exp[-Cᵢ (δ-1)^2 - Dᵢ (τ - 1)^2]
         θ = (1 - τ) + Aᵢ ((δ-1)^2)^(1/(2 βᵢ))
         Δ = θ^2 + Bᵢ ((δ-1)^2)^aᵢ
         dΔbidτ = -2 θ bᵢ Δ^(bᵢ-1)
         dΨdτ = -2 Dᵢ (τ-1) Ψ
         d2Ψdτ = (2 Dᵢ (τ-1)^2 - 1) 2 Dᵢ Ψ
         dΨdδ = -2 Cᵢ (δ-1) Ψ
         d2Δbidτ = 2 bᵢ Δ^(bᵢ-1) + 4 θ^2 bᵢ (bᵢ-1) Δ^(bᵢ-2)
         dΔdδ = (δ-1) (Aᵢ θ 2 / βᵢ ((δ-1)^2)^(1/(2 βᵢ)-1) + 2 Bᵢ aᵢ((δ-1)^2)^(aᵢ-1))
         dΔbidδ = bᵢ Δ^(bᵢ-1) dΔdδ
         d2Ψdδdτ = 4 Cᵢ Dᵢ (δ-1)(τ-1) Ψ
         d2Δbidδdτ = -Aᵢ bᵢ 2 / βᵢ Δ^(bᵢ-1) (δ-1)((δ-1)^2)^(1/(2 βᵢ) - 1) - 2 θ bᵢ (bᵢ-1) Δ^(bᵢ-2) dΔdδ
         sum = sum + nᵢ (Δ^(bᵢ) (dΨdτ + δ d2Ψdδdτ) + δ dΔbidδ dΨdτ + dΔbidτ (Ψ + δ dΨdδ) + d2Δbidδdτ δ Ψ)
      }

      return sum
   }


      /** Numerical values of the coefficients and parameters of the ideal-gas
       part of the dimensionless Helmholtz free energy, Eq. 5.  This is
       table 1 in IAPWS-95.

       The coluns are [ni0, gammai0]
   */
   class var Table1 = [undef,      // No element 0
       -8.3204464837497,
        6.6832105275932,
        3.00632,
       [0.012436, 1.28728967],
       [0.97315,  3.53734222],
       [1.27950,  7.74073708],
       [0.96956,  9.24437796],
       [0.24873, 27.5075105]]
       
   /** Numerical values of the coefficients and parameters of the
       residual part of the dimensionless Helmholtz free energy, Eq. (6) */
   class var Table2 = [undef,    // No element 0
    //[cᵢ, dᵢ,  tᵢ,           nᵢ]    
       [x,  1, -1/2,   0.12533547935523e-1],
       [x,  1, 7/8,    0.78957634722828e1],
       [x,  1,  1,    -0.87803203303561e1],
       [x,  2, 1/2,    0.31802509345418],
       [x,  2, 3/4,   -0.26145533859358],
       [x,  3, 3/8,   -0.78199751687981e-2],
       [x,  4,  1,     0.88089493102134e-2],
       [1,  1,  4,    -0.66856572307965],
       [1,  1,  6,     0.20433810950965],
       [1,  1, 12,    -0.66212605039687e-4],
       [1,  2,  1,    -0.19232721156002],
       [1,  2,  5,    -0.25709043003438],
       [1,  3,  4,     0.16074868486251],
       [1,  4,  2,    -0.40092828925807e-1],
       [1,  4, 13,     0.39343422603254e-6],
       [1,  5,  9,    -0.75941377088144e-5],
       [1,  7,  3,     0.56250979351888e-3],
       [1,  9,  4,    -0.15608652257135e-4],
       [1, 10, 11,     0.11537996422951e-8],
       [1, 11,  4,     0.36582165144204e-6],
       [1, 13, 13,    -0.13251180074668e-11],
       [1, 15,  1,    -0.62639586912454e-9],
       [2,  1,  7,    -0.10793600908932],
       [2,  2,  1,     0.17611491008752e-1],
       [2,  2,  9,     0.22132295167546],
       [2,  2, 10,    -0.40247669763528],
       [2,  3, 10,     0.58083399985759],
       [2,  4,  3,     0.49969146990806e-2],
       [2,  4,  7,    -0.31358700712549e-1],
       [2,  4, 10,    -0.74315929710341],
       [2,  5, 10,     0.47807329915480],
       [2,  6,  6,     0.20527940895948e-1],
       [2,  6, 10,    -0.13636435110343],
       [2,  7, 10,     0.14180634400617e-1],
       [2,  9,  1,     0.83326504880713e-2],
       [2,  9,  2,    -0.29052336009585e-1],
       [2,  9,  3,     0.38615085574206e-1],
       [2,  9,  4,    -0.20393486513704e-1],
       [2,  9,  8,    -0.16554050063734e-2],
       [2, 10,  6,     0.19955571979541e-2],
       [2, 10,  9,     0.15870308324157e-3],
       [2, 12,  8,    -0.16388568342530e-4],
       [3,  3, 16,     0.43613615723811e-1],
       [3,  4, 22,     0.34994005463765e-1],
       [3,  4, 23,    -0.76788197844621e-1],
       [3,  5, 23,     0.22446277332006e-1],
       [4, 14, 10,    -0.62689710414685e-4],
       [6,  3, 50,    -0.55711118565645e-9],
       [6,  6, 44,    -0.19905718354408],
       [6,  6, 46,     0.31777497330738],
       [6,  6, 50,    -0.11841182425981],
    // [cᵢ, dᵢ, tᵢ,         nᵢ,          αᵢ, βᵢ, θᵢ, εᵢ]
       [x,  3,  0,  -0.31306260323435e2,   20,    150,   1.21,      1],
       [x,  3,  1,   0.31546140237781e2,   20,    150,   1.21,      1],
       [x,  3,  4,  -0.25213154341695e4,   20,    250,   1.25,      1],
    // [aᵢ,   bᵢ,    Bᵢ,         nᵢ,         Cᵢ, Dᵢ,  Aᵢ,   βᵢ]
       [7/2, 85/100, 0.2, -0.14874640856724, 28, 700, 0.32, 3/10],
       [7/2, 95/100, 0.2,  0.31806110878444, 32, 800, 0.32, 3/10]]

       
   /** Numerical values for IAPWS97 table 34.  These are used to solve
       saturation-pressure and saturation-temperature equations. */
   class var Table34 = [undef,  // No n0
        0.11670521452767e4, // n1
       -0.72421316703206e6, // n2
       -0.17073846940092e2, // n3
        0.12020824702470e5, // n4
       -0.32325550322333e7, // n5
        0.14915108613530e2, // n6
       -0.48232657361591e4, // n7
        0.40511340542057e6, // n8
       -0.23855557567849,   // n9
        0.65017534844798e3] // n10
}