#include <cmath>

//Define ln(gamma_ftn(xx))
inline double gammaln(double XX)
{
        double STP = 2.5066282746310005;
        const double HALF = 0.5, ONE = 1, FPF = 5.5;
        double X, Y, TMP, SER, ga;

        static double COF[] = {
        76.18009172947146,
        -86.50532032941677,
        24.01409824083091,
        -1.231739572450155,
        0.001208650973866179,
        -0.000005395239384953
        };

        X=XX;
        Y=X;
        TMP = X + FPF;
        TMP = (X+HALF)*std::log(TMP)-TMP;
        SER = 1.000000000190015;

        for(int j = 0; j <= 5 ; j++)
        {
                Y = Y + ONE;
                SER = SER+COF[j] / Y ;
        }

        ga = TMP + std::log(STP*SER/X) ;

        return ga;
}

//Define exp(-2*(erf^-1(2*eta - 1))^2)
inline double AMC(double eta) //from kiva
{

        const double pi = 3.141592;
        const double spi = std::sqrt(pi);
        double a1=0.5, a0=(2*eta - 1), slope, da=0, result=0;
        if(eta < 0.000001 || eta > 0.999999)
        {
                result = 1e-30;
        }
        else
        {
AMC_while:;
                slope = 2/spi*std::exp(-1* std::pow(a1,2) );
                da = (a0 - Foam::erf(a1))/slope;
                a1 = a1+da;
                if( std::fabs(da) > 1e-12 ) //abs function for double variable
                {
                        goto AMC_while;
                }
                result = std::exp(-2* std::pow(a1,2) );
        }
        return result;
}

//Define TDMA function
inline void TDMA(scalarField& a, scalarField& b, scalarField& c, scalarField& Qi, scalarField& QiNew, scalar Numofetaspace)
{
        scalar bet;
        scalarField gam(Numofetaspace,0);

        bet = b[0];
        QiNew[0] = Qi[0]/bet;

        for(label j=1 ; j< Numofetaspace ; )
        {
                gam[j] = c[j-1]/bet;
                bet = b[j] - a[j]*gam[j];
                QiNew[j] = ( Qi[j] - a[j]*QiNew[j-1] ) / bet;
                j = j+1;
        }

        for(label j = Numofetaspace - 2 ; j >= 0 ;)
        {
                QiNew[j] = QiNew[j] - gam[j+1] * QiNew[j+1];
                j = j-1;
        }
}

inline double simpson(scalar min, scalar max, scalar interval, scalarField& ftn)
{
        double result_odd=0, result_even=0, result=0;
        for(label i=min+1 ; i < max ; i=i+2)
        {
                //odd points
                result_odd += ftn[i];
        }
        result_odd = result_odd*(4.0/3.0)*interval;

        for(label i=min+2 ; i < max ; i=i+2)
        {
                //even points
                result_even += ftn[i];
        }
        result_even = result_even*(2.0/3.0)*interval;

        //start and end points
        result = result_odd + result_even + (ftn[min]+ftn[max])*interval*(1.0/3.0);

        return result;
}

inline double integration(scalar deltaftn, scalarField& MFcut, scalarField& Neta, scalarField& pdf, scalarField& f)
{
    double result = 0.0;
    double min, max, interval;
    if(deltaftn < 0.5) //one delta function (deltaftn == 0)
    {
        for(label i=0; i<= (Neta[MFcut.size()-1]);i++)
        {
            if(pdf[i] > 1e-30)
            {
                result = f[i];
            }
        }
    }
    else if(deltaftn > 0.5 && deltaftn < 1.5) //two delta function (deltaftn == 1)
    {
        for(label i=0; i<= (Neta[MFcut.size()-1]);i++)
        {
            if(pdf[i] > 1e-30)
            {
                result += 0.5*f[i];
            }
        }
    }
    else if(deltaftn > 1.5) //not a delta function (deltaftn == 2)
    {
        scalarField ftn = pdf*f;
        for(label i=0; i<(MFcut.size()-1); i++ )
        {
            min = Neta[i];
            max = Neta[i+1];
            interval = (MFcut[i+1]-MFcut[i])/(max-min);
            result += simpson(min, max, interval, ftn);
        }
    }
    else
    {
        //Integration error
    }
    return result;
}

//inline double trapz()
//{
//}

//inline double rectang()
//{
//}

//Returns lognormal PDF of the Nst
//Ref) S.H.Kim and Kang Y. Huh, Combust Flame, vol.138, pp.336-352, 2004
//     Equation 17
inline double lognormalPDF(scalar Nst, scalar Nst_local)
{
    double result, numerator;
    const double pi = 3.141592;
    const double sigmaN = 1.0;

    numerator = sqr(std::log(Nst/Nst_local)+0.5*sqr(sigmaN));
    result = std::exp(-0.5*numerator/sqr(sigmaN));
    result = result/std::sqrt(2.0*pi)/Nst/sigmaN;

    return result;
}

//Returns the integral of the lognormalPDF function between a and b, by 
//ten-point Gauss-Legendre integration, Numerical Recipe 2nd Ed. p.140
inline double qgaus(scalar a, scalar b, scalar Nst_local)
{
    double dx, xm, xr, result;
    static double w[] = {
        0.2955242247,
        0.2692667193,
        0.2190863625,
        0.1494513491,
        0.0666713443
    };
    static double x[] = {
        0.1488743389,
        0.4333953941,
        0.6794095682,
        0.8650633666,
        0.9739065285
    };

    xm = 0.5*(b+a);
    xr = 0.5*(b-a);
    result = 0.0;

    for(label i=0 ; i<5 ; i++)
    {
        dx = xr*x[i];
        result = result + w[i]*(lognormalPDF(xm+dx, Nst_local)
                + lognormalPDF(xm-dx, Nst_local));
    }
    result = xr*result;

    return result;
}

//Returns the P_N(Nst)*dNst
//Ref) S.H.Kim and Kang Y. Huh, Combust Flame, vol.138, pp.336-352, 2004
//     Equation 16
inline void calculate_PdNst(scalar NumofNst, scalarField& Nst, scalar Nst_local, scalarField& result)
{
    const double alarge = 1000.0;
    scalar a,b;

    if(Nst_local > 0.5*Nst[0])
    {
        a = 0.0;
        b = 0.5*(Nst[0]+Nst[1]);
        result[0] = qgaus(a, b, Nst_local);

        for(label i=1 ; i<(NumofNst-1) ; i++)
        {
            a = 0.5*(Nst[i]+Nst[i-1]);
            b = 0.5*(Nst[i]+Nst[i+1]);
            result[i] = qgaus(a, b, Nst_local);
        }

        a = 0.5*(Nst[NumofNst-1]+Nst[NumofNst-2]);
        b = alarge;
        result[NumofNst-1] = qgaus(a, b, Nst_local);
    }
    else
    {
        result[0] = 1.0;
        for(label i=1 ; i<NumofNst ; i++)
        {
            result[i] = 0.0;
        }
    }
}

inline void Interpolation(scalarField& eta, scalarField& eta_temp, scalar NumofNst, scalar Nspecies, scalarField& Y_SLFM_temp, scalarField& result)
{
    scalar eta_size(eta.size()), eta_temp_size(eta_temp.size());
    scalar lower_eta(0), higher_eta(0);
    scalar lindex(0), hindex(0);
    scalar frac(0), diff(0);

    for(label j=0; j<eta_size ; j++)
    {
        for(label jj=0; jj<eta_temp_size ; jj++)
        {
            if(eta_temp[jj] <= eta[j])
            {
                lower_eta = eta_temp[jj];
                lindex = jj;
            }
        }
        for(label jj=eta_temp_size-1 ; jj>=0 ; jj--)
        {
            if(eta_temp[jj] > eta[j])
            {
                higher_eta = eta_temp[jj];
                hindex = jj;
            }
        }

        frac = (eta[j]-lower_eta)/(higher_eta-lower_eta);

        for(label m=0 ; m<NumofNst ; m++)
        {
            for(label i=0 ; i<Nspecies ; i++)
            {
                diff = Y_SLFM_temp[(m*eta_temp_size+hindex)*Nspecies+i]
                    - Y_SLFM_temp[(m*eta_temp_size+lindex)*Nspecies+i];

                result[(m*eta_size+j)*Nspecies+i] = Y_SLFM_temp[(m*eta_temp_size+lindex)*Nspecies+i] + frac*diff;
            }
        }
    }
}

inline void select_ij(label ij, scalar& i, scalar& j)
{
    if(ij == 0) //variance for H, O2, CO, OH, CH4, T
    {
        i = 0;
        j = 0;
    }
    else if(ij == 1)
    {
        i = 1;
        j = 1;
    }
    else if(ij == 2)
    {
        i = 2;
        j = 2;
    }
    else if(ij == 3)
    {
        i = 3;
        j = 3;
    }
    else if(ij == 4)
    {
        i = 4;
        j = 4;
    }
    else if(ij == 5)
    {
        i = 5;
        j = 5;
    }
    else if(ij == 6) //covariance between H and O2
    {
        i = 0;
        j = 1;
    }
    else if(ij == 7) //covariance between CO and OH
    {
        i = 2;
        j = 3;
    }
    else if(ij == 8) //covariance between H and CH4
    {
        i = 0;
        j = 4;
    }
    else if(ij == 9) //covariance between SPECIES and T
    {
        i = 0;
        j = 5;
    }
    else if(ij == 10)
    {
        i = 1;
        j = 5;
    }
    else if(ij == 11)
    {
        i = 2;
        j = 5;
    }
    else if(ij == 12)
    {
        i = 3;
        j = 5;
    }
    else if(ij == 13)
    {
        i = 4;
        j = 5;
    }
    else
    {
        abort();
    }
}

inline void select_species(label i, word& speciei)
{
    if(i == 0)
    {
        speciei = "H";
    }
    else if(i == 1)
    {
        speciei = "O2";
    }
    else if(i == 2)
    {
        speciei = "CO";
    }
    else if(i == 3)
    {
        speciei = "OH";
    }
    else if(i == 4)
    {
        speciei = "CH4";
    }
}
inline scalar select_CiCj(scalar i)
{
    if(i == 1 || i == 2 || i == 4 || i == 5)
    {
        return 1.0;
    }
    else if(i == 0 || i == 3)
    {
        return 2.0;
    }
    else
    {
        abort();
    }
}