diff --git a/solvers_post/SLFMFoam/BetaPDF.H b/solvers_post/SLFMFoam/BetaPDF.H
index d23dd56..0879187 100644
--- a/solvers_post/SLFMFoam/BetaPDF.H
+++ b/solvers_post/SLFMFoam/BetaPDF.H
@@ -3,15 +3,27 @@ class BetaPDF
 
     //private variable 
 
+    //raw parameter mixture fraction and mf variance
+    scalar mf_, mfVar_;
+
     //beta-pdf parameter alpha and beta
     scalar alpha_, beta_;
 
-    //cutting point and number of eta-space
-    scalarField etaCut_, N_;
+    //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_;
@@ -59,6 +71,9 @@ public:
         }
 
         AMCfine_ = AMC(etaSpace_);
+
+        pdfNum_ = scalarField(etaSpace_.size(), 0.0);
+        pdfDen_ = 1.0;
     }
 
     // Destructor
@@ -67,8 +82,11 @@ public:
     // Member functions
 
     // set beta-pdf parameter alpha_ and beta_
-    void setParameter(scalar& mf, scalar& mfVar)
+    void setParameter(const scalar mf, const scalar mfVar)
     {
+        mf_ = mf;
+        mfVar_ = mfVar;
+
         fdelta_ = false;
         delta_ox = false;
         delta_fu = false;
@@ -92,10 +110,17 @@ public:
             {
                 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(scalar& mf, scalar& mfVar, scalarField& etaValue, scalarField& f)
+    scalar evaluate(scalar mf, scalar mfVar, scalarField& etaValue, scalarField& f)
     {
         scalar result(0);
 
@@ -103,8 +128,6 @@ public:
 
         if(fdelta_ == false)
         {
-            limit_ab();
-
             scalar num(0), den(0);
             for(label i=0 ; i<etaCut_.size()-1 ; i++)
             {
@@ -130,7 +153,45 @@ public:
 
         return result;
     }
-    scalar value(scalar& mf, scalar& mfVar, scalar& etaValue)
+
+    // beta-pdf weighted integration for given mf, mfVar and f
+    scalar betaIntegrate(const scalarField& f)
+    {
+        scalar result(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)
+    {
+        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(scalar mf, scalar mfVar, scalar etaValue)
     {
         scalar result(0);
 
@@ -138,8 +199,6 @@ public:
         
         if(fdelta_ == false)
         {
-            limit_ab();
-
             scalar den(0);
             for(label i=0 ; i<etaCut_.size()-1 ; i++)
             {
@@ -172,21 +231,21 @@ public:
             alpha_ = (1.0+fmax*(beta_-2.0))/(1.0-fmax);
         }
     }
-    scalarField etaFunc(scalar& a, scalar& b, scalarField& eta)
+    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)
+    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 simps(const scalar xl, const scalar xh, const label N, const UList<scalar>& fx)
     {
         scalar evensum(0.0), oddsum(0.0), sum(0.0);
 
-        scalar h = (xh - xl)/N;
+        scalar h = (xh - xl)/scalar(N);
 
         for(label i=0 ; i<fx.size() ; i++)
         {
@@ -228,11 +287,11 @@ public:
 
         return result;
     }
-    scalar AMC(scalar& eta)
+    scalar AMC(scalar eta)
     {
         return interpolateXY(eta, etaSpace_, AMCfine_);
     }
-    void C1coeff(scalar& mf, scalarField& varValue, scalarField& C1table)
+    void C1coeff(scalar mf, scalarField& varValue, scalarField& C1table)
     {
         scalar maxVar = mf*(1.0-mf);
         scalarField x, fx;
diff --git a/solvers_post/SLFMFoam/preIntegration.H b/solvers_post/SLFMFoam/preIntegration.H
index 8930a49..a5c4b05 100644
--- a/solvers_post/SLFMFoam/preIntegration.H
+++ b/solvers_post/SLFMFoam/preIntegration.H
@@ -47,18 +47,21 @@ for(label j=1 ; j<etamax-1; j++)
 Info<<"Construct Temp and Rgas table"<<endl;
 for(label j=0 ; j<etamax ; j++)
 {
-    scalar mf = etaValue[j];
+    const scalar mf = etaValue[j];
     for(label v=0 ; v<=NVar ; v++)
     {
-        scalar mfVar = varValue[v]*mf*(1.0-mf);
+        const scalar mfVar = varValue[v]*mf*(1.0-mf);
+
+        bpdf.setParameter(mf, mfVar);
+
         for(label n=0 ; n<NstList.size() ; n++)
         {
-            Ttable[j][v][n] = bpdf.evaluate(mf, mfVar, etaValue, T_SLFM[n]);
-            Rtable[j][v][n] = bpdf.evaluate(mf, mfVar, etaValue, Rgas_SLFM[n]);
+            Ttable[j][v][n] = bpdf.betaIntegrate(T_SLFM[n]);
+            Rtable[j][v][n] = bpdf.betaIntegrate(Rgas_SLFM[n]);
             forAll(postSpecieIndices, yi)
             {
                 const label y = postSpecieIndices[yi];
-                Ytable[yi][j][v][n] = bpdf.evaluate(mf, mfVar, etaValue, Y_SLFM[y][n]);
+                Ytable[yi][j][v][n] = bpdf.betaIntegrate(Y_SLFM[y][n]);
             }
         }
     }
diff --git a/solvers_post/SLFMFoam/readSLFMlib.H b/solvers_post/SLFMFoam/readSLFMlib.H
index fd5a706..78f9acc 100644
--- a/solvers_post/SLFMFoam/readSLFMlib.H
+++ b/solvers_post/SLFMFoam/readSLFMlib.H
@@ -12,7 +12,7 @@ typedef List<scalarField> scalarFieldArray1d;
 typedef List<scalarFieldArray1d> scalarFieldArray2d;
 typedef List<scalarFieldArray2d> scalarFieldArray3d;
 
-CMC::FlameStructure slfmLibrary(SLFMdict, etaValue.size(), Y.size());
+CMC::FlameStructure slfmLibrary(SLFMdict, etaGrid.size(), Y.size());
 
 scalarFieldArray2d& Y_SLFM (slfmLibrary.Y());
 scalarFieldArray2d& W_SLFM (slfmLibrary.W());
@@ -88,7 +88,7 @@ for(label n=0 ; n<NstList.size() ; n++)
             Y_temp[j] = max(0.0, Yj);
         }
 
-        Y_SLFM[spI][n] = interpolateXY(etaValue, etaValue_SLFM, Y_temp);
+        Y_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp);
 
         //read reaction rate (from SLFM library)
         for(label j=0 ; j<etamax_SLFM ; j++)
@@ -96,7 +96,7 @@ for(label n=0 ; n<NstList.size() ; n++)
             fin>>W_temp[j];
         }
 
-        W_SLFM[spI][n] = interpolateXY(etaValue, etaValue_SLFM, W_temp);
+        W_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp);
 
         //blank line
         fin.getLine(gbg);
@@ -115,7 +115,7 @@ for(label n=0 ; n<NstList.size() ; n++)
             Y_temp[j] = max(0.0, Yj);
         }
 
-        T_SLFM[n] = interpolateXY(etaValue, etaValue_SLFM, Y_temp);
+        T_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp);
 
         //read heat source
         for(label j=0 ; j<etamax_SLFM ; j++)
@@ -123,7 +123,7 @@ for(label n=0 ; n<NstList.size() ; n++)
             fin>>W_temp[j];
         }
 
-        Q_SLFM[n] = interpolateXY(etaValue, etaValue_SLFM, W_temp);
+        Q_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp);
 
         //blank line
         fin.getLine(gbg);
@@ -140,7 +140,7 @@ for(label n=0 ; n<NstList.size() ; n++)
     }
     //Info<<rho_temp<<endl;
 
-    rho_SLFM[n] = interpolateXY(etaValue, etaValue_SLFM, rho_temp);
+    rho_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, rho_temp);
 
     //blank line
     fin.getLine(gbg);
@@ -155,7 +155,7 @@ for(label n=0 ; n<NstList.size() ; n++)
         fin>>h_temp[j];
     }
 
-    h_SLFM[n] = interpolateXY(etaValue, etaValue_SLFM, h_temp);
+    h_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, h_temp);
 
     for(label i=0 ; i<Y.size() ; i++)
     {