#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"

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.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=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);
}

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);
}

int setAlgebraicVariables(N_Vector* id, UserData data){
	double *iddata;
  	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);
	printf("Oxidizer index: %lu\n",data->k_oxidizer);
	printf("Bath gas index:%lu\n",data->k_bath);
	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;
	}	
	Rid(1)=ONE;
	Tid(1)=ZERO;
	Tid(data->npts)=ZERO;
	if(data->constantPressure){
		Pid(data->npts)=ONE;
	}
	else{
		Pid(data->npts)=ZERO;
		Rid(data->npts)=ONE;
	}
	for (size_t k = 1; k <=data->nsp; k++) {
		Yid(1,k)=ZERO;
		Yid(data->npts,k)=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){

	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];

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

	double dx=xOld[nRows-1]/((double)(nPts)-1.0e0);
	for (size_t j = 0; j < nPts; j++) {
		xNew[j]=(double)j*dx;
	}

	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 = 1; k < nPts-1; 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;
	/*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();                                           
		//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);
	}

	/*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;
	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;

	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 Rmin=1e-03*data->domainLength;
	double Rmin=0.0e0;
	double dR=(data->domainLength-Rmin)/((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);
	for (size_t i = 1; i <=data->npts; i++) {
		if(data->metric==0){
			R(i)=Rmin+(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;

	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);
			ier=reGrid(data->grid, x);
			if(ier==-1)return(-1);
			updateSolution(ydata, ydotdata, data->nvar,
			               data->grid->xOld,data->grid->x,data->npts);
			storeGrid(data->grid->x,data->grid->xOld,data->npts);
		}
	}
	else{
		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);
			fclose(input);
		}
		else{
			printf("file initialCondition.dat not found!\n");
			return(-1);
		}
	}
	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:
	if(data->adaptiveGrid){
		data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
				  		data->grid->x, data->npts);
		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);
	}

	/*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;

	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],   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 sum, sum1, sum2, sum3;

	size_t j,k;
	int m;
	m=data->metric;					//Unitless
	mass=data->mass;				//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
	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);                                 
	areamhalfsq= areamhalf*areamhalf;                                   
                                                                             
	/*******************************************************************/

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

	/*Mass:*/                                                            
	Rres(1)=Rdot(1);

	/*Energy:*/                                                     

	if(data->dirichletInner){
		//Tres(1)=Tdot(1);
		Tres(1)=T(1)-data->innerTemperature;
	}
	else{
		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));
	}
                                                                        
	/*Species:*/                                                    
	sum=ZERO;
	for (k = 1; k <=nsp; k++) {                                        
		if(k!=k_bath){
			if(fabs(mdotIn)>1e-14){
		      		Yres(1,k)=innerMassFractionsData[k-1]-
					  Y(1,k)-
					  (YVmhalf(k)*areamhalf)/mdotIn;
			}
			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);
			}
		      sum=sum+Y(1,k);
		}
	}                                                                 
	Yres(1,k_bath)=ONE-sum-Y(1,k_bath);                                     

                                                                             
	/*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;		    
		/**************************************************///     

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

		/*Energy:*/
		/* ∂T/∂t = - ṁ(∂T/∂ψ) 
		 * 	   + (1/cₚ)(∂/∂ψ)(λρ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));
		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))*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);
		advTerm   = (mdotIn*(T(j)-T(j-1))/dpsim);
		Tres(j)= tranTerm
			+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   = (mdotIn*(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;
	}
	else{
		Tres(npts)=T(npts)-T(npts-1);
	}
                                                                
	/*Species:*/                                            
	sum=ZERO;
	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);                     
	}


	/*******************************************************************/

	//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+2; k++) {
		fprintf((data->output), "%15lu\t",k+1);
	}
	fprintf((data->output), "\n");

	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\n","Pressure(Pa)");
}

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.6e\t%15.6e\t",psi(i+1),t);
		if(i==0){
			fprintf(output, "%15.6e\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++) {
			fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
		}
		fprintf(output, "\n");
	}
	fprintf(output, "\n\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);
////	}
////}