CubicBackend.h

Go to the documentation of this file.

/*
 
The goal of this backend is to allow the Helmholtz-based structure for cubics and to replace the entire
multi-fluid model with a one-fluid model.  The other changes are relatively trivial.  The primary
change is to replace the core residual Helmholtz energy derivatives from HelmholtzEOSMixtureBackend
with the derivatives from this class.
 
The core code for the Helmholtz translations is from the publication
"Helmholtz energy translations for common cubic equations of state for use in one-fluid and multi-fluid mixture models"
by Ian H. Bell and Andreas Jaeger, J. Res. NIST, 2016
 
*/
 
#ifndef CUBICBACKEND_H_
#define CUBICBACKEND_H_
 
#include <array>
 
#include "CoolProp/detail/tools.h"
#include "CoolProp/DataStructures.h"
#include "GeneralizedCubic.h"
#include "CubicsLibrary.h"
#include "CoolProp/Configuration.h"
#include "CoolProp/AbstractState.h"
#include "Backends/Helmholtz/HelmholtzEOSMixtureBackend.h"
#include "CoolProp/Exceptions.h"
#include "CoolProp/superancillary/cubicsuperancillary.h"
#include <algorithm>
#include <vector>
 
namespace CoolProp {
 
// Forward declaration for use in initialization of AbstractCubicBackend
class CubicResidualHelmholtz;
 
class AbstractCubicBackend : public HelmholtzEOSMixtureBackend
{
   protected:
    shared_ptr<AbstractCubic> cubic;
    std::vector<CubicLibrary::CubicsValues> components;  
   public:
    void setup(bool generate_SatL_and_SatV = true);
 
    void set_alpha_from_components();
 
    void set_alpha0_from_components();
 
    shared_ptr<AbstractCubic>& get_cubic() {
        return cubic;
    };
 
    std::vector<std::string> calc_fluid_names() override;
 
    bool using_mole_fractions() override {
        return true;
    };
    bool using_mass_fractions() override {
        return false;
    };
    bool using_volu_fractions() override {
        return false;
    };
 
    void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) override {
        throw NotImplementedError("Mass composition has not been implemented.");
    };
    void set_volu_fractions(const std::vector<CoolPropDbl>& volu_fractions) override {
        throw NotImplementedError("Volume composition has not been implemented.");
    };
    const std::vector<CoolPropDbl>& get_mole_fractions() override {
        return this->mole_fractions;
    };
 
    const double get_fluid_constant(std::size_t i, parameters param) const override {
        switch (param) {
            case iP_critical:
                return cubic->get_pc()[i];
            case iT_reducing:
            case iT_critical:
                return cubic->get_Tc()[i];
            case iacentric_factor:
                return cubic->get_acentric()[i];
            case imolar_mass:
                return components[i].molemass;
            case iT_triple:
                return HelmholtzEOSMixtureBackend::get_components()[i].EOS().sat_min_liquid.T;  // From the base class data structure
            case iP_triple:
                return HelmholtzEOSMixtureBackend::get_components()[i].EOS().sat_min_liquid.p;  // From the base class data structure
            case irhomolar_reducing:
            case irhomolar_critical:
                return components[i].rhomolarc;
            case igas_constant:
                return get_config_double(R_U_CODATA);
            default:
                throw ValueError(format("I don't know what to do with this fluid constant: %s", get_parameter_information(param, "short").c_str()));
        }
    }
 
    std::string fluid_param_string(const std::string&) override;
    CoolPropDbl calc_gas_constant() override {
        return cubic->get_R_u();
    };
    SimpleState calc_reducing_state_nocache(const std::vector<CoolPropDbl>& mole_fractions) override {
        SimpleState reducing;
        reducing.T = cubic->get_Tr();
        reducing.rhomolar = cubic->get_rhor();
        return reducing;
    };
    CoolPropDbl calc_rhomolar_max_bound() override {
        return 0.9 / get_cubic()->bm_term(mole_fractions);
    };
    CoolPropDbl calc_reduced_density() override {
        return _rhomolar / get_cubic()->get_rhor();
    };
    CoolPropDbl calc_reciprocal_reduced_temperature() override {
        return get_cubic()->get_Tr() / _T;
    };
    std::vector<double> spinodal_densities();
 
    CoolPropDbl calc_T_critical() override {
        if (is_pure_or_pseudopure) {
            return cubic->get_Tc()[0];
        } else {
            return HelmholtzEOSMixtureBackend::calc_T_critical();
        }
    };
    CoolPropDbl calc_p_critical() override {
        if (is_pure_or_pseudopure) {
            return cubic->get_pc()[0];
        } else {
            return HelmholtzEOSMixtureBackend::calc_p_critical();
        }
    };
    CoolPropDbl calc_acentric_factor() override {
        if (is_pure_or_pseudopure) {
            return cubic->get_acentric()[0];
        } else {
            throw ValueError("acentric factor cannot be calculated for mixtures");
        }
    }
    CoolPropDbl calc_rhomolar_critical() override {
        if (is_pure_or_pseudopure) {
            // Curve fit from all the pure fluids in CoolProp (thanks to recommendation of A. Kazakov)
            double v_c_Lmol = 2.14107171795 * (cubic->get_Tc()[0] / cubic->get_pc()[0] * 1000) + 0.00773144012514;  // [L/mol]
            return 1 / (v_c_Lmol / 1000.0);
        } else {
            return HelmholtzEOSMixtureBackend::calc_rhomolar_critical();
        }
    };
    CoolPropDbl calc_pmax() override {
        const std::vector<double>& pc = cubic->get_pc();
        if (pc.empty()) {
            throw ValueError("CubicBackend::calc_pmax: Critical pressure vector is empty. Missing cubic parameters.");
        }
        double pc_max = *std::max_element(pc.begin(), pc.end());
        return 100.0 * pc_max;
    };
    CoolPropDbl calc_p_triple() override {
        const std::vector<double>& pc = cubic->get_pc();
        if (pc.empty()) {
            throw ValueError("CubicBackend::calc_p_triple: Critical pressure vector is empty. Missing cubic parameters.");
        }
        double pc_min = *std::min_element(pc.begin(), pc.end());
        return 0.01 * pc_min;
    };
    CoolPropDbl calc_Tmin() override {
        const std::vector<double>& Tc = cubic->get_Tc();
        if (Tc.empty()) {
            throw ValueError("CubicBackend::calc_Tmin: Critical temperature vector is empty. Missing cubic parameters.");
        }
        double Tc_min = *std::min_element(Tc.begin(), Tc.end());
        return 0.3 * Tc_min;
    };
    CoolPropDbl calc_Tmax() override {
        const std::vector<double>& Tc = cubic->get_Tc();
        if (Tc.empty()) {
            throw ValueError("CubicBackend::calc_Tmax: Critical temperature vector is empty. Missing cubic parameters.");
        }
        double Tc_max = *std::max_element(Tc.begin(), Tc.end());
        return 10.0 * Tc_max;
    };
 
    void get_linear_reducing_parameters(double& rhomolar, double& T);
 
    void get_critical_point_starting_values(double& delta0, double& tau0) override;
 
    void get_critical_point_search_radii(double& R_delta, double& R_tau) override;
 
    bool get_critical_is_terminated(double& delta, double& tau) override;
 
    CoolPropDbl calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
                                          const CoolPropDbl& delta) override;
 
    CoolPropDbl calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar);
 
    virtual void update_DmolarT();
 
    void update(CoolProp::input_pairs input_pair, double value1, double value2) override;
 
    void rho_Tp_cubic(CoolPropDbl T, CoolPropDbl p, int& Nsolns, double& rho0, double& rho1, double& rho2);
 
    CoolPropDbl solver_rho_Tp_SRK(CoolPropDbl T, CoolPropDbl p, phases phase) override {
        // Catch only the expected "no root" CoolProp error (a ValueError, which
        // derives from CoolPropBaseError) -- let std::bad_alloc / programming
        // errors propagate rather than silently degrading them to a -1 sentinel.
        try {
            return solver_rho_Tp(T, p, phase);
        } catch (const CoolPropBaseError&) {
            return -1;
        }
    };
    CoolPropDbl solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rho_guess = -1) override;
 
    CoolPropDbl solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, phases phase);
 
    CoolPropDbl solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomax) override;
 
    void update_TPD_state() {
        TPD_state.reset(get_copy());
    };
 
    void saturation(CoolProp::input_pairs inputs);
 
    CoolPropDbl calc_molar_mass() override;
 
    void set_binary_interaction_double(const std::size_t i1, const std::size_t i2, const std::string& parameter, const double value) override;
    double get_binary_interaction_double(const std::size_t i1, const std::size_t i2, const std::string& parameter) override;
 
    void set_binary_interaction_double(const std::string& CAS1, const std::string& CAS2, const std::string& parameter, const double value) override {
        throw ValueError("set_binary_interaction_double not defined for AbstractCubic not defined for CAS #");
    }
    double get_binary_interaction_double(const std::string& CAS1, const std::string& CAS2, const std::string& parameter) override {
        throw ValueError("get_binary_interaction_double not defined for AbstractCubic not defined for CAS #");
    };
 
    // Return a 1-1 copy of this class
    HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true) override = 0;
 
    // Copy the entire kij matrix from another instance in one shot
    void copy_k(AbstractCubicBackend* donor);
 
    //
    void copy_all_alpha_functions(AbstractCubicBackend* donor);
 
    void copy_internals(AbstractCubicBackend& donor);
 
    // Set the cubic alpha function's constants:
    void set_cubic_alpha_C(const size_t i, const std::string& parameter, const double c1, const double c2, const double c3) override;
 
    // Set fluid parameter (currently the volume translation parameter)
    void set_fluid_parameter_double(const size_t i, const std::string& parameter, const double value) override;
 
    // Get fluid parameter (currently the volume translation parameter)
    double get_fluid_parameter_double(const size_t i, const std::string& parameter) override;
 
    virtual int get_superanc_eos_code() const {
        return CubicSuperAncillary::UNKNOWN_CODE;
    }
 
    CoolPropDbl calc_saturation_ancillary(parameters param, int Q, parameters given, double value) override;
 
    void update_QT_pure_superanc(CoolPropDbl Q, CoolPropDbl T) override;
 
    double calc_superanc_Tmax();
};
 
class SRKBackend : public AbstractCubicBackend
{
 
   public:
    SRKBackend(const std::vector<double>& Tc, const std::vector<double>& pc, const std::vector<double>& acentric, double R_u,
               bool generate_SatL_and_SatV = true) {
        cubic = std::make_shared<SRK>(Tc, pc, acentric, R_u);
        setup(generate_SatL_and_SatV);
    };
    SRKBackend(double Tc, double pc, double acentric, double R_u, bool generate_SatL_and_SatV = true) {
        cubic = std::make_shared<SRK>(Tc, pc, acentric, R_u);
        setup(generate_SatL_and_SatV);
    }
    SRKBackend(const std::vector<std::string>& fluid_identifiers, const double R_u = get_config_double(R_U_CODATA),
               bool generate_SatL_and_SatV = true) {
        std::vector<double> Tc, pc, acentric;
        N = fluid_identifiers.size();
        components.resize(N);
        for (std::size_t i = 0; i < fluid_identifiers.size(); ++i) {
            components[i] = CubicLibrary::get_cubic_values(fluid_identifiers[i]);
            Tc.push_back(components[i].Tc);
            pc.push_back(components[i].pc);
            acentric.push_back(components[i].acentric);
        }
        cubic = std::make_shared<SRK>(Tc, pc, acentric, R_u);
        setup(generate_SatL_and_SatV);
    }
    HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true) override {
        AbstractCubicBackend* ACB = new SRKBackend(cubic->get_Tc(), cubic->get_pc(), cubic->get_acentric(), cubic->get_R_u(), generate_SatL_and_SatV);
        ACB->copy_internals(*this);
        return static_cast<HelmholtzEOSMixtureBackend*>(ACB);
    }
    std::string backend_name() override {
        return get_backend_string(SRK_BACKEND);
    }
    int get_superanc_eos_code() const override {
        return CubicSuperAncillary::SRK_CODE;
    }
};
 
class PengRobinsonBackend : public AbstractCubicBackend
{
 
   public:
    PengRobinsonBackend() = default;  // Default constructor (make sure you know what you are doing)
    PengRobinsonBackend(const std::vector<double>& Tc, const std::vector<double>& pc, const std::vector<double>& acentric, double R_u,
                        bool generate_SatL_and_SatV = true) {
        cubic = std::make_shared<PengRobinson>(Tc, pc, acentric, R_u);
        setup(generate_SatL_and_SatV);
    };
    PengRobinsonBackend(double Tc, double pc, double acentric, double R_u, bool generate_SatL_and_SatV = true) {
        cubic = std::make_shared<PengRobinson>(Tc, pc, acentric, R_u);
        setup(generate_SatL_and_SatV);
    };
    PengRobinsonBackend(const std::vector<std::string>& fluid_identifiers, const double R_u = get_config_double(R_U_CODATA),
                        bool generate_SatL_and_SatV = true) {
        std::vector<double> Tc, pc, acentric;
        N = fluid_identifiers.size();
        components.resize(N);
        for (std::size_t i = 0; i < fluid_identifiers.size(); ++i) {
            components[i] = CubicLibrary::get_cubic_values(fluid_identifiers[i]);
            Tc.push_back(components[i].Tc);
            pc.push_back(components[i].pc);
            acentric.push_back(components[i].acentric);
        }
        cubic = std::make_shared<PengRobinson>(Tc, pc, acentric, R_u);
        setup(generate_SatL_and_SatV);
    };
    HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true) override {
        AbstractCubicBackend* ACB =
          new PengRobinsonBackend(cubic->get_Tc(), cubic->get_pc(), cubic->get_acentric(), cubic->get_R_u(), generate_SatL_and_SatV);
        ACB->copy_internals(*this);
        return static_cast<HelmholtzEOSMixtureBackend*>(ACB);
    }
    std::string backend_name() override {
        return get_backend_string(PR_BACKEND);
    }
    int get_superanc_eos_code() const override {
        return CubicSuperAncillary::PR_CODE;
    }
};
 
class CubicResidualHelmholtz : public ResidualHelmholtz
{
 
   protected:
    AbstractCubicBackend* ACB;
 
   public:
    CubicResidualHelmholtz()
      : ACB(nullptr) {
 
        };
    CubicResidualHelmholtz(AbstractCubicBackend* ACB) : ACB(ACB) {};
 
    // copy assignment
    CubicResidualHelmholtz& operator=(CubicResidualHelmholtz& other) {
        ACB = other.ACB;
        return *this;
    }
 
    HelmholtzDerivatives all(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& mole_fractions, double tau, double delta,
                             bool cache_values = false) override {
        HelmholtzDerivatives a;
        // The cubic mixing functions take a std::vector<double>; mole_fractions is a
        // std::vector<CoolPropDbl>.  When CoolPropDbl == double (the default config) the two types
        // are identical, so bind a const reference and skip the per-call heap allocation in the
        // residual hot path; the long-double config still needs the converting copy.
#ifdef COOLPROPDBL_MAPS_TO_DOUBLE
        const std::vector<double>& z = mole_fractions;
#else
        std::vector<double> z = std::vector<double>(mole_fractions.begin(), mole_fractions.end());
#endif
        shared_ptr<AbstractCubic>& cubic = ACB->get_cubic();
 
        // AbstractCubic::alphar(itau, idelta) is
        //     psi_minus(delta, itau, idelta) - tau_times_a(tau, itau)/(R_u*T_r) * psi_plus(delta, idelta),
        // and psi_minus is identically zero for itau > 0.  So all 15 derivatives factor over three
        // small sets of intermediates, each filled below at every derivative order n in [0, 4] and
        // indexed by the order it carries -- note n means a delta-derivative for the two psi arrays
        // but a tau-derivative for tau_times_a:
        //     psi_minus[n]   = psi_minus(delta, itau=0, idelta=n)  -- n-th delta-derivative of psi_minus
        //     psi_plus[n]    = psi_plus(delta, idelta=n)           -- n-th delta-derivative of psi_plus
        //     tau_times_a[n] = tau_times_a(tau, itau=n)            -- n-th tau-derivative of tau*a
        // The 15 separate alphar() calls would otherwise each rebuild these from scratch (the
        // am/bm/psi mixing terms dominated the residual hot path).  Compute each once, then assemble.
        // The assembly is term-for-term identical to alphar() (same operands, same grouping, with the
        // literal 0.0 standing in for psi_minus at itau > 0), so results are bit-for-bit unchanged.
        std::array<double, 5> psi_minus{}, psi_plus{}, tau_times_a{};
        for (std::size_t n = 0; n <= 4; ++n) {
            psi_minus[n] = cubic->psi_minus(delta, z, 0, n);
            psi_plus[n] = cubic->psi_plus(delta, z, n);
            tau_times_a[n] = cubic->tau_times_a(tau, z, n);
        }
        const double den = cubic->get_R_u() * cubic->get_Tr();
        a.alphar = psi_minus[0] - tau_times_a[0] / den * psi_plus[0];
        a.dalphar_dtau = 0.0 - tau_times_a[1] / den * psi_plus[0];
        a.dalphar_ddelta = psi_minus[1] - tau_times_a[0] / den * psi_plus[1];
        a.d2alphar_dtau2 = 0.0 - tau_times_a[2] / den * psi_plus[0];
        a.d2alphar_ddelta_dtau = 0.0 - tau_times_a[1] / den * psi_plus[1];
        a.d2alphar_ddelta2 = psi_minus[2] - tau_times_a[0] / den * psi_plus[2];
        a.d3alphar_dtau3 = 0.0 - tau_times_a[3] / den * psi_plus[0];
        a.d3alphar_ddelta_dtau2 = 0.0 - tau_times_a[2] / den * psi_plus[1];
        a.d3alphar_ddelta2_dtau = 0.0 - tau_times_a[1] / den * psi_plus[2];
        a.d3alphar_ddelta3 = psi_minus[3] - tau_times_a[0] / den * psi_plus[3];
        a.d4alphar_dtau4 = 0.0 - tau_times_a[4] / den * psi_plus[0];
        a.d4alphar_ddelta_dtau3 = 0.0 - tau_times_a[3] / den * psi_plus[1];
        a.d4alphar_ddelta2_dtau2 = 0.0 - tau_times_a[2] / den * psi_plus[2];
        a.d4alphar_ddelta3_dtau = 0.0 - tau_times_a[1] / den * psi_plus[3];
        a.d4alphar_ddelta4 = psi_minus[4] - tau_times_a[0] / den * psi_plus[4];
        return a;
    }
    CoolPropDbl dalphar_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 0, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 0, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 1, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 2, 0, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 1, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 2, i, xN_flag == XN_INDEPENDENT);
    }
 
    CoolPropDbl d2alphardxidxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 0, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 0, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 1, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
 
    CoolPropDbl d4alphar_dxi_dTau3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 3, 0, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dDelta2_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 2, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dDelta_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 2, 1, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dDelta3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 3, i, xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 2, 0, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 1, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 2, i, j,
                                                  xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d3alphardxidxjdxk(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
                                  x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d3_alphar_dxidxjdxk(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 0, i, j, k,
                                                     xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dxj_dxk_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
                                          x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d3_alphar_dxidxjdxk(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 1, 0, i, j, k,
                                                     xN_flag == XN_INDEPENDENT);
    }
    CoolPropDbl d4alphar_dxi_dxj_dxk_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
                                            x_N_dependency_flag xN_flag) override {
        return ACB->get_cubic()->d3_alphar_dxidxjdxk(HEOS.tau(), HEOS.delta(), HEOS.get_mole_fractions_doubleref(), 0, 1, i, j, k,
                                                     xN_flag == XN_INDEPENDENT);
    }
};
 
} /* namespace CoolProp */
#endif /* CUBICBACKEND_H_ */