class BetaPDF
{

    //private variable 

    //raw parameter mixture fraction and mf variance
    scalar mf_, mfVar_;

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

    //cutting points
    scalarField etaCut_;

    //number of eta-space
    labelList N_;

    //detailed integration space
    scalarField etaSpace_;

    //pdf numerators
    scalarField pdfNum_;

    //pdf denominator
    scalar pdfDen_;

    //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.
        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_);

        pdfNum_ = scalarField(etaSpace_.size(), 0.0);
        pdfDen_ = 1.0;
    }

    // Destructor
    ~BetaPDF(){}

    // Member functions

    // set beta-pdf parameter alpha_ and beta_
    void setParameter(const scalar mf, const scalar mfVar)
    {
        mf_ = mf;
        mfVar_ = 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;
            }

            limit_ab();
        }

        pdfNum_ = etaFunc(alpha_, beta_, etaSpace_);
        pdfDen_ = integrate(pdfNum_)
            + Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL)
            + Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
    }
    // beta-pdf weighted integration for given mf, mfVar and f
    scalar evaluate(const scalar mf, const scalar mfVar, const scalarField& etaValue, const scalarField& f)
    {
        scalar result(0);

        setParameter(mf, mfVar);

        if(fdelta_ == false)
        {
            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;
    }

    // beta-pdf weighted integration for given mf, mfVar and f
    scalar betaIntegrate(const scalarField& f) const
    {
        scalar result = 0.0;

        if(fdelta_)
        {
            result = interpolateXY(mf_, etaSpace_, f);
        }
        else
        {
            scalar num = integrate (f * pdfNum_)
                + f.first()*Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL)
                + f.last()*Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
            result = num/pdfDen_;
        }

        return result;
    }

    scalar integrate(const scalarField& f) const
    {
        scalar total = 0.0;

        forAll(N_, i)
        {
            const labelList::subList prev(N_, i);
            const label baseI = sum(prev);

            const scalarField::subField subInterval(f, N_[i]+1, baseI);

            total += simps(etaCut_[i], etaCut_[i+1], N_[i], subInterval);
        }

        return total;
    }

    scalar value(const scalar mf, const scalar mfVar, const scalar etaValue)
    {
        scalar result(0);

        setParameter(mf, mfVar);
        
        if(fdelta_ == false)
        {
            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(const scalar a, const scalar b, const scalarField& eta) const
    {
        return pow(eta, a-1.0)*pow(1.0-eta, b-1.0);
    }
    scalar etaFunc(const scalar a, const scalar b, const scalar eta) const
    {
        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(const scalar xl, const scalar xh, const label N, const UList<scalar>& fx) const
    {
        scalar evensum(0.0), oddsum(0.0), sum(0.0);

        scalar h = (xh - xl)/scalar(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(const scalarField& eta) const
    {
        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(const scalar eta) const
    {
        return interpolateXY(eta, etaSpace_, AMCfine_);
    }
    void C1coeff(const scalar mf, const 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);
        }
    }
};