class BetaPDF
{

    //private variable 

    //beta-pdf parameter alpha and beta
    scalar alpha_, beta_;

    //cutting point and number of eta-space
    scalarField etaCut_, N_;

    //detailed integration space
    scalarField etaSpace_;

    //detailed integration space, part
    typedef List<scalarField> scalarFieldArray1d;
    scalarFieldArray1d etaPart_;

    //AMC, exp(-2*(erf^-1(2*eta - 1))^2) field
    //for detailed integration space
    scalarField AMCfine_;

    //flag to check forced-delta ftn 
    //or delta ftn at oxidizer or fuel
    bool fdelta_ = false, delta_ox = false, delta_fu = false;

public:

    // Constructor
    BetaPDF(IOdictionary& SLFMdict)
    :
    alpha_(0),
    beta_(0),
    etaCut_(SLFMdict.lookup("detailedEta")),
    N_(SLFMdict.lookup("detailedN")),
    etaPart_(N_.size())
    {
        //Ref. F.Liu et al., INT. J. THERM. SCI. 41 (2002) 763-772.
        Info<<"Construct Beta-PDF"<<endl;

        scalar del(0);
        label cnt(0);

        etaSpace_.append(etaCut_[0]);
        for(label i=0 ; i<etaCut_.size()-1 ; i++)
        {
            scalarField tmp;

            tmp.append(etaCut_[i]);

            del = (etaCut_[i+1] - etaCut_[i])/N_[i];
            for(label j=0 ; j<N_[i] ; j++)
            {
                etaSpace_.append(etaSpace_[cnt] + del);
                tmp.append(etaSpace_[cnt]+del);
                cnt++;
            }
            etaPart_[i].append(tmp);
        }

        AMCfine_ = AMC(etaSpace_);
    }

    // Destructor
    ~BetaPDF(){}

    // Member functions

    // set beta-pdf parameter alpha_ and beta_
    void setParameter(scalar& mf, scalar& mfVar)
    {
        fdelta_ = false;
        delta_ox = false;
        delta_fu = false;

        scalar gamma = mf*(1.0-mf)/(mfVar+SMALL) - 1.0;

        if(gamma<SMALL || (mfVar)/(mf*(1-mf)+SMALL)<0.001)
        {
            fdelta_ = true;
        }
        else
        {
            alpha_ = max(mf*gamma, 0.0);
            beta_ = max((1.0-mf)*gamma, 0.0);

            if(alpha_ < 1.0)
            {
                delta_ox = true;
            }
            if(beta_ < 1.0)
            {
                delta_fu = true;
            }
        }
    }
    // beta-pdf weighted integration for given mf, mfVar and f
    scalar evaluate(scalar& mf, scalar& mfVar, scalarField& etaValue, scalarField& f)
    {
        scalar result(0);

        setParameter(mf, mfVar);

        if(fdelta_ == false)
        {
            limit_ab();

            scalar num(0), den(0);
            for(label i=0 ; i<etaCut_.size()-1 ; i++)
            {
                scalarField fPart = interpolateXY(etaPart_[i], etaValue, f);

                scalarField fx = fPart * etaFunc(alpha_, beta_, etaPart_[i]);
                scalarField gx = etaFunc(alpha_, beta_, etaPart_[i]);

                num += simps(etaCut_[i], etaCut_[i+1], N_[i], fx);
                den += simps(etaCut_[i], etaCut_[i+1], N_[i], gx);
            }
            num += f[0]*Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL);
            den += Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL);
            num += f[etaValue.size()-1]*Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
            den += Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);

            result = num/den;
        }
        else
        {
            result = interpolateXY(mf, etaValue, f);
        }

        return result;
    }
    scalar value(scalar& mf, scalar& mfVar, scalar& etaValue)
    {
        scalar result(0);

        setParameter(mf, mfVar);
        
        if(fdelta_ == false)
        {
            limit_ab();

            scalar den(0);
            for(label i=0 ; i<etaCut_.size()-1 ; i++)
            {
                scalarField gx = etaFunc(alpha_, beta_, etaPart_[i]);

                den += simps(etaCut_[i], etaCut_[i+1], N_[i], gx);
            }
            den += Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL);
            den += Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);

            scalarField pdf = etaFunc(alpha_, beta_, etaSpace_)/den;
            result = interpolateXY(etaValue, etaSpace_, pdf);
        }

        return result;
    }
    void limit_ab()
    {
        scalar fmax = 1.0+(beta_-1.0)/(alpha_-1.0);
        fmax = 1.0/fmax;

        if(alpha_ > 500.0)
        {
            alpha_ = 500.0;
            beta_ = (alpha_-1.0-fmax*(alpha_-2.0))/fmax;
        }
        else if(beta_ > 500.0)
        {
            beta_ = 500.0;
            alpha_ = (1.0+fmax*(beta_-2.0))/(1.0-fmax);
        }
    }
    scalarField etaFunc(scalar& a, scalar& b, scalarField& eta)
    {
        return pow(eta, a-1.0)*pow(1.0-eta, b-1.0);
    }
    scalar etaFunc(scalar& a, scalar& b, scalar& eta)
    {
        return Foam::pow(eta, a-1.0)*Foam::pow(1.0-eta, b-1.0);
    }
    //extended Simpson's rule (Numerical recipes, 2nd Ed. p.128)
    //for equally spaced and even N intervals (or odd N+1 points) 
    scalar simps(scalar& xl, scalar& xh, scalar& N, scalarField& fx)
    {
        scalar evensum(0.0), oddsum(0.0), sum(0.0);

        scalar h = (xh - xl)/N;

        for(label i=0 ; i<fx.size() ; i++)
        {
            if(i%2 == 0)
            {
                evensum += fx[i];
            }
            else
            {
                oddsum += fx[i];
            }
        }
 
        sum = -1.0*fx.first() + 2.0*evensum + 4.0*oddsum - 1.0*fx.last();

        return sum*(h/3.0);
    }
    //Amplitude Mapping Closure (from KIVA)
    //Define exp(-2*(erf^-1(2*eta - 1))^2)
    scalarField AMC(scalarField& eta)
    {
        const scalar pi = 3.141592;
        const scalar spi = Foam::sqrt(pi);

        scalarField result(eta.size(), 0.0);

        for(label i=0 ; i<eta.size() ; i++)
        {
            scalar a1=0.5, a0=(2.0*eta[i]-1.0), slope, da=GREAT;
            while(mag(da) > SMALL)
            {
                slope = 2.0/spi*Foam::exp(-1.0*Foam::pow(a1,2.0));
                da = (a0 - Foam::erf(a1))/slope;
                a1 = a1+da;
            }

            result[i] = Foam::exp(-2.0*Foam::pow(a1,2.0));
        }

        return result;
    }
    scalar AMC(scalar& eta)
    {
        return interpolateXY(eta, etaSpace_, AMCfine_);
    }
    void C1coeff(scalar& mf, scalarField& varValue, scalarField& C1table)
    {
        scalar maxVar = mf*(1.0-mf);
        scalarField x, fx;

        x.append(0.0);
        x.append(etaSpace_);
        x.append(1.0);

        fx.append(0.0);
        fx.append(AMCfine_);
        fx.append(0.0);

        for(label v=1 ; v<varValue.size()-1 ; v++)
        {
            scalar var = maxVar*varValue[v];
            C1table[v] = 1.0/evaluate(mf,var,x,fx);
        }
    }
};