#include "NGCal.h"
#include "Detail.h"
#include "Therm.h"
#include <math.h>
#include <stdbool.h> 
static void CprCvrHS(NGParSTRUCT *ptNGPar, Detail *ptD);
static double H(NGParSTRUCT *ptNGPar, Detail *ptD);
static double S(NGParSTRUCT *ptNGPar, Detail *ptD);
static void HS_Mode(NGParSTRUCT *ptNGPar, Detail *ptD, double H_target, double S_target, bool bGuess);
double  CalTH = 4.1840;
  double GK_root[5] = {
    0.14887433898163121088,
    0.43339539412924719080,
    0.67940956829902440263,
    0.86506336668898451073, 
    0.97390652851717172008 
};
  double GK_weight[5] = {
    0.29552422471475286217, 
    0.26926671930999634918, 
    0.21908636251598204295, 
    0.14945134915058059038, 
    0.066671344308688137179 
};
  int GK_points = 5;
  double ThermConstants[NUMBEROFCOMPONENTS][11] = {
    {-29776.4, 7.95454, 43.9417, 1037.09, 1.56373, 813.205, -24.9027, 1019.98,-10.1601, 1070.14,-20.0615},
    {-3495.34, 6.95587, 0.272892, 662.738,-0.291318,-680.562, 1.78980, 1740.06, 0.0, 100.0, 4.49823},
    { 20.7307, 6.96237, 2.68645, 500.371,-2.56429,-530.443, 3.91921, 500.198, 2.13290, 2197.22, 5.81381},
    {-37524.4, 7.98139, 24.3668, 752.320, 3.53990, 272.846, 8.44724, 1020.13,-13.2732, 869.510,-22.4010},
    {-56072.1, 8.14319, 37.0629, 735.402, 9.38159, 247.190, 13.4556, 1454.78,-11.7342, 984.518,-24.0426},
    {-13773.1, 7.97183, 6.27078, 2572.63, 2.05010, 1156.72, 0.0, 100.0, 0.0, 100.0, -3.24989},
    {-10085.4, 7.94680,-0.08380, 433.801, 2.85539, 843.792, 6.31595, 1481.43,-2.88457, 1102.23,-0.51551},
    {-5565.60, 6.66789, 2.33458, 2584.98, .749019, 559.656, 0.0, 100.0, 0.0, 100.0, -7.94821},
    {-2753.49, 6.95854, 2.02441, 1541.22, .096774, 3674.81, 0.0, 100.0, 0.0, 100.0, 6.23387},
    {-3497.45, 6.96302, 2.40013, 2522.05, 2.21752, 1154.15, 0.0, 100.0, 0.0, 100.0, 9.19749},
    {-72387.0, 17.8143, 58.2062, 1787.39, 40.7621, 808.645, 0.0, 100.0, 0.0, 100.0, -44.1341},
    {-72674.8, 18.6383, 57.4178, 1792.73, 38.6599, 814.151, 0.0, 100.0, 0.0, 100.0, -46.1938},
    {-91505.5, 21.3861, 74.3410, 1701.58, 47.0587, 775.899, 0.0, 100.0, 0.0, 100.0, -60.2474},
    {-83845.2, 22.5012, 69.5789, 1719.58, 46.2164, 802.174, 0.0, 100.0, 0.0, 100.0, -62.2197},
    {-94982.5, 26.6225, 80.3819, 1718.49, 55.6598, 802.069, 0.0, 100.0, 0.0, 100.0, -77.5366},
    {-103353., 30.4029, 90.6941, 1669.32, 63.2028, 786.001, 0.0, 100.0, 0.0, 100.0, -92.0164},
    {-109674., 34.0847, 100.253, 1611.55, 69.7675, 768.847, 0.0, 100.0, 0.0, 100.0, -106.149},
    {-122599., 38.5014, 111.446, 1646.48, 80.5015, 781.588, 0.0, 100.0, 0.0, 100.0, -122.444},
    {-133564., 42.7143, 122.173, 1654.85, 90.2255, 785.564, 0.0, 100.0, 0.0, 100.0, -138.006},
    { 0.0, 4.9680, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0 },
    { 0.0, 4.9680, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0 }
};
static double dSi = 0.0;
static double dTold = 0.0;
static double dMrxold = 0.0;

void ThermInit(void)
{
        dSi = 0.0;
    dTold = 0.0;
    dMrxold = 0.0;
}

void ThermDestroy(void)
{
    }

double coth(double x)
{
    return cosh(x)/sinh(x);
}

void Run(NGParSTRUCT *ptNGPar, Detail *ptD)
{
        double c, x, y, z;
        Detail_Run(ptD, ptNGPar);
        Detail_dZdD(ptD, ptNGPar->dDf);
        CprCvrHS(ptNGPar, ptD);
        ptNGPar->dk = ptNGPar->dCp / ptNGPar->dCv;
        x = ptNGPar->dk * RGAS * 1000.0 * ptNGPar->dTf;
    y = ptNGPar->dMrx;
    z = ptNGPar->dZf + ptNGPar->dDf * ptD->ddZdD;
        c = (x / y) * z;
        ptNGPar->dSOS = sqrt(c);
            ptNGPar->dKappa = (c * ptNGPar->dRhof) / ptNGPar->dPf;
}

double CpiMolar(NGParSTRUCT *ptNGPar)
{
    double Cp = 0.0;
    double Cpx;
    double DT, FT, HT, JT;
    double Dx, Fx, Hx, Jx;
    double T;
    int i;
        T = ptNGPar->dTf;
        for (i = 0; i < NUMBEROFCOMPONENTS; i++)
    {
                if (ptNGPar->adMixture[i] <= 0.0) continue;
                Cpx = 0.0;
                DT = ThermConstants[i][coefD] / T;
        FT = ThermConstants[i][coefF] / T;
        HT = ThermConstants[i][coefH] / T;
        JT = ThermConstants[i][coefJ] / T;
                Dx = DT/sinh(DT);
        Fx = FT/cosh(FT);
        Hx = HT/sinh(HT);
        Jx = JT/cosh(JT);
        Cpx += ThermConstants[i][coefB];
        Cpx += ThermConstants[i][coefC] * Dx * Dx;
        Cpx += ThermConstants[i][coefE] * Fx * Fx;
        Cpx += ThermConstants[i][coefG] * Hx * Hx;
        Cpx += ThermConstants[i][coefI] * Jx * Jx;
                Cpx *= ptNGPar->adMixture[i];
                Cp += Cpx;
    }
        Cp *= CalTH;
    return Cp;
}

double Ho(NGParSTRUCT *ptNGPar)
{
    double H = 0.0;
    double Hx;
    double DT, FT, HT, JT;
    double cothDT, tanhFT, cothHT, tanhJT;
    double T;
    int i;
        T = ptNGPar->dTf;
    for (i = 0; i < NUMBEROFCOMPONENTS; i++)
    {
                if (ptNGPar->adMixture[i] <= 0.0) continue;
        Hx = 0.0;
                DT = ThermConstants[i][coefD] / T;
        FT = ThermConstants[i][coefF] / T;
        HT = ThermConstants[i][coefH] / T;
        JT = ThermConstants[i][coefJ] / T;
        cothDT = coth(DT);
        tanhFT = tanh(FT);
        cothHT = coth(HT);
        tanhJT = tanh(JT);
        Hx += ThermConstants[i][coefA];
        Hx += ThermConstants[i][coefB] * T;
        Hx += ThermConstants[i][coefC] * ThermConstants[i][coefD] * cothDT;
        Hx -= ThermConstants[i][coefE] * ThermConstants[i][coefF] * tanhFT;
        Hx += ThermConstants[i][coefG] * ThermConstants[i][coefH] * cothHT;
        Hx -= ThermConstants[i][coefI] * ThermConstants[i][coefJ] * tanhJT;
                Hx *= ptNGPar->adMixture[i];
                H += Hx;
    }
        H *= CalTH;
        H /= ptNGPar->dMrx;
        return H * 1.e3;
}

double So(NGParSTRUCT *ptNGPar)
{
    double S = 0.0;
    double Sx;
    double DT, FT, HT, JT;
    double cothDT, tanhFT, cothHT, tanhJT;
    double sinhDT, coshFT, sinhHT, coshJT;
    double T;
    int i;
        T = ptNGPar->dTf;
    for (i = 0; i < NUMBEROFCOMPONENTS; i++)
    {
                if (ptNGPar->adMixture[i] <= 0.0) continue;
        Sx = 0.0;
                DT = ThermConstants[i][coefD] / T;
        FT = ThermConstants[i][coefF] / T;
        HT = ThermConstants[i][coefH] / T;
        JT = ThermConstants[i][coefJ] / T;
        cothDT = coth(DT);
        tanhFT = tanh(FT);
        cothHT = coth(HT);
        tanhJT = tanh(JT);
        sinhDT = sinh(DT);
        coshFT = cosh(FT);
        sinhHT = sinh(HT);
        coshJT = cosh(JT);
        Sx += ThermConstants[i][coefK];
        Sx += ThermConstants[i][coefB] * log(T);
        Sx += ThermConstants[i][coefC] * (DT * cothDT - log(sinhDT));
        Sx -= ThermConstants[i][coefE] * (FT * tanhFT - log(coshFT));
        Sx += ThermConstants[i][coefG] * (HT * cothHT - log(sinhHT));
        Sx -= ThermConstants[i][coefI] * (JT * tanhJT - log(coshJT));
                Sx *= ptNGPar->adMixture[i];
                S += Sx;
    }
        S *= CalTH;
        S /= ptNGPar->dMrx;
        return S * 1.e3;
}

void CprCvrHS(NGParSTRUCT *ptNGPar, Detail *ptD)
{
    double Cvinc, Cvr, Cpr;
    double Hinc;
    double Sinc;
    double Smixing;
    double Cp, Si;
    double a, b, x;
    int i;
        Cvinc = 0.0;
    Hinc = 0.0;
    Sinc = 0.0;
        Smixing = 0.0;
        Cp = CpiMolar(ptNGPar);
        ptNGPar->dHo = Ho(ptNGPar);
        Si = So(ptNGPar);
        ptNGPar->dCpi = (Cp * 1000.0) / ptNGPar->dMrx;
        for (i = 0; i < GK_points; i++)
    {
                x = ptNGPar->dDf * (1.0 + GK_root[i]) / 2.0;
                Detail_zdetail(ptD, x);
        Detail_dZdT(ptD, x);
        Detail_d2ZdT2(ptD, x);
                Hinc += GK_weight[i] * ptD->ddZdT / x;
        Cvinc += GK_weight[i] * (2.0 * ptD->ddZdT + ptNGPar->dTf * ptD->dd2ZdT2) / x;
        Sinc += GK_weight[i] * (ptD->dZ + ptNGPar->dTf * ptD->ddZdT - 1.0) / x;
                x = ptNGPar->dDf * (1.0 - GK_root[i]) / 2.0;
                Detail_zdetail(ptD, x);
                Detail_dZdT(ptD, x);
        Detail_d2ZdT2(ptD, x);
                Hinc += GK_weight[i] * ptD->ddZdT / x;
        Cvinc += GK_weight[i] * (2.0 * ptD->ddZdT + ptNGPar->dTf * ptD->dd2ZdT2) / x;
        Sinc += GK_weight[i] * (ptD->dZ + ptNGPar->dTf * ptD->ddZdT - 1.0) / x;
    }
        Detail_zdetail(ptD, ptNGPar->dDf);
    Detail_dZdT(ptD, ptNGPar->dDf);
    Detail_d2ZdT2(ptD, ptNGPar->dDf);
        Cvr = Cp - RGAS * (1.0 + ptNGPar->dTf * Cvinc * 0.5 * ptNGPar->dDf);
        a = (ptNGPar->dZf + ptNGPar->dTf * ptD->ddZdT);
    b = (ptNGPar->dZf + ptNGPar->dDf * ptD->ddZdD);
        Cpr = Cvr + RGAS * ((a * a)/b);
        Cpr /= ptNGPar->dMrx;
    Cvr /= ptNGPar->dMrx;
        ptNGPar->dCv = Cvr * 1000.0;     ptNGPar->dCp = Cpr * 1000.0;
        ptNGPar->dH = ptNGPar->dHo + 1000.0 * RGAS * ptNGPar->dTf *
                 (ptNGPar->dZf - 1.0 - ptNGPar->dTf * Hinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx;
        for (i = 0; i < NUMBEROFCOMPONENTS; i++)
    {
        if (ptNGPar->adMixture[i] > 0.0)              Smixing += ptNGPar->adMixture[i] * log(ptNGPar->adMixture[i]);
    }
    Smixing *= RGAS;
        ptNGPar->dS = Si - Smixing - 1000.0 * RGAS * 
                 (log(ptNGPar->dPf/101325.0) - log(ptNGPar->dZf) + Sinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx;
}

void HS_Mode(NGParSTRUCT *ptNGPar, Detail *ptD, double H_target, double S_target, bool bGuess)
{
    double s0, s1, s2, t0, t1, t2, tmin, tmax;
    double h0, h1, h2, p0, p1, p2, px, pmin, pmax;
    double delta1, delta2;
    double tolerance = 0.001;     int i, j;
        s0 = S_target;
    h0 = H_target;
        if (bGuess)
    {
        t1 = ptNGPar->dTf;
        px = ptNGPar->dPf;
        pmax = px * 2.0;
        pmin = px * 0.1;
        tmax = t1 * 1.5;
        tmin = t1 * 0.67;
    }
    else     {
        t1 = 273.15;
        px = 1013250.0;         pmax = P_MAX;
        pmin = 10000.0;         tmax = T_MAX;
        tmin = T_MIN;
    }
        t2 = t1 + 10.0;
            Detail_Run(ptD, ptNGPar);
        h1 = H(ptNGPar, ptD) - h0;
        for (i = 0; i < MAX_NUM_OF_ITERATIONS; i++)
    {
        ptNGPar->dTf = t2;
        p1 = px;         p2 = px * 0.1;         ptNGPar->dPf = p1;
        Detail_Run(ptD, ptNGPar);
        s1 = S(ptNGPar, ptD) - s0;
                for (j = 0; j < MAX_NUM_OF_ITERATIONS; j++)
        {
            ptNGPar->dPf = p2;
            Detail_Run(ptD, ptNGPar);
            s2 = S(ptNGPar, ptD) - s0;
                        delta2 = fabs(s1 - s2) / s0;
                        if (delta2 < tolerance) break;
                        p0 = p2;
            p2 = (p1 * s2 - p2 * s1) / (s2 - s1);
                        if (p2 <= pmin)
            {
                p2 = pmin;
            }
                        if (p2 >= pmax) p2 = pmax;
                        p1 = p0;
            s1 = s2;
        }
                if (j >= MAX_NUM_OF_ITERATIONS)
            ptNGPar->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED;
                h2 = H(ptNGPar, ptD) - h0;
                delta1 = fabs(h1 - h2) / h0;
                if (delta1 < tolerance && i > 0) break;
                t0 = t2;
        t2 = (t1 * h2 - t2 * h1) / (h2 - h1);
                if (t2 >= tmax) t2 = tmax;
                if (t2 <= tmin)
        {
            t2 = t0 + 10.0;
            ptNGPar->dTf = t2;
            Detail_Run(ptD, ptNGPar);
            h2 = H(ptNGPar, ptD) - h0;
        }
        t1 = t0;
        h1 = h2;
    }
        if (i >= MAX_NUM_OF_ITERATIONS)
        ptNGPar->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED;
}

double H(NGParSTRUCT *ptNGPar, Detail *ptD)
{
    double Hinc;
    double x;
    int i;
        Hinc = 0.0;
        ptNGPar->dHo = Ho(ptNGPar);
        for (i = 0; i < GK_points; i++)
    {
                x = ptNGPar->dDf * (1.0 + GK_root[i]) / 2.0;
        Detail_zdetail(ptD, x);
        Detail_dZdT(ptD, x);
        Detail_d2ZdT2(ptD, x);
        Hinc += GK_weight[i] * ptD->ddZdT / x;
        if (i == 10) break;
        x = ptNGPar->dDf * (1.0 - GK_root[i]) / 2.0;
        Detail_zdetail(ptD, x);
        Detail_dZdT(ptD, x);
        Detail_d2ZdT2(ptD, x);
        Hinc += GK_weight[i] * ptD->ddZdT / x;
    }
    Detail_zdetail(ptD, ptNGPar->dDf);
    Detail_dZdT(ptD, ptNGPar->dDf);
    Detail_d2ZdT2(ptD, ptNGPar->dDf);
        ptNGPar->dH = ptNGPar->dHo + 1000.0 * RGAS * ptNGPar->dTf *
                 (ptNGPar->dZf - 1.0 - ptNGPar->dTf * Hinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx;
    return ptNGPar->dH;
}

double S(NGParSTRUCT *ptNGPar, Detail *ptD)
{
    double Sinc;
    double Smixing;
    double x;
    int i;
        Sinc = 0.0;
        Smixing = 0.0;
        for (i = 0; i < GK_points; i++)
    {
                x = ptNGPar->dDf * (1.0 + GK_root[i]) / 2.0;
        Detail_zdetail(ptD, x);
        Detail_dZdT(ptD, x);
        Detail_d2ZdT2(ptD, x);
        Sinc += GK_weight[i] * (ptD->dZ + ptNGPar->dTf * ptD->ddZdT - 1.0) / x;
        if (i == 10) break;
        x = ptNGPar->dDf * (1.0 - GK_root[i]) / 2.0;
        Detail_zdetail(ptD, x);
        Detail_dZdT(ptD, x);
        Detail_d2ZdT2(ptD, x);
        Sinc += GK_weight[i] * (ptD->dZ + ptNGPar->dTf * ptD->ddZdT - 1.0) / x;
    }
        Detail_zdetail(ptD, ptNGPar->dDf);
    Detail_dZdT(ptD, ptNGPar->dDf);
    Detail_d2ZdT2(ptD, ptNGPar->dDf);
        if (ptNGPar->dTf != dTold || ptNGPar->dMrx != dMrxold)
    {
        dSi = So(ptNGPar);
        dTold = ptNGPar->dTf;
        dMrxold = ptNGPar->dMrx;
    }
        for (i = 0; i < NUMBEROFCOMPONENTS; i++)
    {
        if (ptNGPar->adMixture[i] > 0.0)              Smixing += ptNGPar->adMixture[i] * log(ptNGPar->adMixture[i]);
    }
    Smixing *= RGAS;
        ptNGPar->dS = dSi - Smixing - 1000.0 * RGAS * 
                 (log(ptNGPar->dPf/101325.0) - log(ptNGPar->dZf) + Sinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx;
    return ptNGPar->dS;
}