#ifndef ANCILLARIES_H
#define ANCILLARIES_H

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

namespace CoolProp{

/**
The surface tension correlation class uses correlations for the surface tension that are all
of the form

\f[
\sigma = \sum_i a_i\left(1-\frac{T}{\tilde T_c}\right)^{n_i}
\f]

where \f$ \tilde T_c \f$ is the critical temperature used for the correlation which is
almost always equal to the critical temperature of the equation of state.  Result for
surface tension is in N/m
*/
class SurfaceTensionCorrelation
{
public:
    std::vector<CoolPropDbl> a, ///< the leading coefficients a_i
                             n, ///< the powers n_i
                             s; ///< a summation buffer
    CoolPropDbl Tc; ///< critical temperature in K
    std::size_t N; ///< number of a_i, n_i pairs
    std::string BibTeX; ///< The BiBTeX key for the surface tension curve in use

    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;
    };
    /// Actually evaluate the surface tension equation
    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);
    }
};
/**
 * 
 * This is generalized class that can be used to manage an ancillary curve, 
 * here they are ancillary curves for saturation pressure, density, enthalpy, entropy.
 * 
 * The form of the ancillary equation can take one of a number of forms:
 * 
 * a) So-called "exponential" form (type = TYPE_EXPONENTIAL) that has a form like
 * 
 * \f[ y = y_c\exp\left(\frac{T_c}{T}\sum(n_i \theta^{t_i})\right) \f] 
 * or 
 * \f[ y = y_c\exp\left(\sum(n_i \theta^{t_i})\right) \f]
 * 
 * b) A non-exponential form (type = TYPE_NOT_EXPONENTIAL) that has a form of
 *
 * \f[ y = y_c\left(1+\sum_i(n_i\theta^t_i)\right) \f]
 * with
 * \f[ \theta = \left(1-\frac{T}{T_c}\right) \f]
 * which is conveniently equal to zero at the critical point
 * 
 * c) Rational polynomial form (type = TYPE_RATIONAL_POLYNOMIAL) that has a form of
 * \f[ y = \frac{\sum_iA_iT^i}{\sum_iB_iT^i}\f]
 * where i is an integer, and the coefficients are in increasing order in both numerator and denominator
*/
class SaturationAncillaryFunction
{
private:
    Eigen::MatrixXd num_coeffs, ///< Coefficients for numerator in rational polynomial 
                    den_coeffs; ///< Coefficients for denominator in rational polynomial
    std::vector<double> n, t, s; // For TYPE_NOT_EXPONENTIAL & TYPE_EXPONENTIAL
    union{
        CoolPropDbl max_abs_error; ///< For TYPE_RATIONAL_POLYNOMIAL
        struct{                    // For TYPE_NOT_EXPONENTIAL & TYPE_EXPONENTIAL
            bool using_tau_r; ///< Whether the term \f$ \frac{T_c}{T} \f$ is included in the 
            CoolPropDbl reducing_value, ///< The value used to reduce the output variable
                        T_r; ///< The temperature in K used to reduce the temperature (usually the critical temperature)
            std::size_t N; ///< The number of values in the arrays
        };
    };
    CoolPropDbl Tmax, ///< The maximum temperature in K
                Tmin; ///< The minimum temperature in K
    enum ancillaryfunctiontypes{TYPE_NOT_SET = 0, 
                                TYPE_NOT_EXPONENTIAL, ///< It is a non-exponential type of equation
                                TYPE_EXPONENTIAL, ///< It is an exponential type equation, with or without the T_c/T term
                                TYPE_RATIONAL_POLYNOMIAL ///< It is a rational polynomial equation
                                };
    ancillaryfunctiontypes type; ///< The type of ancillary curve being used
public:

    SaturationAncillaryFunction(){type = TYPE_NOT_SET; Tmin = _HUGE; Tmax = _HUGE;};
    SaturationAncillaryFunction(rapidjson::Value &json_code);
    
    /// Return true if the ancillary is enabled (type is not TYPE_NOT_SET)
    bool enabled(void){return type != TYPE_NOT_SET;}
    
    /// Get the maximum absolute error for this fit
    /// @returns max_abs_error the maximum absolute error for ancillaries that are characterized by maximum absolute error
    CoolPropDbl get_max_abs_error(){return max_abs_error;};
    
    /// Evaluate this ancillary function, yielding for instance the saturated liquid density
    /// @param T The temperature in K
    /// @returns y the value of the ancillary function at temperature T
    double evaluate(double T);
    
    /// Invert this ancillary function, and calculate the temperature given the output the value of the function
    /// @param value The value of the output
    /// @param min_bound (optional) The minimum value for T; ignored if < 0
    /// @param max_bound (optional) The maximum value for T; ignored if < 0
    /// @returns T The temperature in K
    double invert(double value, double min_bound = -1, double max_bound = -1);
    
    /// Get the minimum temperature in K
    double get_Tmin(void){return Tmin;};
    
    /// Get the maximum temperature in K
    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;
};

/**
\brief The evaluator class for a melting curve formed of segments in the form

\f[
a_i((\frac{T}{T_0})^{t_i}-1)
\f]
*/
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;
};

/** 
 \brief The evaluator class for a melting curve formed of segments in the form

 \f[
 a_i(\frac{T}{T_0}-1)^t_i
 \f]
*/
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, ///< A simon-type curve is in use
        MELTING_LINE_POLYNOMIAL_IN_TR_TYPE, ///< a polynomial in \f$ T/T_c \f$ is in use
        MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE, ///< a polynomial in \f$ \theta \f$ is in use
    };
    CoolPropDbl Tmin, ///< Minimum temperature in K
                Tmax, ///< Maximum temperature in K
                pmin, ///< Minimum pressure in Pa
                pmax; ///< Maximum pressure in Pa

    std::string BibTeX; ///< BibTeX key for the melting curve in use
    CoolPropDbl T_m; ///< Melting temperature at 1 atmosphere
    MeltingLinePiecewiseSimonData simon; /// The data used for a Simon-style curve
    MeltingLinePiecewisePolynomialInTrData polynomial_in_Tr; /// The data needed for a melting curve formed of segments that are polynomials in \f$ T/T_c \f$
    MeltingLinePiecewisePolynomialInThetaData polynomial_in_Theta; /// The data needed for a melting curve formed of segments that are polynomials in \f$ \theta \f$
    int type;
    
    MeltingLineVariables() :Tmin(_HUGE), Tmax(_HUGE), pmin(_HUGE), pmax(_HUGE), T_m(_HUGE), type(MELTING_LINE_NOT_SET){};

    /** 
     * \brief Evaluate the melting line
     * @param OF The output variable
     * @param GIVEN The given variable
     * @param value The value of the given variable
     */
    CoolPropDbl evaluate(int OF, int GIVEN, CoolPropDbl value);
    
    /// Evaluate the melting line to calculate the limits of the curve (Tmin/Tmax and pmin/pmax)
    void set_limits();
    
    /// Return true if the ancillary is enabled (type is not the default value of MELTING_LINE_NOT_SET)
    bool enabled(){return type != MELTING_LINE_NOT_SET;};
};

} /* namespace CoolProp */
#endif