//get conditional SDR
        scalar C2(0); //After 22.Oct.2012
        scalar Airside(0);
        scalar Fuelside(0);

        forAll(rho , celli)  //Before 22.Oct.2012
        {
                Airside = Peta[0][celli]*deta[0]; 
                Fuelside = Peta[etamax][celli]*deta[etamax]; 
                if( (1-Airside) < 0.000001 || (1-Fuelside) < 0.000001 )//if there exist delta function at eta = 0 or 1
                {
                        for(label j = 0 ; j<= etamax ; j++)
                        {
                                CSDRCELL[j][celli] = 0;
                        }
                }
                else
                {
                        scalar C1(0);
			for(label j=0 ; j<=etamax ;j++)
			{
				f[j] = AMC(etaValue[j]);
				pdf[j] = Peta[j][celli];
			}
			C1 = integration(deltaftn[celli], MFcut, Neta, pdf, f);
                        //for(label j = 0 ; j <= etamax ; j++)
                        //{
                        //        scalar eta = etaValue[j];
                        //        C1 += AMC(eta) * Peta[j][celli] * deta[j];
                        //}

                        C2 = SDR[celli] / C1; //Eqn (13) from SH Kim C&F paper 120:75-90(20)

                        for(label j = 0 ; j <= etamax ; j++)
                        {
                                scalar eta = etaValue[j];
                                CSDRCELL[j][celli] = C2 * AMC(eta);
                        }
                }
		if(F_total[celli] > 0)
               	{
			for(label k=0 ; k<=group ; k++)
			{
                		Ffrac[k][celli] = F[k][celli]/F_total[celli];
			}
              	}
		else if(F_total[celli] <= 0)
		{
                        for(label k=0 ; k<=group ; k++)
                        {
                                Ffrac[k][celli] = 0.0;
                        }			
		}
        }

//get density weighted averaged value of conditional SDR for each eta and flame-group 
        for(label k=0 ; k <= group ; k++)
        {
                for(label j = 0 ; j<=etamax ; j++) //density weighted average of conditional SDR 
                {
                        scalar SDRTemp(0), DENOMTemp(0);

                        SDRTemp = fvc::domainIntegrate(Ffrac[k]*rho*Peta[j]*CSDRCELL[j]).value();
                        DENOMTemp = fvc::domainIntegrate(Ffrac[k]*rho*Peta[j]).value();
			
                        if(SDRTemp == 0 || DENOMTemp == 0)
                        {
                                SDRCMC[j+k*(etamax+1)] = 0;
                        }
                        else
                        {
                                SDRCMC[j+k*(etamax+1)] = SDRTemp / DENOMTemp ;
                        }
//Info<<"eta space"<<tab<<j<<tab<<etaValue[j]<<tab<<SDRCMC[j+k*(etamax+1)]<<endl;
                }
        }
/*
forAll(mf, celli)
{
	if(Ffrac[1][celli] > 0.05)
	{
		vector axialV(0, 1, 0);
		scalar axialC = mesh.C()[celli] & axialV;
		if(axialC < 0.15 && axialC > 0.05)
		{
			Info<<celli<<tab<<mesh.C()[celli]<<endl;
			for(label j = 0 ; j<=etamax ; j++)
			{
				Info<<etaValue[j]<<tab<<CSDRCELL[j][celli]<<tab<<Ffrac[1][celli]*rho[celli]*mesh.V()[celli]<<tab<<Peta[j][celli]<<endl;
			}
		}
		if(axialC > 0.3 && axialC < 0.35)
                {
                        Info<<celli<<tab<<mesh.C()[celli]<<endl;
                        for(label j = 0 ; j<=etamax ; j++)
                        {
                                Info<<etaValue[j]<<tab<<CSDRCELL[j][celli]<<tab<<Ffrac[1][celli]*rho[celli]*mesh.V()[celli]<<tab<<Peta[j][celli]<<endl;
                        }
                }
	}
}
*/
