thermo/eos.py
def solve_T(self, P, V, quick=True):
        r'''Method to calculate `T` from a specified `P` and `V` for the PR
        EOS. Uses `Tc`, `a`, `b`, and `kappa` as well, obtained from the 
        class's namespace.

        Parameters
        ----------
        P : float
            Pressure, [Pa]
        V : float
            Molar volume, [m^3/mol]
        quick : bool, optional
            Whether to use a SymPy cse-derived expression (3x faster) or 
            individual formulas

        Returns
        -------
        T : float
            Temperature, [K]
        
        Notes
        -----
        The exact solution can be derived as follows, and is excluded for 
        breviety.
        
        >>> from sympy import *
        >>> P, T, V = symbols('P, T, V')
        >>> Tc, Pc, omega = symbols('Tc, Pc, omega')
        >>> R, a, b, kappa = symbols('R, a, b, kappa')
        
        >>> a_alpha = a*(1 + kappa*(1-sqrt(T/Tc)))**2
        >>> PR_formula = R*T/(V-b) - a_alpha/(V*(V+b)+b*(V-b)) - P
        >>> #solve(PR_formula, T)
        '''
        Tc, a, b, kappa = self.Tc, self.a, self.b, self.kappa
        if quick:
            x0 = V*V
            x1 = R*Tc
            x2 = x0*x1
            x3 = kappa*kappa
            x4 = a*x3
            x5 = b*x4
            x6 = 2.*V*b
            x7 = x1*x6
            x8 = b*b
            x9 = x1*x8
            x10 = V*x4
            x11 = (-x10 + x2 + x5 + x7 - x9)**2
            x12 = x0*x0
            x13 = R*R
            x14 = Tc*Tc
            x15 = x13*x14
            x16 = x8*x8
            x17 = a*a
            x18 = x3*x3
            x19 = x17*x18
            x20 = x0*V
            x21 = 2.*R*Tc*a*x3
            x22 = x8*b
            x23 = 4.*V*x22
            x24 = 4.*b*x20
            x25 = a*x1
            x26 = x25*x8
            x27 = x26*x3
            x28 = x0*x25
            x29 = x28*x3
            x30 = 2.*x8
            x31 = 6.*V*x27 - 2.*b*x29 + x0*x13*x14*x30 + x0*x19 + x12*x15 + x15*x16 - x15*x23 + x15*x24 - x19*x6 + x19*x8 - x20*x21 - x21*x22
            x32 = V - b
            x33 = 2.*(R*Tc*a*kappa)
            x34 = P*x2
            x35 = P*x5
            x36 = x25*x3
            x37 = P*x10
            x38 = P*R*Tc
            x39 = V*x17
            x40 = 2.*kappa*x3
            x41 = b*x17
            x42 = P*a*x3
            return -Tc*(2.*a*kappa*x11*sqrt(x32**3*(x0 + x6 - x8)*(P*x7 - P*x9 + x25 + x33 + x34 + x35 + x36 - x37))*(kappa + 1.) - x31*x32*((4.*V)*(R*Tc*a*b*kappa) + x0*x33 - x0*x35 + x12*x38 + x16*x38 + x18*x39 - x18*x41 - x20*x42 - x22*x42 - x23*x38 + x24*x38 + x25*x6 - x26 - x27 + x28 + x29 + x3*x39 - x3*x41 + x30*x34 - x33*x8 + x36*x6 + 3*x37*x8 + x39*x40 - x40*x41))/(x11*x31)
        else:
            return Tc*(-2*a*kappa*sqrt((V - b)**3*(V**2 + 2*V*b - b**2)*(P*R*Tc*V**2 + 2*P*R*Tc*V*b - P*R*Tc*b**2 - P*V*a*kappa**2 + P*a*b*kappa**2 + R*Tc*a*kappa**2 + 2*R*Tc*a*kappa + R*Tc*a))*(kappa + 1)*(R*Tc*V**2 + 2*R*Tc*V*b - R*Tc*b**2 - V*a*kappa**2 + a*b*kappa**2)**2 + (V - b)*(R**2*Tc**2*V**4 + 4*R**2*Tc**2*V**3*b + 2*R**2*Tc**2*V**2*b**2 - 4*R**2*Tc**2*V*b**3 + R**2*Tc**2*b**4 - 2*R*Tc*V**3*a*kappa**2 - 2*R*Tc*V**2*a*b*kappa**2 + 6*R*Tc*V*a*b**2*kappa**2 - 2*R*Tc*a*b**3*kappa**2 + V**2*a**2*kappa**4 - 2*V*a**2*b*kappa**4 + a**2*b**2*kappa**4)*(P*R*Tc*V**4 + 4*P*R*Tc*V**3*b + 2*P*R*Tc*V**2*b**2 - 4*P*R*Tc*V*b**3 + P*R*Tc*b**4 - P*V**3*a*kappa**2 - P*V**2*a*b*kappa**2 + 3*P*V*a*b**2*kappa**2 - P*a*b**3*kappa**2 + R*Tc*V**2*a*kappa**2 + 2*R*Tc*V**2*a*kappa + R*Tc*V**2*a + 2*R*Tc*V*a*b*kappa**2 + 4*R*Tc*V*a*b*kappa + 2*R*Tc*V*a*b - R*Tc*a*b**2*kappa**2 - 2*R*Tc*a*b**2*kappa - R*Tc*a*b**2 + V*a**2*kappa**4 + 2*V*a**2*kappa**3 + V*a**2*kappa**2 - a**2*b*kappa**4 - 2*a**2*b*kappa**3 - a**2*b*kappa**2))/((R*Tc*V**2 + 2*R*Tc*V*b - R*Tc*b**2 - V*a*kappa**2 + a*b*kappa**2)**2*(R**2*Tc**2*V**4 + 4*R**2*Tc**2*V**3*b + 2*R**2*Tc**2*V**2*b**2 - 4*R**2*Tc**2*V*b**3 + R**2*Tc**2*b**4 - 2*R*Tc*V**3*a*kappa**2 - 2*R*Tc*V**2*a*b*kappa**2 + 6*R*Tc*V*a*b**2*kappa**2 - 2*R*Tc*a*b**3*kappa**2 + V**2*a**2*kappa**4 - 2*V*a**2*b*kappa**4 + a**2*b**2*kappa**4))
def solve_T(self, P, V, quick=True):
        r'''Method to calculate `T` from a specified `P` and `V` for the PR
        EOS. Uses `Tc`, `a`, `b`, and `kappa` as well, obtained from the 
        class's namespace.

        Parameters
        ----------
        P : float
            Pressure, [Pa]
        V : float
            Molar volume, [m^3/mol]
        quick : bool, optional
            Whether to use a SymPy cse-derived expression (3x faster) or 
            individual formulas

        Returns
        -------
        T : float
            Temperature, [K]
        
        Notes
        -----
        The exact solution can be derived as follows, and is excluded for 
        breviety.
        
        >>> from sympy import *
        >>> P, T, V = symbols('P, T, V')
        >>> Tc, Pc, omega = symbols('Tc, Pc, omega')
        >>> R, a, b, kappa = symbols('R, a, b, kappa')
        
        >>> a_alpha = a*(1 + kappa*(1-sqrt(T/Tc)))**2
        >>> PR_formula = R*T/(V-b) - a_alpha/(V*(V+b)+b*(V-b)) - P
        >>> #solve(PR_formula, T)
        '''
        Tc, a, b, kappa = self.Tc, self.a, self.b, self.kappa
        if quick:
            x0 = V*V
            x1 = R*Tc
            x2 = x0*x1
            x3 = kappa*kappa
            x4 = a*x3
            x5 = b*x4
            x6 = 2.*V*b
            x7 = x1*x6
            x8 = b*b
            x9 = x1*x8
            x10 = V*x4
            x11 = (-x10 + x2 + x5 + x7 - x9)**2
            x12 = x0*x0
            x13 = R*R
            x14 = Tc*Tc
            x15 = x13*x14
            x16 = x8*x8
            x17 = a*a
            x18 = x3*x3
            x19 = x17*x18
            x20 = x0*V
            x21 = 2.*R*Tc*a*x3
            x22 = x8*b
            x23 = 4.*V*x22
            x24 = 4.*b*x20
            x25 = a*x1
            x26 = x25*x8
            x27 = x26*x3
            x28 = x0*x25
            x29 = x28*x3
            x30 = 2.*x8
            x31 = 6.*V*x27 - 2.*b*x29 + x0*x13*x14*x30 + x0*x19 + x12*x15 + x15*x16 - x15*x23 + x15*x24 - x19*x6 + x19*x8 - x20*x21 - x21*x22
            x32 = V - b
            x33 = 2.*(R*Tc*a*kappa)
            x34 = P*x2
            x35 = P*x5
            x36 = x25*x3
            x37 = P*x10
            x38 = P*R*Tc
            x39 = V*x17
            x40 = 2.*kappa*x3
            x41 = b*x17
            x42 = P*a*x3
            return -Tc*(2.*a*kappa*x11*sqrt(x32**3*(x0 + x6 - x8)*(P*x7 - P*x9 + x25 + x33 + x34 + x35 + x36 - x37))*(kappa + 1.) - x31*x32*((4.*V)*(R*Tc*a*b*kappa) + x0*x33 - x0*x35 + x12*x38 + x16*x38 + x18*x39 - x18*x41 - x20*x42 - x22*x42 - x23*x38 + x24*x38 + x25*x6 - x26 - x27 + x28 + x29 + x3*x39 - x3*x41 + x30*x34 - x33*x8 + x36*x6 + 3*x37*x8 + x39*x40 - x40*x41))/(x11*x31)
        else:
            return Tc*(-2*a*kappa*sqrt((V - b)**3*(V**2 + 2*V*b - b**2)*(P*R*Tc*V**2 + 2*P*R*Tc*V*b - P*R*Tc*b**2 - P*V*a*kappa**2 + P*a*b*kappa**2 + R*Tc*a*kappa**2 + 2*R*Tc*a*kappa + R*Tc*a))*(kappa + 1)*(R*Tc*V**2 + 2*R*Tc*V*b - R*Tc*b**2 - V*a*kappa**2 + a*b*kappa**2)**2 + (V - b)*(R**2*Tc**2*V**4 + 4*R**2*Tc**2*V**3*b + 2*R**2*Tc**2*V**2*b**2 - 4*R**2*Tc**2*V*b**3 + R**2*Tc**2*b**4 - 2*R*Tc*V**3*a*kappa**2 - 2*R*Tc*V**2*a*b*kappa**2 + 6*R*Tc*V*a*b**2*kappa**2 - 2*R*Tc*a*b**3*kappa**2 + V**2*a**2*kappa**4 - 2*V*a**2*b*kappa**4 + a**2*b**2*kappa**4)*(P*R*Tc*V**4 + 4*P*R*Tc*V**3*b + 2*P*R*Tc*V**2*b**2 - 4*P*R*Tc*V*b**3 + P*R*Tc*b**4 - P*V**3*a*kappa**2 - P*V**2*a*b*kappa**2 + 3*P*V*a*b**2*kappa**2 - P*a*b**3*kappa**2 + R*Tc*V**2*a*kappa**2 + 2*R*Tc*V**2*a*kappa + R*Tc*V**2*a + 2*R*Tc*V*a*b*kappa**2 + 4*R*Tc*V*a*b*kappa + 2*R*Tc*V*a*b - R*Tc*a*b**2*kappa**2 - 2*R*Tc*a*b**2*kappa - R*Tc*a*b**2 + V*a**2*kappa**4 + 2*V*a**2*kappa**3 + V*a**2*kappa**2 - a**2*b*kappa**4 - 2*a**2*b*kappa**3 - a**2*b*kappa**2))/((R*Tc*V**2 + 2*R*Tc*V*b - R*Tc*b**2 - V*a*kappa**2 + a*b*kappa**2)**2*(R**2*Tc**2*V**4 + 4*R**2*Tc**2*V**3*b + 2*R**2*Tc**2*V**2*b**2 - 4*R**2*Tc**2*V*b**3 + R**2*Tc**2*b**4 - 2*R*Tc*V**3*a*kappa**2 - 2*R*Tc*V**2*a*b*kappa**2 + 6*R*Tc*V*a*b**2*kappa**2 - 2*R*Tc*a*b**3*kappa**2 + V**2*a**2*kappa**4 - 2*V*a**2*b*kappa**4 + a**2*b**2*kappa**4))
def solve_T ( self , P , V , quick = True ) : Tc , a , b , kappa = self . Tc , self . a , self . b , self . kappa if quick : x0 = V * V x1 = R * Tc x2 = x0 * x1 x3 = kappa * kappa x4 = a * x3 x5 = b * x4 x6 = 2. * V * b x7 = x1 * x6 x8 = b * b x9 = x1 * x8 x10 = V * x4 x11 = ( - x10 + x2 + x5 + x7 - x9 ) ** 2 x12 = x0 * x0 x13 = R * R x14 = Tc * Tc x15 = x13 * x14 x16 = x8 * x8 x17 = a * a x18 = x3 * x3 x19 = x17 * x18 x20 = x0 * V x21 = 2. * R * Tc * a * x3 x22 = x8 * b x23 = 4. * V * x22 x24 = 4. * b * x20 x25 = a * x1 x26 = x25 * x8 x27 = x26 * x3 x28 = x0 * x25 x29 = x28 * x3 x30 = 2. * x8 x31 = 6. * V * x27 - 2. * b * x29 + x0 * x13 * x14 * x30 + x0 * x19 + x12 * x15 + x15 * x16 - x15 * x23 + x15 * x24 - x19 * x6 + x19 * x8 - x20 * x21 - x21 * x22 x32 = V - b x33 = 2. * ( R * Tc * a * kappa ) x34 = P * x2 x35 = P * x5 x36 = x25 * x3 x37 = P * x10 x38 = P * R * Tc x39 = V * x17 x40 = 2. * kappa * x3 x41 = b * x17 x42 = P * a * x3 return - Tc * ( 2. * a * kappa * x11 * sqrt ( x32 ** 3 * ( x0 + x6 - x8 ) * ( P * x7 - P * x9 + x25 + x33 + x34 + x35 + x36 - x37 ) ) * ( kappa + 1. ) - x31 * x32 * ( ( 4. * V ) * ( R * Tc * a * b * kappa ) + x0 * x33 - x0 * x35 + x12 * x38 + x16 * x38 + x18 * x39 - x18 * x41 - x20 * x42 - x22 * x42 - x23 * x38 + x24 * x38 + x25 * x6 - x26 - x27 + x28 + x29 + x3 * x39 - x3 * x41 + x30 * x34 - x33 * x8 + x36 * x6 + 3 * x37 * x8 + x39 * x40 - x40 * x41 ) ) / ( x11 * x31 ) else : return Tc * ( - 2 * a * kappa * sqrt ( ( V - b ) ** 3 * ( V ** 2 + 2 * V * b - b ** 2 ) * ( P * R * Tc * V ** 2 + 2 * P * R * Tc * V * b - P * R * Tc * b ** 2 - P * V * a * kappa ** 2 + P * a * b * kappa ** 2 + R * Tc * a * kappa ** 2 + 2 * R * Tc * a * kappa + R * Tc * a ) ) * ( kappa + 1 ) * ( R * Tc * V ** 2 + 2 * R * Tc * V * b - R * Tc * b ** 2 - V * a * kappa ** 2 + a * b * kappa ** 2 ) ** 2 + ( V - b ) * ( R ** 2 * Tc ** 2 * V ** 4 + 4 * R ** 2 * Tc ** 2 * V ** 3 * b + 2 * R ** 2 * Tc ** 2 * V ** 2 * b ** 2 - 4 * R ** 2 * Tc ** 2 * V * b ** 3 + R ** 2 * Tc ** 2 * b ** 4 - 2 * R * Tc * V ** 3 * a * kappa ** 2 - 2 * R * Tc * V ** 2 * a * b * kappa ** 2 + 6 * R * Tc * V * a * b ** 2 * kappa ** 2 - 2 * R * Tc * a * b ** 3 * kappa ** 2 + V ** 2 * a ** 2 * kappa ** 4 - 2 * V * a ** 2 * b * kappa ** 4 + a ** 2 * b ** 2 * kappa ** 4 ) * ( P * R * Tc * V ** 4 + 4 * P * R * Tc * V ** 3 * b + 2 * P * R * Tc * V ** 2 * b ** 2 - 4 * P * R * Tc * V * b ** 3 + P * R * Tc * b ** 4 - P * V ** 3 * a * kappa ** 2 - P * V ** 2 * a * b * kappa ** 2 + 3 * P * V * a * b ** 2 * kappa ** 2 - P * a * b ** 3 * kappa ** 2 + R * Tc * V ** 2 * a * kappa ** 2 + 2 * R * Tc * V ** 2 * a * kappa + R * Tc * V ** 2 * a + 2 * R * Tc * V * a * b * kappa ** 2 + 4 * R * Tc * V * a * b * kappa + 2 * R * Tc * V * a * b - R * Tc * a * b ** 2 * kappa ** 2 - 2 * R * Tc * a * b ** 2 * kappa - R * Tc * a * b ** 2 + V * a ** 2 * kappa ** 4 + 2 * V * a ** 2 * kappa ** 3 + V * a ** 2 * kappa ** 2 - a ** 2 * b * kappa ** 4 - 2 * a ** 2 * b * kappa ** 3 - a ** 2 * b * kappa ** 2 ) ) / ( ( R * Tc * V ** 2 + 2 * R * Tc * V * b - R * Tc * b ** 2 - V * a * kappa ** 2 + a * b * kappa ** 2 ) ** 2 * ( R ** 2 * Tc ** 2 * V ** 4 + 4 * R ** 2 * Tc ** 2 * V ** 3 * b + 2 * R ** 2 * Tc ** 2 * V ** 2 * b ** 2 - 4 * R ** 2 * Tc ** 2 * V * b ** 3 + R ** 2 * Tc ** 2 * b ** 4 - 2 * R * Tc * V ** 3 * a * kappa ** 2 - 2 * R * Tc * V ** 2 * a * b * kappa ** 2 + 6 * R * Tc * V * a * b ** 2 * kappa ** 2 - 2 * R * Tc * a * b ** 3 * kappa ** 2 + V ** 2 * a ** 2 * kappa ** 4 - 2 * V * a ** 2 * b * kappa ** 4 + a ** 2 * b ** 2 * kappa ** 4 ) )
r'''Method to calculate `T` from a specified `P` and `V` for the PR
        EOS. Uses `Tc`, `a`, `b`, and `kappa` as well, obtained from the 
        class's namespace.

        Parameters
        ----------
        P : float
            Pressure, [Pa]
        V : float
            Molar volume, [m^3/mol]
        quick : bool, optional
            Whether to use a SymPy cse-derived expression (3x faster) or 
            individual formulas

        Returns
        -------
        T : float
            Temperature, [K]
        
        Notes
        -----
        The exact solution can be derived as follows, and is excluded for 
        breviety.
        
        >>> from sympy import *
        >>> P, T, V = symbols('P, T, V')
        >>> Tc, Pc, omega = symbols('Tc, Pc, omega')
        >>> R, a, b, kappa = symbols('R, a, b, kappa')
        
        >>> a_alpha = a*(1 + kappa*(1-sqrt(T/Tc)))**2
        >>> PR_formula = R*T/(V-b) - a_alpha/(V*(V+b)+b*(V-b)) - P
        >>> #solve(PR_formula, T)
r Method to calculate T from a specified P and V for the PR EOS . Uses Tc a b and kappa as well obtained from the class s namespace .