Ancillaries.h

#ifndef ANCILLARIES_H
#define ANCILLARIES_H

#include "Exceptions.h"
#include <vector>
#include "rapidjson_include.h"
#include "Eigen/Core"
#include "PolyMath.h"

namespace CoolProp {

class SurfaceTensionCorrelation
{
   public:
    std::vector<CoolPropDbl> a,  
      n,                         
      s;                         
    CoolPropDbl Tc;              
    std::size_t N;               
    std::string BibTeX;          

    SurfaceTensionCorrelation() : Tc(_HUGE), N(0) {}
    SurfaceTensionCorrelation(rapidjson::Value& json_code) {
        a = cpjson::get_long_double_array(json_code["a"]);
        n = cpjson::get_long_double_array(json_code["n"]);

        Tc = cpjson::get_double(json_code, "Tc");
        BibTeX = cpjson::get_string(json_code, "BibTeX");

        this->N = n.size();
        s = n;
    };
    CoolPropDbl evaluate(CoolPropDbl T) {
        if (a.empty()) {
            throw NotImplementedError(format("surface tension curve not provided"));
        }
        if (T > Tc) {
            throw ValueError(format("Must be saturated state : T <= Tc"));
        }
        CoolPropDbl THETA = 1 - T / Tc;
        for (std::size_t i = 0; i < N; ++i) {
            s[i] = a[i] * pow(THETA, n[i]);
        }
        return std::accumulate(s.begin(), s.end(), 0.0);
    }
};
class SaturationAncillaryFunction
{
   private:
    Eigen::MatrixXd num_coeffs,   
      den_coeffs;                 
    std::vector<double> n, t, s;  // For TYPE_NOT_EXPONENTIAL & TYPE_EXPONENTIAL
    union
    {
        CoolPropDbl max_abs_error;  
        struct
        {                                // For TYPE_NOT_EXPONENTIAL & TYPE_EXPONENTIAL
            bool using_tau_r;            
            CoolPropDbl reducing_value,  
              T_r;                       
            std::size_t N;               
        };
    };
    CoolPropDbl Tmax,  
      Tmin;            
    enum ancillaryfunctiontypes
    {
        TYPE_NOT_SET = 0,
        TYPE_NOT_EXPONENTIAL,     
        TYPE_EXPONENTIAL,         
        TYPE_RATIONAL_POLYNOMIAL  
    };
    ancillaryfunctiontypes type;  
   public:
    SaturationAncillaryFunction() {
        type = TYPE_NOT_SET;
        Tmin = _HUGE;
        Tmax = _HUGE;
    };
    SaturationAncillaryFunction(rapidjson::Value& json_code);

    bool enabled(void) {
        return type != TYPE_NOT_SET;
    }

    CoolPropDbl get_max_abs_error() {
        return max_abs_error;
    };

    double evaluate(double T);

    double invert(double value, double min_bound = -1, double max_bound = -1);

    double get_Tmin(void) {
        return Tmin;
    };

    double get_Tmax(void) {
        return Tmax;
    };
};

// ****************************************************************************
// ****************************************************************************
//                                 MELTING LINE
// ****************************************************************************
// ****************************************************************************

struct MeltingLinePiecewiseSimonSegment
{
    CoolPropDbl T_0, a, c, p_0, T_max, T_min, p_min, p_max;
};
struct MeltingLinePiecewiseSimonData
{
    std::vector<MeltingLinePiecewiseSimonSegment> parts;
};

class MeltingLinePiecewisePolynomialInTrSegment
{
   public:
    std::vector<CoolPropDbl> a, t;
    CoolPropDbl T_0, p_0, T_max, T_min, p_min, p_max;
    CoolPropDbl evaluate(CoolPropDbl T) {
        CoolPropDbl summer = 0;
        for (std::size_t i = 0; i < a.size(); ++i) {
            summer += a[i] * (pow(T / T_0, t[i]) - 1);
        }
        return p_0 * (1 + summer);
    }
};
struct MeltingLinePiecewisePolynomialInTrData
{
    std::vector<MeltingLinePiecewisePolynomialInTrSegment> parts;
};

class MeltingLinePiecewisePolynomialInThetaSegment
{
   public:
    std::vector<CoolPropDbl> a, t;
    CoolPropDbl T_0, p_0, T_max, T_min, p_min, p_max;

    CoolPropDbl evaluate(CoolPropDbl T) {
        CoolPropDbl summer = 0;
        for (std::size_t i = 0; i < a.size(); ++i) {
            summer += a[i] * pow(T / T_0 - 1, t[i]);
        }
        return p_0 * (1 + summer);
    }
};
struct MeltingLinePiecewisePolynomialInThetaData
{
    std::vector<MeltingLinePiecewisePolynomialInThetaSegment> parts;
};

class MeltingLineVariables
{
   public:
    enum MeltingLineVariablesEnum
    {
        MELTING_LINE_NOT_SET = 0,
        MELTING_LINE_SIMON_TYPE,                
        MELTING_LINE_POLYNOMIAL_IN_TR_TYPE,     
        MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE,  
    };
    CoolPropDbl Tmin,  
      Tmax,            
      pmin,            
      pmax;            

    std::string BibTeX;                   
    CoolPropDbl T_m;                      
    MeltingLinePiecewiseSimonData simon;  
    MeltingLinePiecewisePolynomialInTrData
      polynomial_in_Tr;  
    MeltingLinePiecewisePolynomialInThetaData
      polynomial_in_Theta;  
    int type;

    MeltingLineVariables() : Tmin(_HUGE), Tmax(_HUGE), pmin(_HUGE), pmax(_HUGE), T_m(_HUGE), type(MELTING_LINE_NOT_SET){};

    CoolPropDbl evaluate(int OF, int GIVEN, CoolPropDbl value);

    void set_limits();

    bool enabled() {
        return type != MELTING_LINE_NOT_SET;
    };
};

} /* namespace CoolProp */
#endif