Implement the PIP parameter of Venkatharathnam for HEOS backend through to python and B and C virial coefficients

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>

include/AbstractState.h

@@ -145,7 +145,8 @@ protected:
     virtual CoolPropDbl calc_gas_constant(void){throw NotImplementedError("calc_gas_constant is not implemented for this backend");};
     /// Using this backend, calculate the fugacity coefficient (dimensionless)
     virtual CoolPropDbl calc_fugacity_coefficient(int i){throw NotImplementedError("calc_fugacity_coefficient is not implemented for this backend");};
-
+    /// Using this backend, calculate the phase identification parameter (PIP)
+    virtual CoolPropDbl calc_PIP(void){throw NotImplementedError("calc_PIP is not implemented for this backend");};
 
     // Derivatives of residual helmholtz energy
     /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r\f$ (dimensionless)
@@ -326,6 +327,7 @@ protected:
     /// Using this backend, get the molar density in mol/m^3
     virtual long double calc_rhomolar(void){return _rhomolar;}
 
+
 public:
 
     AbstractState():_fluid_type(FLUID_TYPE_UNDEFINED),_phase(iphase_unknown),_rhospline(-_HUGE),_dsplinedp(-_HUGE),_dsplinedh(-_HUGE){clear();}
@@ -491,9 +493,13 @@ public:
     double acentric_factor(void);
     /// Return the mole-fraction weighted gas constant in J/mol/K
     double gas_constant(void);
+    /// Return the B virial coefficient
     double Bvirial(void);
+    /// Return the derivative of the B virial coefficient with respect to temperature
     double dBvirial_dT(void);
+    /// Return the C virial coefficient
     double Cvirial(void);
+    /// Return the derivative of the C virial coefficient with respect to temperature
     double dCvirial_dT(void);
     /// Return the compressibility factor \f$ Z = p/(rho R T) \f$
     double compressibility_factor(void);
@@ -530,7 +536,10 @@ public:
     /// Return the isobaric expansion coefficient \f$ \beta = \frac{1}{v}\left.\frac{\partial v}{\partial T}\right|_p = -\frac{1}{\rho}\left.\frac{\partial \rho}{\partial T}\right|_p\f$  in 1/K
     double isobaric_expansion_coefficient(void);
     double fugacity_coefficient(int i);
-    //double fundamental_derivative_of_gas_dynamics(void);
+    /// Return the fundamental derivative of gas dynamics
+    //double fundamental_derivative_of_gas_dynamics(void){return this->second_partial_deriv(iP, iDmolar, iSmolar, iDmolar, iSmolar)/pow(speed_sound(), 2)/2/pow(this->rhomolar(),3);};
+    /// Return the phase identification parameter (PIP) of G. Venkatarathnam and L.R. Oellrich, "Identification of the phase of a fluid using partial derivatives of pressure, volume, and temperature without reference to saturation properties: Applications in phase equilibria calculations"
+    double PIP(){return calc_PIP();};
     
     std::vector<CoolPropDbl> mole_fractions_liquid(void){return calc_mole_fractions_liquid();};
     std::vector<CoolPropDbl> mole_fractions_vapor(void){return calc_mole_fractions_vapor();};

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -262,6 +262,9 @@ public:
     CoolPropDbl calc_rhomolar_reducing(void){return get_reducing_state().rhomolar;};
     CoolPropDbl calc_p_reducing(void){return get_reducing_state().p;};
 
+    // Calculate the phase identification parameter of Venkatarathnam et al, Fluid Phase Equilibria
+    CoolPropDbl calc_PIP(void){return 2-rhomolar()*(second_partial_deriv(iP,iDmolar, iT, iT, iDmolar)/first_partial_deriv(iP, iT, iDmolar)-second_partial_deriv(iP, iDmolar, iT, iDmolar, iT)/first_partial_deriv(iP, iDmolar, iT));};
+
     std::string calc_name(void);
 	std::vector<std::string> calc_fluid_names(void);
 

wrappers/Python/CoolProp/AbstractState.pxd

@@ -67,6 +67,7 @@ cdef class AbstractState:
     cpdef double Prandtl(self) except *
     cpdef double Bvirial(self) except *
     cpdef double Cvirial(self) except *
+    cpdef double PIP(self) except *
     
     cpdef double molar_mass(self) except *
     cpdef double acentric_factor(self) except*

wrappers/Python/CoolProp/AbstractState.pyx

@@ -178,6 +178,9 @@ cdef class AbstractState:
     cpdef double Cvirial(self) except *: 
         """ Get the C virial coefficient - wrapper of c++ function :cpapi:`CoolProp::AbstractState::Cvirial(void)` """
         return self.thisptr.Cvirial()
+    cpdef double PIP(self) except *: 
+        """ Get the phase identification parameter - wrapper of c++ function :cpapi:`CoolProp::AbstractState::PIP(void)` """
+        return self.thisptr.PIP()
     
     cpdef mole_fractions_liquid(self):
         """ Get the mole fractions of the liquid phase - wrapper of c++ function :cpapi:`CoolProp::AbstractState::mole_fractions_liquid(void)` """

wrappers/Python/CoolProp/cAbstractState.pxd

@@ -82,6 +82,7 @@ cdef extern from "AbstractState.h" namespace "CoolProp":
         double Prandtl() except +ValueError
         double Bvirial() except +ValueError
         double Cvirial() except +ValueError
+        double PIP() except +ValueError
         
         double keyed_output(constants_header.parameters) except+ValueError
         double trivial_keyed_output(constants_header.parameters) except+ValueError