//CMC P(eta) equation 
// Return probability by Beta function from sampling mixture fraction, eta (NOT mf and mfVar !!)
// eta = (1/etamax)*etacount

{
        scalar betaftn=0.0, gammaCMC=0.0, alphaCMC=0.0, betaCMC=0.0, PTOTAL=0.0; 
        scalar loop_counter = 0, mf_diff=0.0, mfnew=0.0;

        forAll(mf, celli) //cell number loop
        {
		mf_temp[celli] = mf[celli]; //karam, test

Re_calculate_pdf:;

                //Adjust mixture fraction and its variance
                if(mfVar[celli] < 1.0e-30)
                {
                        mfVar[celli] = 1.0e-30;
                }

                if(mf_temp[celli] < 0.0)
                {
                        mf_temp[celli] = 0.0;
                }

                if(mf_temp[celli] < 5.0e-4) //karam, 5e-4 -> 1e-4 )
                {
			pdf = 0.0;
                        pdf[0] = 1.0 / deta[0];//1;
			deltaftn[celli] = 0;
                        goto Peta_NextCell;
                }

                gammaCMC = mf_temp[celli]*(1.0-mf_temp[celli]) / mfVar[celli] - 1.0;

                if(gammaCMC <= 0.0)
                {
                //      Info<<"GammaCMC is lower than zero -- two delta function";
			pdf = 0.0;
                        pdf[0] = 0.5 * 1.0/deta[0]; //0.5;
                        pdf[etamax] = 0.5 * 1.0/deta[etamax]; //0.5;
			deltaftn[celli] = 1;
                        goto Peta_NextCell;
                }

                //calculate alpha and beta for PDF function
                betaCMC = (1.0-mf_temp[celli]) * gammaCMC;
                alphaCMC = mf_temp[celli] * gammaCMC;

                if(alphaCMC <= 1.0 && betaCMC <= 1.0)
                {
                //      Info<<"AlphaCMC and betaCMC are lower than one -- two delta function";
			pdf = 0.0;
                        pdf[0] = 0.5 * 1.0/deta[0]; //0.5;
                        pdf[etamax] = 0.5 * 1.0/deta[etamax]; //0.5;
			deltaftn[celli] = 1;
                        goto Peta_NextCell;
                }

                else if(alphaCMC < 1.0 && betaCMC >= 1.0)
                {
			pdf = 0.0;
                        pdf[0] = 1.0/deta[0]; //1;
			deltaftn[celli] = 0;
                        goto Peta_NextCell;
                }

                else if(alphaCMC >= 1.0 && betaCMC < 1.0)
                {
			pdf = 0.0;
                        pdf[etamax] = 1.0/deta[etamax]; //1;
			deltaftn[celli] = 0;
                        goto Peta_NextCell;
                }

                betaftn = std::exp(  gammaln(alphaCMC) + gammaln(betaCMC) - gammaln(alphaCMC + betaCMC) );       //Numerial recipe 2nd Edition 206page  

                if(betaftn < 1.0e-50)
                {
set_onedelta:;
			pdf = 0.0;
                        for(label etacount = 0; etacount<=etamax ; etacount++) //set one delta function
                        {
                                if( etaValue[etacount] >= mf_temp[celli])
                                {
                                        pdf[etacount] = 1.0/deta[etacount]; //1;
					deltaftn[celli] = 0;
                                        goto Peta_NextCell;
                                }
                        }

                }

                for(label etacount=0 ; etacount <= etamax ; etacount++)
                {
	                if(etacount == 0 || etacount == etamax)
                        {
         	             	pdf[etacount] = 0.0;
                        }
                        else
                        {
                 	     	pdf[etacount] = Foam::pow( etaValue[etacount] , alphaCMC - 1.0 ) * Foam::pow ( 1.0 - etaValue[etacount] ,  betaCMC  -1.0 ) / betaftn ;
                        }
                }
		deltaftn[celli] = 2;

                PTOTAL = 0.0; //PTOTAL initialization in each cell
		
		f = 1.0;

		PTOTAL = integration(deltaftn[celli], MFcut, Neta, pdf, f);

		for(label etacount=0; etacount<=etamax ; etacount++)
		{
			pdf[etacount] = pdf[etacount]/PTOTAL;
		}

		mfnew = integration(deltaftn[celli], MFcut, Neta, pdf, etaValue);

                mf_diff = (mf[celli] - mfnew)/mf[celli];

                if(std::fabs(mf_diff) > 0.00001)//1e-5)
                {
                        mf_temp[celli] = mf_temp[celli] + 0.5*mf_diff*mf[celli];
                        if(loop_counter > 100)
                        {
                                Info<<"Force set pdf to one delta"<<endl;
                                goto set_onedelta;
                        }
                        loop_counter++;
                        goto Re_calculate_pdf;
                }

Peta_NextCell:; 

		for(label etacount = 0; etacount <= etamax ; etacount++)
		{
			Peta[etacount][celli] = pdf[etacount];
		}
        	loop_counter = 0;
        }
}