Make the examples in the documentation easier to copy with copybutton…

….js, after a tweak to get it to work with recent versions of sphinx

docs/_static/copybutton.js

@@ -0,0 +1,63 @@
+$(document).ready(function() {
+    /* Add a [>>>] button on the top-right corner of code samples to hide
+     * the >>> and ... prompts and the output and thus make the code
+     * copyable. */
+    var div = $('.highlight-python .highlight,' +
+                '.highlight-default .highlight,' +
+                '.highlight-python3 .highlight')
+    var pre = div.find('pre');
+
+    // get the styles from the current theme
+    pre.parent().parent().css('position', 'relative');
+    var hide_text = 'Hide the prompts and output';
+    var show_text = 'Show the prompts and output';
+    var border_width = pre.css('border-top-width');
+    var border_style = pre.css('border-top-style');
+    var border_color = pre.css('border-top-color');
+    var button_styles = {
+        'cursor':'pointer', 'position': 'absolute', 'top': '0', 'right': '0',
+        'border-color': border_color, 'border-style': border_style,
+        'border-width': border_width, 'color': border_color, 'text-size': '75%',
+        'font-family': 'monospace', 'padding-left': '0.2em', 'padding-right': '0.2em',
+        'border-radius': '0 3px 0 0'
+    }
+
+    // create and add the button to all the code blocks that contain >>>
+    div.each(function(index) {
+        var jthis = $(this);
+        if (jthis.find('.gp').length > 0) {
+            var button = $('<span class="copybutton">&gt;&gt;&gt;</span>');
+            button.css(button_styles)
+            button.attr('title', hide_text);
+            button.data('hidden', 'false');
+            jthis.prepend(button);
+        }
+        // tracebacks (.gt) contain bare text elements that need to be
+        // wrapped in a span to work with .nextUntil() (see later)
+        jthis.find('pre:has(.gt)').contents().filter(function() {
+            return ((this.nodeType == 3) && (this.data.trim().length > 0));
+        }).wrap('<span>');
+    });
+
+    // define the behavior of the button when it's clicked
+    $('.copybutton').click(function(e){
+        e.preventDefault();
+        var button = $(this);
+        if (button.data('hidden') === 'false') {
+            // hide the code output
+            button.parent().find('.go, .gp, .gt').hide();
+            button.next('pre').find('.gt').nextUntil('.gp, .go').css('visibility', 'hidden');
+            button.css('text-decoration', 'line-through');
+            button.attr('title', show_text);
+            button.data('hidden', 'true');
+        } else {
+            // show the code output
+            button.parent().find('.go, .gp, .gt').show();
+            button.next('pre').find('.gt').nextUntil('.gp, .go').css('visibility', 'visible');
+            button.css('text-decoration', 'none');
+            button.attr('title', hide_text);
+            button.data('hidden', 'false');
+        }
+    });
+});
+

docs/_templates/layout.html

@@ -0,0 +1,5 @@
+{% extends "!layout.html" %}
+{% block extrahead %}
+    {% if not embedded %}<script type="text/javascript" src="{{ pathto('_static/copybutton.js', 1) }}"></script>{% endif %}
+{{ super() }}
+{% endblock %}

thermo/eos.py

@@ -132,18 +132,25 @@ def set_from_PT(self, Vs):
             Three possible molar volumes, [m^3/mol]
         '''
         # All roots will have some imaginary component; ignore them if > 1E-9
-        imaginary_roots_count = len([True for i in Vs if abs(i.imag) > 1E-9]) 
-        if imaginary_roots_count == 2: 
-            V = [i for i in Vs if abs(i.imag) < 1E-9][0].real
+        good_roots = []
+        bad_roots = []
+        for i in Vs:
+            j = i.real
+            if abs(i.imag) > 1E-9 or j < 0:
+                bad_roots.append(i)
+            else:
+                good_roots.append(j)
+
+        if len(bad_roots) == 2: 
+            V = good_roots[0]
             self.phase = self.set_properties_from_solution(self.T, self.P, V, self.b, self.delta, self.epsilon, self.a_alpha, self.da_alpha_dT, self.d2a_alpha_dT2)
             if self.phase == 'l':
                 self.V_l = V
             else:
                 self.V_g = V
         else:
             # Even in the case of three real roots, it is still the min/max that make sense
-            Vs = [i.real for i in Vs]
-            self.V_l, self.V_g = min(Vs), max(Vs)
+            self.V_l, self.V_g = min(good_roots), max(good_roots)
             [self.set_properties_from_solution(self.T, self.P, V, self.b, self.delta, self.epsilon, self.a_alpha, self.da_alpha_dT, self.d2a_alpha_dT2) for V in [self.V_l, self.V_g]]
             self.phase = 'l/g'
 
@@ -372,6 +379,7 @@ def set_properties_from_solution(self, T, P, V, b, delta, epsilon, a_alpha,
         phase = 'l' if PIP > 1 else 'g' # phase_identification_parameter_phase(PIP)
 
         if phase == 'l':
+            self.Z_l = self.P*V/(R*self.T)
             self.beta_l, self.kappa_l = beta, kappa
             self.PIP_l, self.Cp_minus_Cv_l = PIP, Cp_m_Cv
 
@@ -389,6 +397,7 @@ def set_properties_from_solution(self, T, P, V, b, delta, epsilon, a_alpha,
             self.fugacity_l, self.phi_l = fugacity, phi
             self.Cp_dep_l, self.Cv_dep_l = Cp_dep, Cv_dep
         else:
+            self.Z_g = self.P*V/(R*self.T)
             self.beta_g, self.kappa_g = beta, kappa
             self.PIP_g, self.Cp_minus_Cv_g = PIP, Cp_m_Cv
 

thermo/eos_mix.py

@@ -24,14 +24,18 @@
 __all__ = ['GCEOSMIX', 'PRMIX', 'SRKMIX', 'PR78MIX', 'VDWMIX', 'PRSVMIX', 
 'PRSV2MIX', 'TWUPRMIX', 'TWUSRKMIX', 'APISRKMIX']
 from scipy.optimize import newton
+from scipy.misc import derivative
 from thermo.utils import Cp_minus_Cv, isobaric_expansion, isothermal_compressibility, phase_identification_parameter
 from thermo.utils import R
 from thermo.utils import log, exp, sqrt
 from thermo.eos import *
+import sys
 
 R2 = R*R
 two_root_two = 2*2**0.5
 root_two = sqrt(2.)
+log_min = log(sys.float_info.min)
+
 
 class GCEOSMIX(GCEOS):
     r'''Class for solving a generic pressure-explicit three-parameter cubic 
@@ -224,6 +228,65 @@ def fugacities(self, xs=None, ys=None):
             self.fugacities_g = [phi*y*self.P for phi, y in zip(self.phis_g, ys)]
             self.lnphis_g = [log(i) for i in self.phis_g]
 
+
+    def _dphi_dn(self, zi, i, phase):
+        z_copy = list(self.zs)
+        z_copy.pop(i)
+        z_sum = sum(z_copy) + zi
+        z_copy = [j/z_sum if j else 0 for j in z_copy]
+        z_copy.insert(i, zi)
+
+        eos = self.to_TP_zs(self.T, self.P, z_copy)
+        if phase == 'g':
+            return eos.phis_g[i]
+        elif phase == 'l':
+            return eos.phis_l[i]
+
+    def _dfugacity_dn(self, zi, i, phase):
+        z_copy = list(self.zs)
+        z_copy.pop(i)
+        z_sum = sum(z_copy) + zi
+        z_copy = [j/z_sum if j else 0 for j in z_copy]
+        z_copy.insert(i, zi)
+
+        eos = self.to_TP_zs(self.T, self.P, z_copy)
+        if phase == 'g':
+            return eos.fugacities_g[i]
+        elif phase == 'l':
+            return eos.fugacities_l[i]
+
+
+    def fugacities_partial_derivatives(self, xs=None, ys=None):
+        if self.phase in ['l', 'l/g']:
+            if xs is None:
+                xs = self.zs
+            self.dphis_dni_l = [derivative(self._dphi_dn, xs[i], args=[i, 'l'], dx=1E-7, n=1) for i in self.cmps]
+            self.dfugacities_dni_l = [derivative(self._dfugacity_dn, xs[i], args=[i, 'l'], dx=1E-7, n=1) for i in self.cmps]
+            self.dlnphis_dni_l = [dphi/phi for dphi, phi in zip(self.dphis_dni_l, self.phis_l)]
+        if self.phase in ['g', 'l/g']:
+            if ys is None:
+                ys = self.zs
+            self.dphis_dni_g = [derivative(self._dphi_dn, ys[i], args=[i, 'g'], dx=1E-7, n=1) for i in self.cmps]
+            self.dfugacities_dni_g = [derivative(self._dfugacity_dn, ys[i], args=[i, 'g'], dx=1E-7, n=1) for i in self.cmps]
+            self.dlnphis_dni_g = [dphi/phi for dphi, phi in zip(self.dphis_dni_g, self.phis_g)]
+            # confirmed the relationship of the above 
+            # There should be an easy way to get dfugacities_dn_g but I haven't figured it out
+
+    def fugacities_partial_derivatives_2(self, xs=None, ys=None):
+        if self.phase in ['l', 'l/g']:
+            if xs is None:
+                xs = self.zs
+            self.d2phis_dni2_l = [derivative(self._dphi_dn, xs[i], args=[i, 'l'], dx=1E-5, n=2) for i in self.cmps]
+            self.d2fugacities_dni2_l = [derivative(self._dfugacity_dn, xs[i], args=[i, 'l'], dx=1E-5, n=2) for i in self.cmps]
+            self.d2lnphis_dni2_l = [d2phi/phi  - dphi*dphi/(phi*phi) for d2phi, dphi, phi in zip(self.d2phis_dni2_l, self.dphis_dni_l, self.phis_l)]
+        if self.phase in ['g', 'l/g']:
+            if ys is None:
+                ys = self.zs
+            self.d2phis_dni2_g = [derivative(self._dphi_dn, ys[i], args=[i, 'g'], dx=1E-5, n=2) for i in self.cmps]
+            self.d2fugacities_dni2_g = [derivative(self._dfugacity_dn, ys[i], args=[i, 'g'], dx=1E-5, n=2) for i in self.cmps]
+            self.d2lnphis_dni2_g = [d2phi/phi  - dphi*dphi/(phi*phi) for d2phi, dphi, phi in zip(self.d2phis_dni2_g, self.dphis_dni_g, self.phis_g)]
+        # second derivative lns confirmed
+
     def solve_T(self, P, V, quick=True):
         r'''Generic method to calculate `T` from a specified `P` and `V`.
         Provides SciPy's `newton` solver, and iterates to solve the general
@@ -401,13 +464,17 @@ def fugacity_coefficients(self, Z, zs):
         .. [2] Walas, Stanley M. Phase Equilibria in Chemical Engineering. 
            Butterworth-Heinemann, 1985.
         '''
+        from cmath import log
         A = self.a_alpha*self.P/(R2*self.T*self.T)
         B = self.b*self.P/(R*self.T)
         phis = []
         for i in self.cmps:
-            t1 = self.bs[i]/self.b*(Z - 1.) - log(Z - B)
+            # The two log terms need to use a complex log; typically these are
+            # calculated at "liquid" volume solutions which are unstable
+            # and cannot exist
+            t1 = self.bs[i]/self.b*(Z - 1.) - log(Z - B).real
             t2 = 2./self.a_alpha*sum([zs[j]*self.a_alpha_ijs[i][j] for j in self.cmps])
-            t3 = t1 - A/(two_root_two*B)*(t2 - self.bs[i]/self.b)*log((Z + (root_two + 1.)*B)/(Z - (root_two - 1.)*B))
+            t3 = t1 - A/(two_root_two*B)*(t2 - self.bs[i]/self.b)*log((Z + (root_two + 1.)*B)/(Z - (root_two - 1.)*B)).real
             phis.append(exp(t3))
         return phis
 

thermo/volume.py

@@ -1314,7 +1314,7 @@ def Amgat(xs, Vms):
 
     Parameters
     ----------
-    xs: array
+    xs : array
         Mole fractions of each component, []
     Vms : array
         Molar volumes of each fluids at conditions [m^3/mol]