def fugacity_coefficients(self, Z, zs):
from cmath import log<EOL>A = self.a_alpha*self.P/(R2*self.T*self.T)<EOL>B = self.b*self.P/(R*self.T)<EOL>phis = []<EOL>for i in self.cmps:<EOL><INDENT>t1 = self.bs[i]/self.b*(Z - <NUM_LIT:1.>) - log(Z - B).real<EOL>t2 = <NUM_LIT>/self.a_alpha*sum([zs[j]*self.a_alpha_ijs[i][j] for j in self.cmps])<EOL>t3 = t1 - A/(two_root_two*B)*(t2 - self.bs[i]/self.b)*log((Z + (root_two + <NUM_LIT:1.>)*B)/(Z - (root_two - <NUM_LIT:1.>)*B)).real<EOL>phis.append(exp(t3))<EOL><DEDENT>return phis<EOL>
r'''Literature formula for calculating fugacity coefficients for each
        species in a mixture. Verified numerically. Applicable to most 
        derivatives of the Peng-Robinson equation of state as well.
        Called by `fugacities` on initialization, or by a solver routine 
        which is performing a flash calculation.

        .. math::
            \ln \hat \phi_i = \frac{B_i}{B}(Z-1)-\ln(Z-B) + \frac{A}{2\sqrt{2}B}
            \left[\frac{B_i}{B} - \frac{2}{a\alpha}\sum_i y_i(a\alpha)_{ij}\right]
            \log\left[\frac{Z + (1+\sqrt{2})B}{Z-(\sqrt{2}-1)B}\right]

            A = \frac{(a\alpha)P}{R^2 T^2}

            B = \frac{b P}{RT}

        Parameters
        ----------
        Z : float
            Compressibility of the mixture for a desired phase, [-]
        zs : list[float], optional
            List of mole factions, either overall or in a specific phase, [-]

        Returns
        -------
        phis : float
            Fugacity coefficient for each species, [-]

        References
        ----------
        .. [1] Peng, Ding-Yu, and Donald B. Robinson. "A New Two-Constant  
           Equation of State." Industrial & Engineering Chemistry Fundamentals 
           15, no. 1 (February 1, 1976): 59-64. doi:10.1021/i160057a011.
        .. [2] Walas, Stanley M. Phase Equilibria in Chemical Engineering. 
           Butterworth-Heinemann, 1985.