#ifndef GSL_DEF
#define GSL_DEF
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#endif
#include "residue.h"
#include "macros.h"
#include <cmath>
#include <stdio.h>
#include "timing.hpp"

void REGRID(double* ydata,double* ydotdata,UserData data){
    size_t nPts = data->npts;
    size_t nvar = data->nvar;
    double WORK[nPts],XNEW[nPts],XOLD[nPts]; 
    double P,R,R0,R1,R2,TV1,TV2,XLEN,B1,B2;
    double DX,DDX;
    double DETA,ETAJ,DEL;
    size_t ISTART; 
    P = data->PCAD;
    R = data->RGTC;
    R0 =1.0 - P;
    R1 =P*R/(R+1.0);
    R2 =P - R1 ; 
    TV1 = 0.0; 
    TV2 = 0.0;

    for(size_t i=2;i<=nPts;i++){
        TV1 = TV1 + fabs(T(i)-T(i-1));
    }
    
    for(size_t i=2;i<=nPts-1;i++){
        TV2 = TV2 +  fabs( (T(i+1)-T(i))/(R(i+1)-R(i)) - (T(i)-T(i-1))/(R(i)-R(i-1)) );
    }
    
    XLEN = R(nPts) - R(1); 
    B1 = R1* XLEN/ ( (1.0-P)*TV1);
    B2 = R2* XLEN/ ( (1.0-P)*TV2);

    //Compute partial sum of weight function
    WORK(1) =0.0;
    DX = 0.0;
    for(size_t i=2;i<=nPts-1;i++){
        DX = R(i) - R(i-1) ; 
        WORK(i) = DX + B1*fabs(T(i)-T(i-1)) + WORK(i-1)+B2*fabs( (T(i+1)-T(i))/(R(i+1)-R(i)) - (T(i)-T(i-1))/(R(i)-R(i-1)) );
    }
    
    DDX = R(nPts) - R(nPts-1); 
    WORK(nPts) = WORK(nPts-1) +DDX ; 
    
    for(size_t i=2;i<=nPts;i++){
        WORK(i) = WORK(i)/WORK(nPts); 
    }
        
    XNEW(1)=R(1); 
    XNEW(nPts)=R(nPts); 
    ISTART =2;
    //double DETA;
    DETA = 1.0/(nPts-1) ; 

    //Fill new grid XNEW[nPts],duplicate from PREMIX REGRID SUBROUTINE
    for(size_t j=2;j<=nPts-1;j++){
        ETAJ= (j-1)*DETA; 

        for(size_t i=ISTART;i<=nPts;i++){
            if(ETAJ <= WORK(i)){
                DEL = (ETAJ-WORK(i-1))/(WORK(i)-WORK(i-1)) ; 
                XNEW(j)=R(i-1)+(R(i)-R(i-1))*DEL;
                break;
            }else{
                ISTART =i;
            }

        }
    }

    //Interpolate solution based on XNEW[]&XOLD[]
    for(size_t i=1;i<=nPts;i++){
        XOLD[i-1]=R(i);
    }

    INTERPO(ydata,ydotdata,nvar,nPts,XNEW,XOLD);
}


void INTERPO(double* y,double* ydot,const size_t nvar,size_t nPts,const double XNEW[], const double XOLD[]){
    double ytemp[nPts],ydottemp[nPts] ; 
    
    gsl_interp_accel* acc; 
    gsl_spline* spline; 
    acc = gsl_interp_accel_alloc();
    spline = gsl_spline_alloc(gsl_interp_steffen,nPts);
    
    gsl_interp_accel* accdot; 
    gsl_spline* splinedot; 
    accdot = gsl_interp_accel_alloc();
    splinedot = gsl_spline_alloc(gsl_interp_steffen,nPts);

    for(size_t j=0;j<nvar;j++){
        
        for(size_t i=0;i<nPts;i++){
            ytemp[i]=y[j+i*nvar];
            ydottemp[i]=y[j+i*nvar];
        }

        gsl_spline_init(spline,XOLD,ytemp,nPts);
        gsl_spline_init(splinedot,XOLD,ydottemp,nPts);
        
        for(size_t i=0;i<nPts;i++){
            y[j+i*nvar]=gsl_spline_eval(spline,XNEW[i],acc);
            ydot[j+i*nvar]=gsl_spline_eval(splinedot,XNEW[i],accdot);
        }

    }

    gsl_interp_accel_free(acc);
    gsl_spline_free(spline);

    gsl_interp_accel_free(accdot);
    gsl_spline_free(splinedot);

}

double maxTemperaturePosition(const double* y,const size_t nt,const size_t nvar,const double *x ,size_t nPts){
    double maxT = 0.0e0;
    double TempT = 0.0e0;
    double pos = 0.0e0;
    for (size_t i=1;i<=nPts;i++){
        TempT = y[(i-1)*nvar+nt] ; 
        if(TempT > maxT){
            maxT = TempT ; 
            pos = x[i-1];
        }
    }
    return(pos); 
}

double maxTemperature(const double* y,const size_t nt,const size_t nvar ,size_t nPts){
    double maxT = 0.0e0;
    double TempT = 0.0e0;
    //int index = 0 ; 
    for (size_t i=1;i<=nPts;i++){
        TempT = y[(i-1)*nvar+nt] ; 
        if(TempT > maxT){
            maxT = TempT ; 
        }
    }
    return(maxT); 
}
//Index here is 1-based
int maxTemperatureIndex(const double* y,const size_t nt,const size_t nvar ,size_t nPts){
    double maxT = 0.0e0;
    double TempT = 0.0e0;
    int index = 0 ; 
    for (size_t i=1;i<=nPts;i++){
        TempT = y[(i-1)*nvar+nt] ; 
        if(TempT > maxT){
            maxT = TempT ; 
            index = i; 
        }
    }
    return(index); 
}

double maxCurvPositionR(const double* y, const size_t nt, 
	               const size_t nvar, const size_t nr, size_t nPts){
	double maxCurvT=0.0e0;
	double gradTp=0.0e0;
	double gradTm=0.0e0;
	double curvT=0.0e0;
	double dx=0.0e0;
    double dr=0.0e0;
	double pos=0.0e0;
	size_t j,jm,jp;
    size_t r,rp,rm;
	for (size_t i = 1; i <nPts-1; i++) {
		j=i*nvar+nt;
		jm=(i-1)*nvar+nt;
		jp=(i+1)*nvar+nt;
        r = i*nvar+nr;
        rm = (i-1)*nvar+nr;
        rp = (i+1)*nvar+nr; 
		//gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
		//gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
		//dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
		//curvT=(gradTp-gradTm)/dx;
        gradTp = fabs((y[jp]-y[j])/(y[rp]-y[r]));
        gradTp = fabs((y[j]-y[jm])/(y[r]-y[rm]));
		dr = 0.5e0*(y[rp]-y[rm]);
        curvT=(gradTp-gradTm)/dr; 
        if (curvT>=maxCurvT) {
			maxCurvT=curvT;
			pos=y[r];
		}
	}
	return(pos);
}

int maxCurvIndexR(const double* y, const size_t nt, 
	               const size_t nvar, const size_t nr, size_t nPts){
	double maxCurvT=0.0e0;
	double gradTp=0.0e0;
	double gradTm=0.0e0;
	double curvT=0.0e0;
	double dx=0.0e0;
    double dr=0.0e0;
	int pos=0;
	size_t j,jm,jp;
    size_t r,rm,rp;
	for (size_t i = 1; i <nPts-1; i++) {
		j=i*nvar+nt;
		jm=(i-1)*nvar+nt;
		jp=(i+1)*nvar+nt;
        r = i*nvar+nr;
        rm = (i-1)*nvar+nr;
        rp = (i+1)*nvar+nr; 
		//gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
		//gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
		//dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
		//curvT=(gradTp-gradTm)/dx;
        gradTp = fabs((y[jp]-y[j])/(y[rp]-y[r]));
        gradTp = fabs((y[j]-y[jm])/(y[r]-y[rm]));
		dr = 0.5e0*(y[rp]-y[rm]);
        curvT=(gradTp-gradTm)/dr; 
        if (curvT>=maxCurvT) {
			maxCurvT=curvT;
			pos=i;
		}
	}
	return(pos);
}

double maxGradPosition(const double* y, const size_t nt, 
	               const size_t nvar, const double* x, size_t nPts){
	double maxGradT=0.0e0;
	double gradT=0.0e0;
	double pos=0.0e0;
	size_t j,jm;
	for (size_t i = 1; i <nPts; i++) {
		j=i*nvar+nt;
		jm=(i-1)*nvar+nt;
		gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
		if (gradT>=maxGradT) {
			maxGradT=gradT;
			pos=x[i];
		}
	}
	return(pos);
}

int maxGradIndex(const double* y, const size_t nt, 
	         const size_t nvar, const double* x, size_t nPts){
	double maxGradT=0.0e0;
	double gradT=0.0e0;
	int pos=0;
	size_t j,jm;
	for (size_t i = 1; i <nPts; i++) {
		j=i*nvar+nt;
		jm=(i-1)*nvar+nt;
		gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
		if (gradT>=maxGradT) {
			maxGradT=gradT;
			pos=i;
		}
	}
	return(pos);
}

double maxCurvPosition(const double* y, const size_t nt, 
	               const size_t nvar, const double* x, size_t nPts){
	double maxCurvT=0.0e0;
	double gradTp=0.0e0;
	double gradTm=0.0e0;
	double curvT=0.0e0;
	double dx=0.0e0;
	double pos=0.0e0;
	size_t j,jm,jp;
	for (size_t i = 1; i <nPts-1; i++) {
		j=i*nvar+nt;
		jm=(i-1)*nvar+nt;
		jp=(i+1)*nvar+nt;
		gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
		gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
		dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
		curvT=(gradTp-gradTm)/dx;
		if (curvT>=maxCurvT) {
			maxCurvT=curvT;
			pos=x[i];
		}
	}
	return(pos);
}

int maxCurvIndex(const double* y, const size_t nt, 
	               const size_t nvar, const double* x, size_t nPts){
	double maxCurvT=0.0e0;
	double gradTp=0.0e0;
	double gradTm=0.0e0;
	double curvT=0.0e0;
	double dx=0.0e0;
	int pos=0;
	size_t j,jm,jp;
	for (size_t i = 1; i <nPts-1; i++) {
		j=i*nvar+nt;
		jm=(i-1)*nvar+nt;
		jp=(i+1)*nvar+nt;
		gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
		gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
		dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
		curvT=(gradTp-gradTm)/dx;
		if (curvT>=maxCurvT) {
			maxCurvT=curvT;
			pos=i;
		}
	}
	return(pos);
}

//double isothermPosition(const double* y, const double T, const size_t nt, 
//	                const size_t nvar, const double* x, const size_t nPts){
//	double pos=x[nPts-1];
//	size_t j;
//	for (size_t i = 1; i <nPts; i++) {
//		j=i*nvar+nt;
//		if (y[j]<=T) {
//			pos=x[i];
//			break;
//		}
//	}
//	return(pos);
//}

//******************* scan the temperature from left to right ******************//
double isothermPosition(const double* y, const double T, const size_t nt, 
	                const size_t nvar, const double* x, const size_t nPts){
	double pos=x[0];
	size_t j;
	for (size_t i = 0; i <nPts; i++) {
		j=i*nvar+nt;
		if (y[j]>=T) {
			pos=x[i];
			break;
		}
	}
	return(pos);
}

void updateSolution(double* y, double* ydot, const size_t nvar,
		    const double xOld[],const double xNew[],const size_t nPts){

	double ytemp[nPts],ydottemp[nPts];

	gsl_interp_accel* acc;
	gsl_spline* spline;
	acc = gsl_interp_accel_alloc();
	spline = gsl_spline_alloc(gsl_interp_steffen, nPts);

	gsl_interp_accel* accdot;
	gsl_spline* splinedot;
	accdot    = gsl_interp_accel_alloc();
	splinedot = gsl_spline_alloc(gsl_interp_steffen, nPts);

	for (size_t j = 0; j < nvar; j++) {

		for (size_t i = 0; i < nPts; i++) {
			ytemp[i]=y[j+i*nvar];
			ydottemp[i]=ydot[j+i*nvar];
		}

		gsl_spline_init(spline,xOld,ytemp,nPts);
		gsl_spline_init(splinedot,xOld,ydottemp,nPts);

		for (size_t i = 0; i < nPts; i++) {

		      	y[j+i*nvar]=gsl_spline_eval(spline,xNew[i],acc);
		      	ydot[j+i*nvar]=gsl_spline_eval(splinedot,xNew[i],accdot);
		}
	}

	//Exploring "fixing" boundary conditions:
	//for (size_t j = 1; j < nvar; j++) {
	//	//printf("%15.6e\t%15.6e\n", y[j],y[j+nvar]);
	//	y[j]=y[j+nvar];
	//	//y[j+(nPts-1)*nvar]=y[j+(nPts-2)*nvar];
	//	//ydot[j+nvar]=ydot[j];
	//}
	//y[0]=0.0e0;
	
	gsl_interp_accel_free(acc);
	gsl_spline_free(spline);

	gsl_interp_accel_free(accdot);
	gsl_spline_free(splinedot);
}

//Locate bath gas:
size_t BathGasIndex(UserData data){
	size_t index=0;
	double max1;
	double max=data->gas->massFraction(0);
	for(size_t k=1;k<data->nsp;k++){
		max1=data->gas->massFraction(k);
		if(max1>=max){
			max=max1;
			index=k;
		}
	}
	return(index+1);
}

//Locate Oxidizer:
size_t oxidizerIndex(UserData data){
	size_t index=0;
	for(size_t k=1;k<data->nsp;k++){
		if(data->gas->speciesName(k-1)=="O2"){
			index=k;
		}
	}
	return(index);
}

//Locate OH:
size_t OHIndex(UserData data){
	size_t index=0;
	for(size_t k=1;k<data->nsp;k++){
		if(data->gas->speciesName(k-1)=="OH"){
			index=k;
		}
	}
	return(index);
}

//Locate HO2:
size_t HO2Index(UserData data){
	size_t index=0;
	for(size_t k=1;k<data->nsp;k++){
		if(data->gas->speciesName(k-1)=="HO2"){
			index=k;
		}
	}
	return(index);
}

//Locate species index:
size_t specIndex(UserData data,char *specName){
	size_t index=0;
	for(size_t k=1;k<data->nsp;k++){
		if(data->gas->speciesName(k-1)==specName){
			index=k;
		}
	}
	return(index);
}

int setAlgebraicVariables(N_Vector* id, UserData data){
	double *iddata;
    char* specPtr[2]; 

  	N_VConst(ONE, *id);
  	iddata   = N_VGetArrayPointer_OpenMP(*id);
	data->k_bath=BathGasIndex(data);
	data->k_oxidizer=oxidizerIndex(data);
	data->k_OH=OHIndex(data);
	data->k_HO2=HO2Index(data);
    //use char* pointer to get the address
    for(int i=0;i<2;i++){
        specPtr[i] = data->dropSpec[i];
    }
	data->k_drop[0]=specIndex(data,specPtr[0]);
    data->k_drop[1]=specIndex(data,specPtr[1]);

	printf("Oxidizer index: %lu\n",data->k_oxidizer);
	printf("Bath gas index:%lu\n",data->k_bath);
	printf("Droplet species index:%lu,%lu\n",data->k_drop[0],data->k_drop[1]);
	for (size_t i = 1; i <=data->npts; i++) {
		/*Algebraic variables: indicated by ZERO.*/
		Rid(i)=ZERO;
		Pid(i)=ZERO;
		Yid(i,data->k_bath)=ZERO;
		Mdotid(i)=ZERO;
	}	
	Mdotid(1)=ZERO;
	Rid(1)=ONE;
	Yid(1,data->k_drop[0])=ZERO;
    Yid(1,data->k_drop[1])=ZERO;
	Tid(data->npts)=ZERO;
	if(data->constantPressure){
		Pid(data->npts)=ONE;
	}else{
		Pid(data->npts)=ZERO;
		Rid(data->npts)=ONE;
	}
	if(data->dirichletInner){
		Tid(1)=ZERO; 
	}else{
		Tid(1)=ONE;
	}
	if(data->dirichletOuter){
		for (size_t k = 1; k <=data->nsp; k++) {
			if(k!=data->k_bath){
				Yid(data->npts,k)=ONE;
			}
			Yid(1,k)=ZERO;
			Tid(data->npts)=ONE;
		}
	}else{
		for (size_t k = 1; k <=data->nsp; k++){
			Yid(1,k)=ZERO;
			Yid(data->npts,k)=ZERO;
			Tid(data->npts)=ZERO;
		}
	}
	return(0);
}

inline double calc_area(double x,int* i){
	switch (*i) {
		case 0:
			return(ONE);
		case 1:
			return(x);
		case 2:
			return(x*x);
		default:
			return(ONE);
	}
}

void readInitialCondition(FILE* input, double* ydata, const size_t nvar, const size_t nr, const size_t nPts, double Rg){

	FILE* output;output=fopen("test.dat","w");

	size_t bufLen=10000;
	size_t nRows=0;
	size_t nColumns=nvar;

	char buf[bufLen];
	char buf1[bufLen];
	char comment[1];
	char *ret;

	while (fgets(buf,bufLen, input)!=NULL){
		comment[0]=buf[0];
		if(strncmp(comment,"#",1)!=0){
			nRows++;
		}
	}
	rewind(input);

	printf("nRows: %ld\n", nRows);
	double y[nRows*nColumns];

	size_t i=0;
	while (fgets(buf,bufLen, input)!=NULL){
		comment[0]=buf[0];
		if(strncmp(comment,"#",1)==0){
		}
		else{
			ret=strtok(buf,"\t");
			size_t j=0;
			y[i*nColumns+j]=(double)(atof(ret));
			j++;
			while(ret!=NULL){
				ret=strtok(NULL,"\t");
				if(j<nColumns){
					y[i*nColumns+j]=(double)(atof(ret));
				}
				j++;
			}
			i++;
		}
	}

	for (i = 0; i < nRows; i++) {
		for (size_t j = 0; j < nColumns; j++) {
			fprintf(output, "%15.6e\t",y[i*nColumns+j]);
		}
		fprintf(output, "\n");
	}
	fclose(output);

	//double xOld[nRows],xNew[nPts],ytemp[nPts];
	double xOld[nRows],xNew[nPts],ytemp[nRows];

	for (size_t j = 0; j < nRows; j++) {
		xOld[j]=y[j*nColumns+nr];
	}

	//double dx=(xOld[nRows-1] - xOld[0])/((double)(nPts)-1.0e0);
	//for (size_t j = 0; j < nPts-1; j++) {
	//	xNew[j]=(double)j*dx + xOld[0];
	//}
	//xNew[nPts-1]=xOld[nRows-1];
	double dX0;
	int NN = nPts-2;
	
	dX0 = (xOld[nRows-1]-xOld[0])*(pow(Rg,1.0/NN) - 1.0)/(pow(Rg,(NN+1.0)/NN) - 1.0);
	xNew[0] = xOld[0];
	for (size_t i = 1; i < nPts-1; i++) {
		xNew[i]=xNew[i-1]+dX0*pow(Rg,(i-1.0)/NN);
	}
	xNew[nPts-1] = xOld[nRows-1];

	gsl_interp_accel* acc;
	gsl_spline* spline;
	acc = gsl_interp_accel_alloc();
	spline = gsl_spline_alloc(gsl_interp_steffen, nRows);

	for (size_t j = 0; j < nColumns; j++) {

		for (size_t k = 0; k < nRows; k++) {
			ytemp[k]=y[j+k*nColumns];
		}

		gsl_spline_init(spline,xOld,ytemp,nRows);

		for (size_t k = 0; k < nPts; k++) {
		      	ydata[j+k*nColumns]=gsl_spline_eval(spline,xNew[k],acc);
		}
	}
	
	gsl_interp_accel_free(acc);
	gsl_spline_free(spline);

}

double systemMass(double* ydata,UserData data){
	double mass=0.0e0;
	double rho;
	for (size_t i = 2; i <=data->npts; i++) {
  		data->gas->setState_TPY(T(i), P(i), &Y(i,1));   
		rho=data->gas->density();                                           
		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&m);
		mass+=rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
	}
	return(mass);
}


int initializePsiGrid(double* ydata, double* psidata, UserData data){

	double rho,rhom;
	/*Create a psi grid that corresponds CONSISTENTLY to the spatial grid
	 * "R" created above. Note that the Lagrangian variable psi has  units
	 * of kg. */
	psi(1)=ZERO;
	for (size_t i = 2; i <=data->npts; i++) {
  		data->gas->setState_TPY(T(i), P(i), &Y(i,1));   
		rho=data->gas->density();                                           
  		data->gas->setState_TPY(T(i-1), P(i-1), &Y(i-1,1));   
		rhom=data->gas->density();                                           
		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&data->metric);
		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
		psi(i)=psi(i-1)+(R(i)-R(i-1))*(rho*calc_area(R(i),&data->metric)+rhom*calc_area(R(i-1),&data->metric))/TWO;
	}

	/*The mass of the entire system is the value of psi at the last grid
	 * point. Normalize psi by this mass so that it varies from zero to
	 * one. This makes psi dimensionless. So the mass needs to be
	 * multiplied back in the approporiate places in the governing
	 * equations so that units match.*/
	data->mass=psi(data->npts);
	for (size_t i = 1; i <=data->npts; i++) {
		psi(i)=psi(i)/data->mass;
	}
	return(0);
}


int setInitialCondition(N_Vector* y, 
	    		N_Vector* ydot, 
            		UserData data){
	double* ydata;
	double* ydotdata;
	double* psidata;
	double* innerMassFractionsData, Rd, massDrop;
	double rhomhalf, lambdamhalf, YVmhalf[data->nsp];
	double f=ZERO;
	double g=ZERO;

	double perturb,rho;
	double epsilon=ZERO;
	int m,ier;
  	ydata    = N_VGetArrayPointer_OpenMP(*y);
  	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
	innerMassFractionsData =  data->innerMassFractions;
	Rd = data->Rd;
	massDrop = data->massDrop;

	//Mass of droplet
	//massDrop=1.0/3.0*Rd*Rd*Rd*997.0; //TODO: The density of the droplet should be a user input
    //massDrop=1.0/3.0*Rd*Rd*Rd*684.0; //TODO:The density of the droplet(n-heptane) should be a user input 
	massDrop = 1.0/3.0*Rd*Rd*Rd*data->dropRho;
    if(data->adaptiveGrid){
  		psidata  = data->grid->xOld;
	}
	else{
  		psidata  = data->uniformGrid;
	}

	m=data->metric;

	data->innerTemperature=data->initialTemperature;
	for (size_t k = 1; k <=data->nsp; k++) {
  		innerMassFractionsData[k-1]=data->gas->massFraction(k-1);
	}

	//Define Grid:
	double dR=(data->domainLength)/((double)(data->npts)-1.0e0);
	double dv=(pow(data->domainLength,1+data->metric)-pow(data->firstRadius*data->domainLength,1+data->metric))/((double)(data->npts)-1.0e0);
	
    //Initialize the R(i),T(i)(data->initialTemperature),Y(i,k),P(i)
    for (size_t i = 1; i <=data->npts; i++) {
		if(data->metric==0){
			R(i)=Rd+(double)((i-1)*dR);
		}else{
			if(i==1){
				R(i)=ZERO;
			}else if(i==2){
				R(i)=data->firstRadius*data->domainLength;
			}else{
				R(i)=pow(pow(R(i-1),1+data->metric)+dv,1.0/((double)(1+data->metric)));
			}
		}
		T(i)=data->initialTemperature;
		for (size_t k = 1; k <=data->nsp; k++) {
  	        	Y(i,k)=data->gas->massFraction(k-1);	//Indexing different in Cantera
		}
  	        P(i)=data->initialPressure*Cantera::OneAtm;
	}	
	R(data->npts)=data->domainLength+data->Rd;
    
//    /********** test R(i) and the volumn between grids *****************/
//    printf("Print the first 4 R(i)s and volumn between them: \n") ; 
//    printf("Grids: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e \n",R(1),R(2),R(3),R(4));
//    double v1,v2,v3,v4,v5; 
//    v1 =pow(R(2),1+data->metric) - pow(R(1),1+data->metric);
//    v2 =pow(R(3),1+data->metric) - pow(R(2),1+data->metric);
//    v3 =pow(R(4),1+data->metric) - pow(R(3),1+data->metric);
//    v4 =pow(R(5),1+data->metric) - pow(R(4),1+data->metric); 
//    v5 =pow(R(6),1+data->metric) - pow(R(5),1+data->metric);
//    printf("Volumn: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e,5th:%2.6e\n",v1,v2,v3,v4,v5) ; 


	double Tmax;
	double Tmin=data->initialTemperature;
	double w=data->mixingWidth;
	double YN2=ZERO;
	double YO2=ZERO;
	double YFuel,YOxidizer,sum;

	//if(data->problemType==0){
	//	data->gas->equilibrate("HP");
	//	data->maxTemperature=data->gas->temperature();
	//}
	//else if(data->problemType==1){
	//	/*Premixed Combustion: Equilibrium products comprise ignition
	//	 * kernel at t=0. The width of the kernel is "mixingWidth"
	//	 * shifted by "shift" from the center.*/
	//	data->gas->equilibrate("HP");
	//	Tmax=data->gas->temperature();
	//	for (size_t i = 1; i <=data->npts; i++) {
	//		g=HALF*(tanh((R(i)-data->shift)/w)+ONE);	//increasing function of x
	//		f=ONE-g;					//decreasing function of x
	//		T(i)=(Tmax-Tmin)*f+Tmin;
	//		for (size_t k = 1; k <=data->nsp; k++) {
	//			Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);
	//		}
	//	}
	//	if(data->dirichletOuter){
	//		T(data->npts)=data->wallTemperature;
	//	}
	//}
	//else if(data->problemType==2){
	//	FILE* input;
	//	if(input=fopen("initialCondition.dat","r")){
	//		readInitialCondition(input, ydata, data->nvar, data->nr, data->npts);
	//		fclose(input);
	//	}
	//	else{
	//		printf("file initialCondition.dat not found!\n");
	//		return(-1);
	//	}
	//}

	initializePsiGrid(ydata,psidata,data);

	if(data->adaptiveGrid){
	//	if(data->problemType!=0){
	//		data->grid->position=maxGradPosition(ydata, data->nt, data->nvar,
	//						 data->grid->xOld, data->npts);
	//		//data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
	//		//				 data->grid->xOld, data->npts);
	//	}
	//	else{
	//	}
		if(data->problemType!=3){
			data->grid->position=0.0e0;
			double x=data->grid->position+data->gridOffset*data->grid->leastMove;
			printf("New grid center:%15.6e\n",x);

//            /**************** TEST THE data->grid->xOld *******************/
//            double* ptr3 = data->grid->xOld ; 
//            printf("SetInitialCondition function is called,before reGrid,Start print the first 5 elements of the xOld array : \n");
//            printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr3[0],ptr3[1],ptr3[2],ptr3[3],ptr3[4]);
			
            ier=reGrid(data->grid, x);
			if(ier==-1)return(-1);

//            /**************** TEST THE data->grid->xOld *******************/
//            double* ptr = data->grid->xOld ; 
//            printf("SetInitialCondition function is called,after reGrid,Start print the first 5 elements of the xOld array : \n");
//            printf("1st:%.6f, 2nd:%.6f, 3rd:%.6f, 4th:%.6f, 5th:%.6f.\n",ptr[0],ptr[1],ptr[2],ptr[3],ptr[4]);


			updateSolution(ydata, ydotdata, data->nvar,
			               data->grid->xOld,data->grid->x,data->npts);
			storeGrid(data->grid->x,data->grid->xOld,data->npts);
		}
	}
	else{
		double Rg = data->Rg, dpsi0;
		int NN = data->npts-2;
		double psiNew[data->npts];
		
		dpsi0 = (pow(Rg,1.0/NN) - 1.0)/(pow(Rg,(NN+1.0)/NN) - 1.0);
		psiNew[0] = 0;
		for (size_t i = 1; i < data->npts-1; i++) {
			psiNew[i]=psiNew[i-1]+dpsi0*pow(Rg,(i-1.0)/NN);
		}
		psiNew[data->npts-1] = 1.0;
		printf("Last point:%15.6e\n",psiNew[data->npts-1]);
		updateSolution(ydata, ydotdata, data->nvar,
		               data->uniformGrid,psiNew,data->npts);
		storeGrid(psiNew,data->uniformGrid,data->npts);

		//double psiNew[data->npts];
		//double dpsi=1.0e0/((double)(data->npts)-1.0e0);
		//for (size_t i = 0; i < data->npts; i++) {
		//	psiNew[i]=(double)(i)*dpsi;
		//}
		//printf("Last point:%15.6e\n",psiNew[data->npts-1]);
		//updateSolution(ydata, ydotdata, data->nvar,
		//               data->uniformGrid,psiNew,data->npts);
		//storeGrid(psiNew,data->uniformGrid,data->npts);
	}

	if(data->problemType==0){
		data->gas->equilibrate("HP");
		data->maxTemperature=data->gas->temperature();
	}
	else if(data->problemType==1){
		/*Premixed Combustion: Equilibrium products comprise ignition
		 * kernel at t=0. The width of the kernel is "mixingWidth"
		 * shifted by "shift" from the center.*/
		data->gas->equilibrate("HP");
		Tmax=data->gas->temperature();
		for (size_t i = 1; i <=data->npts; i++) {
			g=HALF*(tanh((R(i)-data->shift)/w)+ONE);	//increasing function of x
			f=ONE-g;					//decreasing function of x
			T(i)=(Tmax-Tmin)*f+Tmin;
			for (size_t k = 1; k <=data->nsp; k++) {
				Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);
			}
		}
		if(data->dirichletOuter){
			T(data->npts)=data->wallTemperature;
		}
	}
	else if(data->problemType==2){
		FILE* input;
		if(input=fopen("initialCondition.dat","r")){
			readInitialCondition(input, ydata, data->nvar, data->nr, data->npts, data->Rg);
			fclose(input);
		}
		else{
			printf("file initialCondition.dat not found!\n");
			return(-1);
		}
		initializePsiGrid(ydata,psidata,data);
	}
	else if(data->problemType==3){
		FILE* input;
		if(input=fopen("restart.bin","r")){
			readRestart(y, ydot, input, data);
			fclose(input);
			printf("Restart solution loaded!\n");
			printf("Problem starting at t=%15.6e\n",data->tNow);
			return(0);
		}
		else{
			printf("file restart.bin not found!\n");
			return(-1);
		}
	}

	if(data->reflectProblem){
		double temp;
		int j=1;
		while (data->npts+1-2*j>=0) {
			temp=T(j);
			T(j)=T(data->npts+1-j);
			T(data->npts+1-j)=temp;
			for (size_t k = 1; k <=data->nsp; k++) {
				temp=Y(j,k);
				Y(j,k)=Y(data->npts+1-j,k);
				Y(data->npts+1-j,k)=temp;
			}
			j=j+1;
		}
	}

	/*Floor small values to zero*/
	for (size_t i = 1; i <=data->npts; i++) {
		for (size_t k = 1; k <=data->nsp; k++) {
			if(fabs(Y(i,k))<=data->massFractionTolerance){
				Y(i,k)=0.0e0;
			}
		}
	}

	//Set grid to location of maximum curvature calculated by r instead of x:
	if(data->adaptiveGrid){
		//data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
		//		  		data->grid->x, data->npts);
        int maxCurvII = 0; 
        maxCurvII = maxCurvIndexR(ydata,data->nt,data->nvar,data->nr,data->npts);
        
        data->grid->position = data->grid->x[maxCurvII];
		ier=reGrid(data->grid, data->grid->position);
		updateSolution(ydata, ydotdata, data->nvar,
		               data->grid->xOld,data->grid->x,data->npts);
		storeGrid(data->grid->x,data->grid->xOld,data->npts);
        
        /******** Test the maxCurvPosition and related variables *******/
        printf("The maxCurvPositition(based on r) is : %.6f,Temperature is : %.3f \n",data->grid->position,T(maxCurvII+1));
	}


	/*Ensure consistent boundary conditions*/	
	//T(1)=T(2);
	//for (size_t k = 1; k <=data->nsp; k++) {
	//	Y(1,k)=Y(2,k);
	//	Y(data->npts,k)=Y(data->npts-1,k);
	//}

	return(0);
}

inline double Qdot(double* t, 
	    	   double* x, 
	    	   double* ignTime, 
	    	   double* kernelSize, 
	    	   double* maxQdot){
	double qdot;
	if(*x<=*kernelSize){
		if((*t)<=(*ignTime)){
			qdot=(*maxQdot);
		}
		else{
			qdot=0.0e0;
		}
	}else{
		qdot=0.0e0;
	}
	return(qdot);
}

inline void setGas(UserData data, 
	    double *ydata, 
	    size_t gridPoint){
  	data->gas->setTemperature(T(gridPoint)); 	        
  	data->gas->setMassFractions_NoNorm(&Y(gridPoint,1)); 	
  	data->gas->setPressure(P(gridPoint));
}

void getTransport(UserData data, 
		  double *ydata, 
		  size_t gridPoint,
		  double *rho,
		  double *lambda,
		  double YV[]){

	double YAvg[data->nsp],
		 XLeft[data->nsp],
		 XRight[data->nsp],
		 gradX[data->nsp];

	setGas(data,ydata,gridPoint);
	data->gas->getMoleFractions(XLeft);                                     
	setGas(data,ydata,gridPoint+1);
	data->gas->getMoleFractions(XRight);                                     
                                                                             
	for (size_t k = 1; k <=data->nsp; k++) {                                          
		YAvg(k)=HALF*(Y(gridPoint,k)+
			      Y(gridPoint+1,k));                                 
		gradX(k)=(XRight(k)-XLeft(k))/
			 (R(gridPoint+1)-R(gridPoint));                       
	}                                                                    
	double TAvg = HALF*(T(gridPoint)+T(gridPoint+1));
	double gradT=(T(gridPoint+1)-T(gridPoint))/
		       (R(gridPoint+1)-R(gridPoint));                                  


  	data->gas->setTemperature(TAvg); 	        
  	data->gas->setMassFractions_NoNorm(YAvg); 	
  	data->gas->setPressure(P(gridPoint));

	*rho=data->gas->density();                                       
	*lambda=data->trmix->thermalConductivity();                      
	data->trmix->getSpeciesFluxes(1,&gradT,data->nsp,                     
				      gradX,data->nsp,YV);                
	//setGas(data,ydata,gridPoint);
}

int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){

	/*Declare and fetch nvectors and user data:*/

  	double *ydata, *ydotdata, *resdata, *psidata, *innerMassFractionsData;
  	
	UserData data;
	data = (UserData)user_data;
	size_t npts=data->npts;
	size_t nsp=data->nsp;
	size_t k_bath = data->k_bath;
	size_t k_drop[2] = {data->k_drop[0],data->k_drop[1]};

	ydata 	= N_VGetArrayPointer_OpenMP(y); 
	ydotdata= N_VGetArrayPointer_OpenMP(ydot); 
	resdata = N_VGetArrayPointer_OpenMP(res);
	if(data->adaptiveGrid==1){
		psidata = data->grid->x;
	}else{
		psidata = data->uniformGrid;
	}

	innerMassFractionsData =  data->innerMassFractions;

	/* Grid stencil:*/

	/*-------|---------*---------|---------*---------|-------*/
	/*-------|---------*---------|---------*---------|-------*/
	/*-------|---------*---------|---------*---------|-------*/
	/*-------m-------mhalf-------j-------phalf-------p-------*/
	/*-------|---------*---------|---------*---------|-------*/
	/*-------|---------*---------|---------*---------|-------*/
	/*-------|<=======dxm=======>|<=======dxp=======>|-------*/
	/*-------|---------*<======dxav=======>*---------|-------*/
	/*-------|<================dxpm=================>|-------*/

	/* Various variables defined for book-keeping and storing previously
	 * calculated values:
	 * rho		: densities at points  m, mhalf, j, p, and phalf.
	 * area		: the matric at points m, mhalf, j, p, and phalf.
	 * m 		: exponent that determines geometry;
	 * lambda	: thermal conductivities at mhalf and phalf.
	 * mdot		: mass flow rate at m, j, and p.
	 * X		: mole fractions at j and p.
	 * YV		: diffusion fluxes at mhalf and phalf.
	 * Tgrad	: temperature gradient at mhalf and phalf.
	 * Tav		: average temperature between two points.
	 * Pav		: average pressure between two points.
	 * Yav		: average mass fractions between two points.
	 * Xgradhalf	: mole fraction gradient at j.
	 * Cpb		: mass based bulk specific heat.
	 * tranTerm	: transient terms.
	 * advTerm	: advection terms.
	 * diffTerm	: diffusion terms.
	 * srcTerm	: source terms.
	 */

	double rhomhalf,  rhom, lambdamhalf, YVmhalf[nsp], 
		 rho, 
		 rhophalf, lambdaphalf, YVphalf[nsp],
		 Cpb, Cvb,  Cp[nsp], Cpl[2], dHvl[2], rhol,  wdot[nsp],   enthalpy[nsp], energy[nsp],
		 tranTerm, diffTerm,    srcTerm, advTerm,
		 area,areamhalf,areaphalf,aream,areamhalfsq,areaphalfsq;
	
	/*Aliases for difference coefficients:*/
	double cendfm, cendfc, cendfp;
	/*Aliases for various grid spacings:*/
	double dpsip, dpsiav, dpsipm, dpsim, dpsimm;
	dpsip=dpsiav=dpsipm=dpsim=dpsimm=ONE;
	//double mass, mdotIn;
	double mass, massDrop;
	double sum, sum1, sum2, sum3;

	size_t j,k;
	int m;
	m=data->metric;					//Unitless
	mass=data->mass;				//Units: kg
	//massDrop=data->massDrop;		//Units: kg
	//mdotIn=data->mdot*calc_area(R(npts),&m);	//Units: kg/s

//	/*evaluate properties at j=1*************************/
	setGas(data,ydata,1);
	rhom=data->gas->density();                                       
    Cpb=data->gas->cp_mass();       	//J/kg/K
    Cvb=data->gas->cv_mass();       	//J/kg/K
	//TODO: Create user input model for these. They should not be constant
	//Cpl=4.182e03;						//J/kg/K for liquid water
	//Cpl=  5.67508633e-07*pow(T(1),4) - 7.78060597e-04*pow(T(1),3) + 4.08310544e-01*pow(T(1),2)  //J/kg/K for liquid water  
    //     - 9.62429538e+01*T(1)        + 1.27131046e+04;
    
    /*Update specific heat:Cpl for liquid proprane&n-heptane&dodecane*/
   /*Based on the existiong data,a linear expression is used*/
   // Cpl= 12.10476*T(1) - 746.60143; // Unit:J/(kg*K), Need to be improved later
    //Cpl[0] = 2.423*T(1)+1661.074; //Unit:J/(kg*K),range:1atm,propane 
    Cpl[0] = 3.336*T(1)+1289.5; //Unit:J/(kg*K),range:1atm,n-Heptane
    Cpl[1] = 3.815*T(1)+1081.8; //Unit:J/(kg*K),range:1atm,n-Dodecane
    double Cp_l = 0.0; 
    for(size_t i=0;i<=1;i++){
        Cp_l=Cp_l+Cpl[i]*data->dropMassFrac[i];
    }

	//dHvl=2.260e6;						//J/kg	heat of vaporization of water
	//double Tr = T(1)/647.1;				//Reduced Temperature: Droplet Temperature divided by critical temperature of water
	//dHvl= 5.2053e07*pow(1-Tr,0.3199 - 0.212*Tr + 0.25795*Tr*Tr)/18.01;	//J/kg	latent heat of vaporization of water
    
    double Tr1 = T(1)/369.8;            //Reduced Temperature;
    dHvl[0] = 6.32716e5*exp(-0.0208*Tr1)*pow((1-Tr1),0.3766);    //Unit:J/kg,Latent Heat of Vaporizaiton of Liquid n-propane,NIST

    double Tr2 = T(1)/540.2;            //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST
    dHvl[1] = 5.366e5*exp(-0.2831*Tr2) * pow((1-Tr2),0.2831);    //Unit:J/kg, latent heat of vaporization of liquid n-heptane, Source: NIST
    
    double dHv_l = 0.0; 
    for(size_t i=0;i<=1;i++){
        dHv_l=dHv_l+dHvl[i]*data->dropMassFrac[i];
    }

    /*Following section is related to the property of water*/
    //rhol = 997.0;	
	//massDrop=1.0/3.0*R(1)*R(1)*R(1)*997.0; //TODO: The density of the droplet should be a user input
    
    /*Following section is related to the property of liquid n-heptane and n-dodecane (mainly density rho)*/
    /*Density of these two species should be temperature dependent*/
    //data->dropDens[0] = -1.046*T(1)+823.794;       //Unit:kg/m^3,density of propane @1atm
    data->dropDens[0] = -0.858*T(1)+933.854;       //Unit:kg/m^3,density of n-Heptane @1atm
    data->dropDens[1] = -0.782*T(1)+979.643;       //Unit:kg/m^3,density of n-Dodecane @1atm
    //rhol = data->dropRho;                          //Unit:kg/m^3
    rhol = 0.0; 
    for(size_t i=0;i<=1;i++){
       rhol = rhol + data->dropMassFrac[i]*data->dropDens[i];
    }
    massDrop = 1.0/3.0*R(1)*R(1)*R(1)*rhol;  //Unit:kg, this is the mass of liquid droplet  

    aream= calc_area(R(1),&m);	    
    
 
	/*******************************************************************/
	/*Calculate values at j=2's m and mhalf*****************************/

	getTransport(data, ydata, 1, &rhomhalf,&lambdamhalf,YVmhalf);
	areamhalf= calc_area(HALF*(R(1)+R(2)),&m);                                 
	aream = calc_area(R(1),&m);                                 
	areamhalfsq= areamhalf*areamhalf;                                   
                                                                             
	/*******************************************************************/

	/*Fill up res with left side (center) boundary conditions:**********/
	/*We impose zero fluxes at the center:*/                        


	/*Mass:*/                                                            
	if(data->quasiSteady){
		Rres(1)=Rdot(1); // Should remain satisfied for quasi-steady droplet
	}else{
		Rres(1)=Rdot(1) + Mdot(1)/rhol/aream; 
	}
                                                                        
	/*Species:*/                                                    
	sum=ZERO;
	
	//TODO: Antoine Parameters should be part of user input, not hardcoded
	//Yres(1,k_drop)=Y(1,k_drop) - 1.0e5*pow(10,4.6543-(1435.264/(T(1)-64.848)))*data->gas->molecularWeight(k_drop-1) 
	//			    / P(1) / data->gas->meanMolecularWeight();
	//Yres(1,k_drop)=Y(1,k_drop) - 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)))*data->gas->molecularWeight(k_drop-1) 
	//			    / P(1) / data->gas->meanMolecularWeight();
	//Yres(1,k_drop)=Y(1,k_drop)-2.566785e-02;
	
    /*Following using the Antoine Parameters to calculate the pressure of volatile component(n-heptane)*/
    //Antoine parameter from NIST website
    double p_i[2]={0.0,0.0} ; 
    
    p_i[0] = 1.0e5*pow(10,4.53678-(1149.36/(T(1)+24.906)))*data->dropMassFrac[0] ;    //unit:Pa,For N-PROPANE
    p_i[1] = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199)))*data->dropMassFrac[1] ;    //unit:Pa,FOR N-HEPTANE 
    //p_i = 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)));     //Unit:Pa,FOR WATER
    //Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1) 
	//			    / P(1) / data->gas->meanMolecularWeight();
    for(size_t i=0;i<=1;i++){
        Yres(1,k_drop[i])=Y(1,k_drop[i])-p_i[i]*data->gas->molecularWeight(k_drop[i]-1)/(P(1)*data->gas->meanMolecularWeight());           
    }

	
    sum=sum+Y(1,k_drop[0])+Y(1,k_drop[1]);

	//Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass
    Mdotres(1)= Mdot(1) - ((YVmhalf(k_drop[0])+YVmhalf(k_drop[1])) *areamhalf) / (1-Y(1,k_drop[0])-Y(1,k_drop[1]) ); //Evaporating mass
	for (k = 1; k <=nsp; k++) {                                        
		if(k!=k_bath && k!=k_drop[0] && k!=k_drop[1]){
			//TODO: May need to include chemical source term
			Yres(1,k)=Y(1,k)*Mdot(1) + YVmhalf(k)*areamhalf;
			//Yres(1,k)=YVmhalf(k)*areamhalf;
		    sum=sum+Y(1,k);
			//if(fabs(Mdot(1))>1e-14){
		      		///Yres(1,k)=innerMassFractionsData[k-1]-
					///  Y(1,k)-
					///  (YVmhalf(k)*areamhalf)/Mdot(1);
					
			//}
			//else{
		    //  		//Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1];
		    //  		Yres(1,k)=Y(2,k)-Y(1,k);
			//	/*The below flux boundary condition makes the
			//	 * problem more prone to diverge. How does one
			//	 * fix this?*/
		    //  		//Yres(1,k)=YVmhalf(k);
			//}
		}
	}                                                                 
	Yres(1,k_bath)=ONE-sum-Y(1,k_bath);                                     
	
	/*Energy:*/                                                     
	if(data->dirichletInner){
		//Tres(1)=Tdot(1);
		//Tres(1)=T(1)-data->innerTemperature;
        /*Following code should be revised in the future*/
		//Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHv_l - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop));
	}else{
		//Tres(1)=Tdot(1) - 
		//	(rhomhalf*areamhalfsq*lambdamhalf*(T(2)-T(1))/((psi(2)-psi(1))*mass) - Mdot(1)*dHvl)
		//	/ (massDrop * Cpl);
		
        //Tres(1)=Tdot(1) - 
		//	(areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl)
		//	/ (massDrop * Cpl);
		Tres(1)=Tdot(1) - 
			(areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHv_l)
			/ (massDrop * Cp_l);



		//Tres(1)=T(2)-T(1);
		//Tres(1)=Tdot(1) - (Pdot(1)/(rhom*Cpb))
		//	+(double)(data->metric+1)*(rhomhalf*lambdamhalf*areamhalfsq*(T(2)-T(1))/psi(2)-psi(1));
	}

	/*Pressure:*/                                                        
	Pres(1)=P(2)-P(1);
	/*Fill up res with governing equations at inner points:*************/
	for (j = 2; j < npts; j++) {                                         
                                                                             
		/*evaluate various mesh differences*///                      
        	dpsip =        (psi(j+1) - psi(j)  )*mass;                            
        	dpsim =        (psi(j)   - psi(j-1))*mass;                            
        	dpsiav =  HALF*(psi(j+1) - psi(j-1))*mass;                            
        	dpsipm =       (psi(j+1) - psi(j-1))*mass;                            
		/***********************************///                      
                                                                             
		/*evaluate various central difference coefficients*/       
        	cendfm = - dpsip / (dpsim*dpsipm);                               
        	cendfc =   (dpsip-dpsim) / (dpsip*dpsim);                         
        	cendfp =   dpsim / (dpsip*dpsipm);                               
		/**************************************************/
                                                                   
                                                                  
		/*evaluate properties at j*************************/
		setGas(data,ydata,j);
		rho=data->gas->density();		//kg/m^3
        	Cpb=data->gas->cp_mass();       	//J/kg/K
        	Cvb=data->gas->cv_mass();       	//J/kg/K
		data->gas->getNetProductionRates(wdot); //kmol/m^3
		data->gas->getEnthalpy_RT(enthalpy);	//unitless
		data->gas->getCp_R(Cp);			//unitless
		area = calc_area(R(j),&m);	        //m^2
                                                                           
		/*evaluate properties at p*************************/
		getTransport(data, ydata, j, &rhophalf,&lambdaphalf,YVphalf);
		areaphalf= calc_area(HALF*(R(j)+R(j+1)),&m);	    
		areaphalfsq= areaphalf*areaphalf;		    
		/**************************************************///     

		/*Evaporating Mass*/
		Mdotres(j)=Mdot(j) - Mdot(j-1); 

		/*Mass:*/
		/* ∂r/∂ψ = 1/ρA */
		Rres(j)=((R(j)-R(j-1))/dpsim)-(TWO/(rhom*aream+rho*area));

		/*Energy:*/
		/* ∂T/∂t = - ṁ(∂T/∂ψ) 
		 * 	   + (∂/∂ψ)(λρA²∂T/∂ψ) 
		 * 	   - (A/cₚ) ∑ YᵢVᵢcₚᵢ(∂T/∂ψ) 
		 * 	   - (1/ρcₚ)∑ ώᵢhᵢ 
		 * 	   + (1/ρcₚ)(∂P/∂t) */
		/*Notes: 
		 * λ has units J/m/s/K.
		 * YᵢVᵢ has units kg/m^2/s.
		 * hᵢ has units J/kmol, so we must multiply the enthalpy
		 * defined above (getEnthalpy_RT) by T (K) and the gas constant
		 * (J/kmol/K) to get the right units.
		 * cₚᵢ has units J/kg/K, so we must multiply the specific heat
		 * defined above (getCp_R) by the gas constant (J/kmol/K) and
		 * divide by the molecular weight (kg/kmol) to get the right
		 * units.
		 * */

		//enthalpy formulation:
		//tranTerm = Tdot(j) - (Pdot(j)/(rho*Cpb));
		tranTerm = Tdot(j);
		sum=ZERO;
		sum1=ZERO;
		for (k = 1; k <=nsp; k++) {
			sum=sum+wdot(k)*enthalpy(k);
			sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho
			     *HALF*(YVmhalf(k)+YVphalf(k));
		}
		sum=sum*Cantera::GasConstant*T(j);
		sum1=sum1*Cantera::GasConstant;
		diffTerm  =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
			      -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
			     /(dpsiav*Cpb) )
			  -(sum1*area*(cendfp*T(j+1)
			              +cendfc*T(j)
			              +cendfm*T(j-1))/Cpb);
		srcTerm   = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cpb); // Qdot is forced heating to cause ignition, sum is the conversion of chemical enthalpy to sensible enthalpy
		advTerm   = (Mdot(j)*(T(j)-T(j-1))/dpsim);
		Tres(j)= tranTerm - (Pdot(j)/(rho*Cpb))
			+advTerm 
			-diffTerm
			+srcTerm;

	//	//energy formulation:
	//	tranTerm = Tdot(j);
	//	sum=ZERO;
	//	sum1=ZERO;
	//	sum2=ZERO;
	//	sum3=ZERO;
	//	for (k = 1; k <=nsp; k++) {
	//		energy(k)=enthalpy(k)-ONE;
	//		sum=sum+wdot(k)*energy(k);
	//		sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho
	//		     *HALF*(YVmhalf(k)+YVphalf(k));
	//		sum2=sum2+(YVmhalf(k)/data->gas->molecularWeight(k-1));
	//		sum3=sum3+(YVphalf(k)/data->gas->molecularWeight(k-1));
	//	}
	//	sum=sum*Cantera::GasConstant*T(j);
	//	sum1=sum1*Cantera::GasConstant;
	//	diffTerm  =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
	//		      -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
	//		     /(dpsiav*Cvb) )
	//		  -(sum1*area*(cendfp*T(j+1)
	//		              +cendfc*T(j)
	//		              +cendfm*T(j-1))/Cvb);
	//	srcTerm   = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cvb);
	//	advTerm   = (mdotIn*(T(j)-T(j-1))/dpsim);
	//	advTerm = advTerm + (Cantera::GasConstant*T(j)*area/Cvb)*((sum3-sum2)/dpsiav);
	//	Tres(j)= tranTerm
	//		+advTerm 
	//		-diffTerm
	//		+srcTerm;

		/*Species:*/
		/* ∂Yᵢ/∂t = - ṁ(∂Yᵢ/∂ψ) 
		 * 	    - (∂/∂ψ)(AYᵢVᵢ) 
		 * 	    + (ώᵢWᵢ/ρ)  */

		sum=ZERO;
		for (k = 1; k <=nsp; k++) {
			if(k!=k_bath){
				tranTerm  = Ydot(j,k);
				diffTerm  = (YVphalf(k)*areaphalf
					    -YVmhalf(k)*areamhalf)/dpsiav;
				srcTerm   = wdot(k)
			     		    *(data->gas->molecularWeight(k-1))/rho;
				advTerm   = (Mdot(j)*(Y(j,k)-Y(j-1,k))/dpsim);

				Yres(j,k)= tranTerm 
					  +advTerm
					  +diffTerm
					  -srcTerm;

		      		sum=sum+Y(j,k);
			}
		}
		Yres(j,k_bath)=ONE-sum-Y(j,k_bath);

		/*Pressure:*/
		Pres(j) = P(j+1)-P(j);

		/*Assign values evaluated at p and phalf to m       
		 * and mhalf to save some cpu cost:****************/
		areamhalf=areaphalf;                                
		areamhalfsq=areaphalfsq;                            
		aream=area;
		rhom=rho;
		rhomhalf=rhophalf;                                  
		lambdamhalf=lambdaphalf;                            
		for (k = 1; k <=nsp; k++) {                         
			YVmhalf(k)=YVphalf(k);                      
		}                                                   
		/**************************************************/
		                                                    
	}                                                           
	/*******************************************************************///


	/*Fill up res with right side (wall) boundary conditions:***********/
	/*We impose zero fluxes at the wall:*/                               
	setGas(data,ydata,npts);
	rho=data->gas->density();                               
	area = calc_area(R(npts),&m);				   

	/*Mass:*/                                                          
	dpsim=(psi(npts)-psi(npts-1))*mass;
	Rres(npts)=((R(npts)-R(npts-1))/dpsim)-(TWO/(rhom*aream+rho*area));

	/*Energy:*/                                             
	if(data->dirichletOuter){
		//Tres(npts)=T(npts)-data->wallTemperature;
		Tres(npts)=Tdot(npts);
	}
	else{
		Tres(npts)=T(npts)-T(npts-1);
	}
                                                                
	/*Species:*/                                            
	sum=ZERO;
	if(data->dirichletOuter){
		for (k = 1; k <=nsp; k++) {                             
			if(k!=k_bath){
			    	Yres(npts,k)=Ydot(npts,k);
			    	sum=sum+Y(npts,k);
			}
		}                                                       
	}
	else{
		for (k = 1; k <=nsp; k++) {                             
			if(k!=k_bath){
			    	Yres(npts,k)=Y(npts,k)-Y(npts-1,k);
			    	//Yres(npts,k)=YVmhalf(k);         
			    	sum=sum+Y(npts,k);
			}
		}                                                       
	}
	Yres(npts,k_bath)=ONE-sum-Y(npts,k_bath);             
                                                                
	
	/*Pressure:*/                                                      
	if(data->constantPressure){
		Pres(npts)=Pdot(npts)-data->dPdt;                                     
	}
	else{
		Pres(npts)=R(npts)-data->domainLength;                     
		//Pres(npts)=Rdot(npts);                     
	}

	/*Evaporating Mass*/
	Mdotres(npts)=Mdot(npts) - Mdot(npts-1); 

	//for (j = 1; j <=npts; j++) {
	//	//for (k = 1; k <=nsp; k++) {
	//	//	Yres(j,k)=Ydot(j,k);
	//	//}
	//	//Tres(j)=Tdot(j);
	//}

	return(0);

}

void printSpaceTimeHeader(UserData data)
{
	fprintf((data->output), "%15s\t","#1");
	for (size_t k = 1; k <=data->nvar+1; k++) {
		fprintf((data->output), "%15lu\t",k+1);
	}
	fprintf((data->output), "%15lu\n",data->nvar+3);

	fprintf((data->output), "%15s\t%15s\t%15s\t","#psi","time(s)","dpsi");
	fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");
	for (size_t k = 1; k <=data->nsp; k++) {
		fprintf((data->output), "%15s\t",data->gas->speciesName(k-1).c_str());
	}
	fprintf((data->output), "%15s\t","Pressure(Pa)");
	fprintf((data->output), "%15s\n","Mdot (kg/s)");
}

void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
{
	double *ydata,*psidata;
  	ydata    = N_VGetArrayPointer_OpenMP(*y);

	if(data->adaptiveGrid){
  		psidata  = data->grid->x;
	}else{
  		psidata  = data->uniformGrid;
	}

	for (size_t i = 0; i < data->npts; i++) {
		fprintf(output, "%15.9e\t%15.9e\t",psi(i+1),t);
		if(i==0){
			fprintf(output, "%15.9e\t",psi(2)-psi(1));
		}
		else{
			fprintf(output, "%15.6e\t",psi(i+1)-psi(i));
		}
		for (size_t j = 0; j < data->nvar; j++) {
			if(j!=data->nvar-1){
				fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
			}else{
				fprintf(output, "%15.9e",ydata[j+i*data->nvar]);
			}
		}
		fprintf(output, "\n");
	}
	fprintf(output, "\n");
}

void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data){
	double *ydata,*psidata, *ydotdata;
	ydata = N_VGetArrayPointer_OpenMP(*y);
	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
	if(data->adaptiveGrid){
  		psidata  = data->grid->x;
	}else{
  		psidata  = data->uniformGrid;
	}
 	fwrite(&t, sizeof(t), 1, output);		//write time
 	fwrite(psidata, data->npts*sizeof(psidata), 1, output);	//write grid
 	fwrite(ydata, data->neq*sizeof(ydata), 1, output);	//write solution
 	fwrite(ydotdata, data->neq*sizeof(ydotdata), 1, output);	//write solutiondot
}

void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data){
	double *ydata,*psidata, *ydotdata;
	double t;
	if(data->adaptiveGrid){
  		psidata  = data->grid->x;
	}else{
  		psidata  = data->uniformGrid;
	}
	ydata = N_VGetArrayPointer_OpenMP(*y);
	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
 	fread(&t, sizeof(t), 1, input);
	data->tNow=t;
 	fread(psidata, data->npts*sizeof(psidata), 1, input);
 	fread(ydata, data->neq*sizeof(ydata), 1, input);
 	fread(ydotdata, data->neq*sizeof(ydotdata), 1, input);
	if(data->adaptiveGrid){
		storeGrid(data->grid->x,data->grid->xOld,data->npts);
	}
}

void printGlobalHeader(UserData data)
{
	fprintf((data->globalOutput), "%8s\t","#Time(s)");
	//fprintf((data->globalOutput), "%15s","S_u(exp*)(m/s)");
	fprintf((data->globalOutput), "%15s","bFlux(kg/m^2/s)");
	fprintf((data->globalOutput), "%16s","  IsothermPos(m)");
	fprintf((data->globalOutput), "%15s\t","Pressure(Pa)");
	fprintf((data->globalOutput), "%15s\t","Pdot(Pa/s)");
	fprintf((data->globalOutput), "%15s\t","gamma");
	fprintf((data->globalOutput), "%15s\t","S_u(m/s)");
	fprintf((data->globalOutput), "%15s\t","Tu(K)");
	fprintf((data->globalOutput), "\n");
}
//
//void printSpaceTimeRates(double t, N_Vector ydot, UserData data)
//{
//	double *ydotdata,*psidata;
//  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);
//  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);
//	for (int i = 0; i < data->npts; i++) {
//		fprintf((data->ratesOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
//		for (int j = 0; j < data->nvar; j++) {
//			fprintf((data->ratesOutput), "%15.6e\t",ydotdata[j+i*data->nvar]);
//		}
//		fprintf((data->ratesOutput), "\n");
//	}
//	fprintf((data->ratesOutput), "\n\n");
//}
//
void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data)
{
	double *ydata,*ydotdata, *innerMassFractionsData, *psidata;
	innerMassFractionsData =  data->innerMassFractions;

	if(data->adaptiveGrid){
  		psidata  = data->grid->x;
	}else{
  		psidata  = data->uniformGrid;
	}

	double TAvg, RAvg, YAvg, psiAvg;
  	ydata    = N_VGetArrayPointer_OpenMP(*y);
  	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
	TAvg=data->isotherm;
	double sum=ZERO;
	double dpsim,area,aream,drdt;
	double Cpb,Cvb,gamma,rho,flameArea,Tu;


	/*Find the isotherm chosen by the user*/
	size_t j=1;
	size_t jj=1;
	size_t jjj=1;
	double wdot[data->nsp];
	double wdotMax=0.0e0;
	double advTerm=0.0e0;
	psiAvg=0.0e0;
	if(T(data->npts)>T(1)){
		while (T(j)<TAvg) {
			j=j+1;
		}
		YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(data->npts,data->k_oxidizer);
		while (fabs((T(jj+1)-T(jj))/T(jj))>1e-08) {
			jj=jj+1;
		}

		setGas(data,ydata,jj);
		Tu=T(jj);
		rho=data->gas->density();
        	Cpb=data->gas->cp_mass();       	//J/kg/K
        	Cvb=data->gas->cv_mass();       	//J/kg/K
		gamma=Cpb/Cvb;
		
	}
	else{
		while (T(j)>TAvg) {
			j=j+1;
		}
		YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(1,data->k_oxidizer);
		while (fabs((T(data->npts-jj-1)-T(data->npts-jj))/T(data->npts-jj))>1e-08) {
			jj=jj+1;
		}

		setGas(data,ydata,data->npts-jj);
		Tu=T(data->npts-jj);
		rho=data->gas->density();
        	Cpb=data->gas->cp_mass();       	//J/kg/K
        	Cvb=data->gas->cv_mass();       	//J/kg/K
		gamma=Cpb/Cvb;
	}
	if(T(j)<TAvg){
		RAvg=((R(j+1)-R(j))/(T(j+1)-T(j)))*(TAvg-T(j))+R(j);
	}
	else{
		RAvg=((R(j)-R(j-1))/(T(j)-T(j-1)))*(TAvg-T(j-1))+R(j-1);
	}

	////Experimental burning speed calculation:	
	//int nMax=0;
	////nMax=maxCurvIndex(ydata, data->nt, data->nvar,
	////     	data->grid->x, data->npts);
	//nMax=maxGradIndex(ydata, data->nt, data->nvar,
	//     	data->grid->x, data->npts);
	//advTerm=(T(nMax)-T(nMax-1))/(data->mass*(psi(nMax)-psi(nMax-1)));
	//aream=calc_area(R(nMax),&data->metric);	    
	////setGas(data,ydata,nMax);
	////rho=data->gas->density();
	//psiAvg=-Tdot(nMax)/(rho*aream*advTerm);
	////if(t>data->ignTime){
	////	for(size_t n=2;n<data->npts;n++){
	////		setGas(data,ydata,n);
	////		data->gas->getNetProductionRates(wdot); //kmol/m^3
	////		advTerm=(T(n)-T(n-1))/(data->mass*(psi(n)-psi(n-1)));
	////		if(fabs(wdot[data->k_oxidizer-1])>=wdotMax){
	////			aream=calc_area(R(n),&data->metric);	    
	////			psiAvg=-Tdot(n)/(rho*aream*advTerm);
	////			wdotMax=fabs(wdot[data->k_oxidizer-1]);
	////		}
	////	}
	////}
	////else{
	////	psiAvg=0.0e0;
	////}


	//drdt=(RAvg-data->flamePosition[1])/(t-data->flameTime[1]);
	//data->flamePosition[0]=data->flamePosition[1];
	//data->flamePosition[1]=RAvg;
	//data->flameTime[0]=data->flameTime[1];
	//data->flameTime[1]=t;
	//flameArea=calc_area(RAvg,&data->metric);

	/*Use the Trapezoidal rule to calculate the mass burning rate based on
	 * the consumption of O2*/
	aream= calc_area(R(1)+1e-03*data->domainLength,&data->metric);	    
	for (j = 2; j <data->npts; j++) {
        	dpsim=(psi(j)-psi(j-1))*data->mass;                            
		area= calc_area(R(j),&data->metric);	    
		sum=sum+HALF*dpsim*((Ydot(j-1,data->k_oxidizer)/aream)
				   +(Ydot(j,data->k_oxidizer)/area));
		aream=area;
	}

	//double maxOH,maxHO2;
	//maxOH=0.0e0;
	//maxHO2=0.0e0;
	//for(j=1;j<data->npts;j++){
	//	if(Y(j,data->k_OH)>maxOH){
	//		maxOH=Y(j,data->k_OH);
	//	}
	//}
	//for(j=1;j<data->npts;j++){
	//	if(Y(j,data->k_HO2)>maxHO2){
	//		maxHO2=Y(j,data->k_HO2);
	//	}
	//}

	fprintf((data->globalOutput), "%15.6e\t",t);
	//fprintf((data->globalOutput), "%15.6e\t",psiAvg);
	fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/YAvg);
	fprintf((data->globalOutput), "%15.6e\t",RAvg);
	fprintf((data->globalOutput), "%15.6e\t",P(data->npts));
	fprintf((data->globalOutput), "%15.6e\t",Pdot(data->npts));
	fprintf((data->globalOutput), "%15.6e\t",gamma);
	fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/(YAvg*rho));
	fprintf((data->globalOutput), "%15.6e\t",Tu);
	fprintf((data->globalOutput), "\n");
}

//
//void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data)
//{
//	double *ydata,*psidata;
//  	ydata    = N_VGetArrayPointer_OpenMP(y);
//  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);
//	for (int i = 0; i < data->npts; i++) {
//		fprintf((data->gridOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
//		for (int j = 0; j < data->nvar; j++) {
//			fprintf((data->gridOutput), "%15.6e\t",ydata[j+i*data->nvar]);
//		}
//		fprintf((data->gridOutput), "\n");
//	}
//	fprintf((data->gridOutput), "\n\n");
//}
//
//
////void repairSolution(N_Vector y, N_Vector ydot, UserData data){
////	int npts=data->npts;
////  	double *ydata;
////  	double *ydotdata;
////	ydata    = N_VGetArrayPointer_OpenMP(y); 
////	ydotdata    = N_VGetArrayPointer_OpenMP(ydot); 
////
////	T(2)=T(1);
////	T(npts-1)=T(npts);
////	Tdot(2)=Tdot(1);
////	Tdot(npts-1)=Tdot(npts);
////	for (int k = 1; k <=data->nsp; k++) {
////		    Y(2,k)=Y(1,k);
////		    Y(npts-1,k)=Y(npts,k);
////
////		    Ydot(2,k)=Ydot(1,k);
////		    Ydot(npts-1,k)=Ydot(npts,k);
////	}
////}


/*Following functions are added to derive the characteristic time scale of species*/
void getTimescale(UserData data, N_Vector* y){
	size_t i, k, nsp,npts ; 
	nsp = data->nsp ; 
	npts = data->npts ; 
	double rho, wdot_mole[nsp], wdot_mass[nsp],MW[nsp],concentra[nsp]; 
	//double time_scale[npts*nsp] ;   

	double* ydata ; 
	ydata = N_VGetArrayPointer_OpenMP(*y); //access the data stored in N_Vector* y 
	
	for(i=1; i<= npts;i++){
		setGas(data,ydata,i); //set the gas state at each grid point
		rho = data->gas->density(); //get the averaged density at each grid point, Unit:kg/m^3
		data->gas->getNetProductionRates(wdot_mole) ; //Unit:kmol/(m^3 * s)
		
		for(k=1; k<= nsp; k++){
			MW(k) = data->gas->molecularWeight(k-1) ; //Unit:kg/kmol
		}

		for(k=1; k<= nsp; k++){
			wdot_mass(k) = wdot_mole(k) * MW(k) ; //Unit:kg/(m^3 * s)
		}

		for(k=1; k<= nsp; k++){
			concentra(k) = Y(i,k) * rho ;  //Unit:kg/m^3
		}

		for(k=1;k<= nsp;k++){
			data->time_scale(i,k) = concentra(k)/(wdot_mass(k)+ 1.00e-16) ; 
		}

	}
}


void printTimescaleHeader(UserData data)
{
	fprintf((data->timescaleOutput), "%15s\t","#1");
	for (size_t k = 1; k <=data->nsp+1; k++) {
		fprintf((data->timescaleOutput), "%15lu\t",k+1);
	}
	fprintf((data->timescaleOutput), "%15lu\n",data->nsp+3);


	fprintf((data->timescaleOutput), "%15s\t%15s\t%15s\t","#time","radius","Temp(K)");
	//fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");
	for (size_t k = 1; k <=(data->nsp); k++) {
		fprintf((data->timescaleOutput), "%15s\t",data->gas->speciesName(k-1).c_str());
	}
	//fprintf((data->output), "%15s\t","Pressure(Pa)");
	//fprintf((data->timescaleOutput), "%15s\n",data->gas->speciesName(data->nsp-1).c_str());
	//fprintf((data->output), "%15s\n","Mdot (kg/s)");
	fprintf((data->timescaleOutput), "\n");
}

void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data)
{
	double* ydata; 
	ydata = N_VGetArrayPointer_OpenMP(*y); 

	for(size_t i=1 ; i<= data->npts; i++){
		fprintf(output, "%15.9e\t%15.9e\t",t,R(i));
		fprintf(output, "%15.9e\t",T(i));
        
		for(size_t k=1; k<=data->nsp;k++){
			fprintf(output, "%15.9e\t",data->time_scale(i,k));
		}
		fprintf(output, "\n");
	}
	fprintf(output, "\n");
}

void floorSmallValue(UserData data, N_Vector* y)
{
    double* ydata;
    ydata = N_VGetArrayPointer_OpenMP(*y);
    //double sum = 0.00; 
    size_t k_bath = data->k_bath;
   
	/*Floor small values to zero*/
	for (size_t i = 1; i <=data->npts; i++) {
		for (size_t k = 1; k <=data->nsp; k++) {
			if(fabs(Y(i,k))<=data->massFractionTolerance){
				Y(i,k)=0.0e0;
			}
		}
	}
    
    /*Dump the error to the bath gas*/
	for (size_t i = 1; i <=data->npts; i++) {
        double sum = 0.00;
		for (size_t k = 1; k <=data->nsp; k++) {
			if(k!=k_bath){
				sum = sum + Y(i,k);
			}
		}
        Y(i,k_bath) = ONE-sum;
	}
    

}

void resetTolerance(UserData data, N_Vector* y,N_Vector* atolv)
{
    double* ydata;
    realtype* atolvdata; 
    ydata = N_VGetArrayPointer_OpenMP(*y);
    atolvdata = N_VGetArrayPointer_OpenMP(*atolv);
    double relTol,radTol,tempTol,presTol,massFracTol,bathGasTol,mdotTol;
    double maxT=0.00, deltaT=400.0;
    relTol = 1.0e-03 ; 
    radTol = 1.0e-05 ; 
    tempTol = 1.0e-03 ; 
    presTol = 1.0e-3 ; 
    massFracTol = 1.0e-06; 
    bathGasTol = 1.0e-05 ; 
    mdotTol = 1.0e-10 ; 


    /*Get the maximum Temperature*/
    for (size_t i =1;i <= data->npts;i++){
        if(T(i) > maxT){
            maxT = T(i);
        }
    }

    /*reset the tolerance when maxT > initialTemperature +deltaT*/
    if(maxT >= (data->initialTemperature + deltaT)){
        data->relativeTolerance = relTol;

        for(size_t i =1; i<=data->npts; i++){
            atolT(i) = tempTol; 
            atolR(i) = radTol ; 
            atolP(i) = presTol ;
            atolMdot(i) = mdotTol ; 
            
            for(size_t k =1; k<= data->nsp; k++){
                if(k!=data->k_bath){
                    atolY(i,k) = massFracTol;
                }else{
                    atolY(i,k) = bathGasTol;
                }
            }
        }
    }

}

void getReactions(UserData data,N_Vector* y,FILE* output){
    double Tmax;
    double* ydata; 
    //double deltaT = 400.0;
    //int i = 0;
    size_t nRxns;
    int index;
    nRxns = data->gas->nReactions();
    //DEBUG
    printf("Total Number of Rxns:%zu \n",nRxns);
    double fwdROP[nRxns],revROP[nRxns],netROP[nRxns];
    std::string* rxnArr = new std::string[nRxns];
    /*DEBUG*/
  	printf("Memory Allocation for Rxns Array Succeed!\n");
    ydata = N_VGetArrayPointer_OpenMP(*y);
    Tmax = maxTemperature(ydata,data->nt,data->nvar,data->npts);
    //get the rxns' equation array
    for(size_t ii=0;ii<nRxns;ii++){
        rxnArr[ii]= data->gas->reactionString(ii);
    }
    if(Tmax>=(data->initialTemperature+data->deltaT)){
        index = maxTemperatureIndex(ydata,data->nt,data->nvar,data->npts); 
        setGas(data,ydata,index);
        /*Get forward/reverse/net rate of progress of rxns*/
        data->gas->getRevRatesOfProgress(revROP);
        data->gas->getFwdRatesOfProgress(fwdROP);
        data->gas->getNetRatesOfProgress(netROP);
	    for(size_t j=0 ; j<nRxns; j++){
            fprintf(output,"%30s\t",rxnArr[j]);
            fprintf(output,"%15.9e\t%15.9e\t%15.9e\t\n",fwdROP[j],revROP[j],netROP[j]);
	    }
	    fprintf(output, "\n");
        fclose(output);
     }
     delete[] rxnArr;
}

void getSpecies(UserData data,N_Vector* y,FILE* output){
    double Tmax;
    double* ydata; 
    //double deltaT = 400.0;
    //int i = 0;
    int index;
    double fwdROP[data->nsp],revROP[data->nsp],netROP[data->nsp];
    std::string* spArr = new std::string[data->nsp];
    /*DEBUG*/
  	printf("Memory Allocation for Species Array Succeed!\n");
    ydata = N_VGetArrayPointer_OpenMP(*y);
    Tmax = maxTemperature(ydata,data->nt,data->nvar,data->npts);
    //get the species name array
    for(size_t ii=0;ii<data->nsp;ii++){
        spArr[ii]= data->gas->speciesName(ii);
    }
   
    if(Tmax>=(data->initialTemperature+data->deltaT)){
        //i = i + 1 ;
        index = maxTemperatureIndex(ydata,data->nt,data->nvar,data->npts); 
        setGas(data,ydata,index);
        /*Get forward/reverse/net rate of progress of rxns*/
        data->gas->getDestructionRates(revROP);
        data->gas->getCreationRates(fwdROP);
        data->gas->getNetProductionRates(netROP);
        /*Print data to the pertinent output file*/
	    for(size_t j=0 ; j<data->nsp; j++){
            fprintf(output,"%15s\t",spArr[j]);
            fprintf(output,"%15.9e\t%15.9e\t%15.9e\t\n",fwdROP[j],revROP[j],netROP[j]);
	    }
	    fprintf(output, "\n");
        fclose(output);
    }
    delete[] spArr;
}