#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();
	}
}	
