Add chemical potential as abstract state method

For example:

```
auto fluids = strsplit("Methane&Ethane", '&');
auto x = { 0.6,0.4 };

shared_ptr<CoolProp::AbstractState> RP(CoolProp::AbstractState::factory("HEOS", fluids));
RP->set_mole_fractions(x);
RP->update(PT_INPUTS, 101325, 300);
auto mu1 = RP->chemical_potential(0);

RP.reset(CoolProp::AbstractState::factory("REFPROP", fluids));
RP->set_mole_fractions(x);
RP->update(PT_INPUTS, 101325, 300);
auto mu2 = RP->chemical_potential(0);
```

include/AbstractState.h

@@ -154,6 +154,8 @@ protected:
     virtual CoolPropDbl calc_fugacity_coefficient(std::size_t i){ throw NotImplementedError("calc_fugacity_coefficient is not implemented for this backend"); };
     /// Using this backend, calculate the fugacity in Pa
     virtual CoolPropDbl calc_fugacity(std::size_t i){ throw NotImplementedError("calc_fugacity is not implemented for this backend"); };
+    /// Using this backend, calculate the chemical potential in J/mol
+    virtual CoolPropDbl calc_chemical_potential(std::size_t i) { throw NotImplementedError("calc_chemical_potential 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"); };
 
@@ -646,6 +648,8 @@ public:
     double fugacity_coefficient(std::size_t i);
     /// Return the fugacity of the i-th component of the mixture
     double fugacity(std::size_t i);
+    /// Return the chemical potential of the i-th component of the mixture
+    double chemical_potential(std::size_t i);
     /// 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"

src/AbstractState.cpp

@@ -504,6 +504,10 @@ double AbstractState::fugacity(std::size_t i){
     // TODO: Cache the fug. coeff for each component
     return calc_fugacity(i);
 }
+double AbstractState::chemical_potential(std::size_t i) {
+    // TODO: Cache the chemical potential for each component
+    return calc_chemical_potential(i);
+}
 void AbstractState::build_phase_envelope(const std::string &type)
 {
     calc_phase_envelope(type);

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -2453,6 +2453,15 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity(std::size_t i)
     x_N_dependency_flag xN_flag = XN_DEPENDENT;
     return MixtureDerivatives::fugacity_i(*this, i, xN_flag);
 }
+CoolPropDbl HelmholtzEOSMixtureBackend::calc_chemical_potential(std::size_t i)
+{
+    x_N_dependency_flag xN_flag = XN_DEPENDENT;
+    double Tci = get_fluid_constant(i, iT_critical);
+    double rhoci = get_fluid_constant(i, irhomolar_critical);
+    double dnar_dni__const_T_V_nj = MixtureDerivatives::dnalphar_dni__constT_V_nj(*this, i, xN_flag);
+    double dna0_dni__const_T_V_nj = components[i].EOS().alpha0.base(tau()*(Tci / T_reducing()), delta()/(rhoci / rhomolar_reducing())) + 1 + log(mole_fractions[i]);
+    return gas_constant()*T()*(dna0_dni__const_T_V_nj + dnar_dni__const_T_V_nj);
+}
 CoolPropDbl HelmholtzEOSMixtureBackend::calc_phase_identification_parameter(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));

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -133,6 +133,7 @@ public:
         switch(param){
             case iP_critical: return fld.crit.p;
             case iT_critical: return fld.crit.T;
+            case irhomolar_critical: return fld.crit.rhomolar;
             case iacentric_factor: return fld.EOS().acentric;
             default:
                 throw ValueError(format("I don't know what to do with this fluid constant: %s", get_parameter_information(param,"short")));
@@ -256,7 +257,8 @@ public:
 
     CoolPropDbl calc_phase_identification_parameter(void);
     CoolPropDbl calc_fugacity(std::size_t i); 
-    CoolPropDbl calc_fugacity_coefficient(std::size_t i);
+    CoolPropDbl calc_fugacity_coefficient(std::size_t i); 
+    CoolPropDbl calc_chemical_potential(std::size_t i);
 
     /// Using this backend, calculate the flame hazard
     CoolPropDbl calc_flame_hazard(void){ return components[0].environment.FH;};

src/Backends/REFPROP/REFPROPMixtureBackend.cpp

@@ -928,6 +928,20 @@ CoolPropDbl REFPROPMixtureBackend::calc_fugacity(std::size_t i)
     //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
     return static_cast<CoolPropDbl>(f[i]*1000);
 }
+CoolPropDbl REFPROPMixtureBackend::calc_chemical_potential(std::size_t i)
+{
+    this->check_loaded_fluid();
+    double rho_mol_L = 0.001*_rhomolar;
+    long ierr = 0;
+    std::vector<double> chem_pot(mole_fractions.size());
+    char herr[255];
+    CHEMPOTdll(&_T, &rho_mol_L, &(mole_fractions[0]),  // Inputs
+        &(chem_pot[0]),                           // Outputs
+        &ierr, herr, errormessagelength);        // Error message
+    if (static_cast<int>(ierr) > 0) { throw ValueError(format("%s", herr).c_str()); }
+    //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
+    return static_cast<CoolPropDbl>(chem_pot[i]);
+}
 
 void REFPROPMixtureBackend::calc_phase_envelope(const std::string &type)
 {

src/Backends/REFPROP/REFPROPMixtureBackend.h

@@ -156,6 +156,7 @@ public:
 
     CoolPropDbl calc_fugacity_coefficient(std::size_t i);
     CoolPropDbl calc_fugacity(std::size_t i);
+    CoolPropDbl calc_chemical_potential(std::size_t i);
     CoolPropDbl calc_melting_line(int param, int given, CoolPropDbl value);
     bool has_melting_line(){return true;};
     double calc_melt_Tmax();