def _phir(self, tau, delta):
        """Residual contribution to the free Helmholtz energy

        Parameters
        ----------
        tau : float
            Inverse reduced temperature Tc/T, [-]
        delta : float
            Reduced density rho/rhoc, [-]

        Returns
        -------
        prop : dict
          Dictionary with residual adimensional helmholtz energy and deriv:
            * fir
            * firt: ∂fir/∂τ|δ,x
            * fird: ∂fir/∂δ|τ,x
            * firtt: ∂²fir/∂τ²|δ,x
            * firdt: ∂²fir/∂τ∂δ|x
            * firdd: ∂²fir/∂δ²|τ,x

        References
        ----------
        IAPWS, Revised Release on the IAPWS Formulation 1995 for the
        Thermodynamic Properties of Ordinary Water Substance for General and
        Scientific Use, September 2016, Table 5
        http://www.iapws.org/relguide/IAPWS-95.html
        """
        fir = fird = firdd = firt = firtt = firdt = 0

        # Polinomial terms
        nr1 = self._constants.get("nr1", [])
        d1 = self._constants.get("d1", [])
        t1 = self._constants.get("t1", [])
        for n, d, t in zip(nr1, d1, t1):
            fir += n*delta**d*tau**t
            fird += n*d*delta**(d-1)*tau**t
            firdd += n*d*(d-1)*delta**(d-2)*tau**t
            firt += n*t*delta**d*tau**(t-1)
            firtt += n*t*(t-1)*delta**d*tau**(t-2)
            firdt += n*t*d*delta**(d-1)*tau**(t-1)

        # Exponential terms
        nr2 = self._constants.get("nr2", [])
        d2 = self._constants.get("d2", [])
        g2 = self._constants.get("gamma2", [])
        t2 = self._constants.get("t2", [])
        c2 = self._constants.get("c2", [])
        for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
            fir += n*delta**d*tau**t*exp(-g*delta**c)
            fird += n*exp(-g*delta**c)*delta**(d-1)*tau**t*(d-g*c*delta**c)
            firdd += n*exp(-g*delta**c)*delta**(d-2)*tau**t * \
                ((d-g*c*delta**c)*(d-1-g*c*delta**c)-g**2*c**2*delta**c)
            firt += n*t*delta**d*tau**(t-1)*exp(-g*delta**c)
            firtt += n*t*(t-1)*delta**d*tau**(t-2)*exp(-g*delta**c)
            firdt += n*t*delta**(d-1)*tau**(t-1)*(d-g*c*delta**c)*exp(
                -g*delta**c)

        # Gaussian terms
        nr3 = self._constants.get("nr3", [])
        d3 = self._constants.get("d3", [])
        t3 = self._constants.get("t3", [])
        a3 = self._constants.get("alfa3", [])
        e3 = self._constants.get("epsilon3", [])
        b3 = self._constants.get("beta3", [])
        g3 = self._constants.get("gamma3", [])
        for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
            fir += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)
            fird += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    d/delta-2*a*(delta-e))
            firdd += n*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    -2*a*delta**d + 4*a**2*delta**d*(delta-e)**2 -
                    4*d*a*delta**(d-1)*(delta-e) + d*(d-1)*delta**(d-2))
            firt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    t/tau-2*b*(tau-g))
            firtt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    (t/tau-2*b*(tau-g))**2-t/tau**2-2*b)
            firdt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    t/tau-2*b*(tau-g))*(d/delta-2*a*(delta-e))

        # Non analitic terms
        nr4 = self._constants.get("nr4", [])
        a4 = self._constants.get("a4", [])
        b4 = self._constants.get("b4", [])
        Ai = self._constants.get("A", [])
        Bi = self._constants.get("B", [])
        Ci = self._constants.get("C", [])
        Di = self._constants.get("D", [])
        bt4 = self._constants.get("beta4", [])
        for n, a, b, A, B, C, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
            Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
            F = exp(-C*(delta-1)**2-D*(tau-1)**2)
            Fd = -2*C*F*(delta-1)
            Fdd = 2*C*F*(2*C*(delta-1)**2-1)
            Ft = -2*D*F*(tau-1)
            Ftt = 2*D*F*(2*D*(tau-1)**2-1)
            Fdt = 4*C*D*F*(delta-1)*(tau-1)

            Delta = Tita**2+B*((delta-1)**2)**a
            Deltad = (delta-1)*(A*Tita*2/bt*((delta-1)**2)**(0.5/bt-1) +
                                2*B*a*((delta-1)**2)**(a-1))
            if delta == 1:
                Deltadd = 0
            else:
                Deltadd = Deltad/(delta-1)+(delta-1)**2*(
                    4*B*a*(a-1)*((delta-1)**2)**(a-2) +
                    2*A**2/bt**2*(((delta-1)**2)**(0.5/bt-1))**2 +
                    A*Tita*4/bt*(0.5/bt-1)*((delta-1)**2)**(0.5/bt-2))

            DeltaBd = b*Delta**(b-1)*Deltad
            DeltaBdd = b*(Delta**(b-1)*Deltadd+(b-1)*Delta**(b-2)*Deltad**2)
            DeltaBt = -2*Tita*b*Delta**(b-1)
            DeltaBtt = 2*b*Delta**(b-1)+4*Tita**2*b*(b-1)*Delta**(b-2)
            DeltaBdt = -A*b*2/bt*Delta**(b-1)*(delta-1)*((delta-1)**2)**(
                0.5/bt-1)-2*Tita*b*(b-1)*Delta**(b-2)*Deltad

            fir += n*Delta**b*delta*F
            fird += n*(Delta**b*(F+delta*Fd)+DeltaBd*delta*F)
            firdd += n*(Delta**b*(2*Fd+delta*Fdd) + 2*DeltaBd*(F+delta*Fd) +
                        DeltaBdd*delta*F)
            firt += n*delta*(DeltaBt*F+Delta**b*Ft)
            firtt += n*delta*(DeltaBtt*F+2*DeltaBt*Ft+Delta**b*Ftt)
            firdt += n*(Delta**b*(Ft+delta*Fdt)+delta*DeltaBd*Ft +
                        DeltaBt*(F+delta*Fd)+DeltaBdt*delta*F)

        prop = {}
        prop["fir"] = fir
        prop["firt"] = firt
        prop["firtt"] = firtt
        prop["fird"] = fird
        prop["firdd"] = firdd
        prop["firdt"] = firdt
        return prop
Residual contribution to the free Helmholtz energy

        Parameters
        ----------
        tau : float
            Inverse reduced temperature Tc/T, [-]
        delta : float
            Reduced density rho/rhoc, [-]

        Returns
        -------
        prop : dict
          Dictionary with residual adimensional helmholtz energy and deriv:
            * fir
            * firt: ∂fir/∂τ|δ,x
            * fird: ∂fir/∂δ|τ,x
            * firtt: ∂²fir/∂τ²|δ,x
            * firdt: ∂²fir/∂τ∂δ|x
            * firdd: ∂²fir/∂δ²|τ,x

        References
        ----------
        IAPWS, Revised Release on the IAPWS Formulation 1995 for the
        Thermodynamic Properties of Ordinary Water Substance for General and
        Scientific Use, September 2016, Table 5
        http://www.iapws.org/relguide/IAPWS-95.html