/*
sensBVP: A program that calculates the sensitivity of ignition delay time to
kinetic rate parameters by using a boundary value problem formulation.
Copyright (C) 2019  Vyaas Gururajan

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/
#include <stdbool.h>
#include <math.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <time.h>

/*Cantera include files*/
#include "IdealGasMix.h"

/*Sundials include files*/
#include <cvode/cvode.h>               /* prototypes for CVODE fcts., consts.  */
#include <nvector/nvector_serial.h>    /* access to serial N_Vector            */
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix            */
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver      */
#include <cvode/cvode_direct.h>        /* access to CVDls interface            */
#include <sundials/sundials_types.h>   /* defs. of realtype, sunindextype      */
#include <sundials/sundials_math.h>

#include <kinsol/kinsol.h>             /* access to KINSOL func., consts. */
//#include <kinsol/kinsol_impl.h>        /* access to KINSOL data structs */
#include <sunmatrix/sunmatrix_band.h>  /* access to band SUNMatrix        */
#include <kinsol/kinsol_direct.h>      /* access to KINDls interface      */
#include <sunlinsol/sunlinsol_band.h>  /* access to band SUNLinearSolver  */
//#include <sunlinsol/sunlinsol_lapackband.h>  /* access to band SUNLinearSolver  */

#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>

static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
        (imaxarg1) : (imaxarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
        (iminarg1) : (iminarg2))

/*quick-sort related routines*/
static unsigned long *lvector(long nl, long nh);
static void free_lvector(unsigned long *v, long nl, long nh);
static void sort(unsigned long n, double arr[], int ind[]);

#define ZERO  RCONST(0.0)
#define HALF  RCONST(0.5)
#define ONE   RCONST(1.0)
#define TWO   RCONST(2.0)
#define THREE RCONST(3.0)
#define TEN   RCONST(10.0)
#define BUFSIZE 1000

/* In order to keep begin the index numbers from 1 instead of 0, we define
 * macros here. Also, we define macros to ease referencing various variables in
 * the sundials nvector.
 */

#define Y0(k)  NV_Ith_S(data->Y0,k-1)
#define x(i)   NV_Ith_S(data->x,i-1)

#define t(i)   NV_Ith_S(data->tArray,i)
#define P(i)   NV_Ith_S(data->PArray,i)
#define dPdt(i)   NV_Ith_S(data->dPdtArray,i)

#define T(i)   NV_Ith_S(y,((i-1)*data->nvar)+data->nt)
#define Y(i,k) NV_Ith_S(y,((i-1)*data->nvar)+data->ny+k-1)

#define Tdot(i)   NV_Ith_S(ydot,((i-1)*data->nvar)+data->nt)
#define Ydot(i,k) NV_Ith_S(ydot,((i-1)*data->nvar)+data->ny+k-1)

#define Ttemp(i)   NV_Ith_S(ytemp,((i-1)*(data->nvar+1))+data->nt)
#define Ytemp(i,k) NV_Ith_S(ytemp,((i-1)*(data->nvar+1))+data->ny+k-1)
#define tautemp(i) NV_Ith_S(ytemp,((i-1)*(data->nvar+1))+data->ntau)

#define Tp(i)   NV_Ith_S(yp,((i-1)*(data->nvar+1))+data->nt)
#define Yp(i,k) NV_Ith_S(yp,((i-1)*(data->nvar+1))+data->ny+k-1)
#define taup(i) NV_Ith_S(yp,((i-1)*(data->nvar+1))+data->ntau)

#define Pp(i)      NV_Ith_S(data->Pp,i-1)
#define dPdtp(i)   NV_Ith_S(data->dPdtp,i-1)

#define Tres(i)   NV_Ith_S(res,((i-1)*(data->nvar+1))+data->nt)
#define Yres(i,k) NV_Ith_S(res,((i-1)*(data->nvar+1))+data->ny+k-1)
#define taures(i) NV_Ith_S(res,((i-1)*(data->nvar+1))+data->ntau)

#define tautmp1(i) NV_Ith_S(tmp1,((i-1)*(data->nvar+1))+data->ntau)
#define tautmp2(i) NV_Ith_S(tmp2,((i-1)*(data->nvar+1))+data->ntau)

#define wdot(i) wdot[i-1]
#define enthalpy(i) enthalpy[i-1]


typedef struct {
	/* This structure contains all information relevant to evaluating the
	 * residue.
	 */
	Cantera::IdealGasMix* gas;//Ideal gas object containing thermodynamic and
				  //kinetic data
				
	realtype P;		//Pressure (in Pascals)
	realtype T0;		//Initial Temperature (in K)
	N_Vector Y0;		//Initial Mass Fractions
	realtype rho0;		//initial density
	realtype TIgnFac;	//factor that determines the ignition
				//temperature: should be anywhere between 0 and 1 where 0 corresponds to the
		     		//initial temperature and 1 corresponds to adiabatic flame
		     		//temperature
	realtype TIgn;		//Ignition Temperature (in K)
	realtype tauIgn;	//Ignition delay time  (in s)
	N_Vector x;		//The grid
	N_Vector tArray;
	N_Vector PArray;
	N_Vector dPdtArray;
	N_Vector Pp;
	N_Vector dPdtp;
	realtype dPdt;
	int tArraySize;
	bool constantVolume;	//boolean to enable/disable constant volume
				//version of the ignition delay problem
	bool imposedPdt;	//boolean to enable/disable manual entry for P and dPdt
	int jguess;		//saved index to accelerate lookup via "hunt"
	bool writeSpeciesSensitivities;	//boolean to enable/disable writing
				//of species sensitivities
	bool IVPSuccess;
	bool BVPSuccess;
	bool sensSuccess;
	bool solveBVP; 		//The BVP does not HAVE to be solved: use this flag to skip solving it!
	int nsp,neq,npts,nvar;	//key quantities required in for loops:
				//nsp-> no: of chemical species
				//neq-> number of governing equations
				//npts-> no: of grid points
				//nvar->no: of variables
	int nt,ny,ntau;		//array indices for various variables: 
				//nt-> temperature (0)
				//ny-> species mass fractions
				//ntau-> ignition delay time
	int KSneq;		//No: of equations in the BVP formulation of
				//the ignition delay problem
	realtype rel;		//the relative perturbation factor used in
				//computing difference quotients
	realtype atol;
	realtype rtol;
	realtype ftol;
	
  	realtype t1;		//Maximum time to run simulation
	int nreac;		
	int pert_index;		//index of reaction whose rate is perturbed
	bool sensitivityAnalysis;// boolean to activate perturbations for
				//sensitivity analysis
	bool sort;		//boolean to activate sorting
	bool sens;		//if true, perform sensitivity analysis
	FILE *output;		//output file for temperature and species profiles
	FILE *ignSensOutput;	//output file for ignition sensitivities
	FILE *speSensOutput;	//output file for species sensitivities
	FILE *dPdtInput;	//output file for species sensitivities

	gsl_interp_accel* acc;
	gsl_interp_accel* accdot;
	gsl_spline* spline;
	gsl_spline* splinedot;
} *UserData;


/* Set the initial conditions using this subroutine: */
static int setInitialProfile(UserData data, N_Vector y);

/* Evaluate the residue for the IVP using this subroutine: */
static int residue(realtype t, N_Vector y, N_Vector ydot, void *user_data);

/* Evaluate the residue for the BVP using this subroutine: */
static int residueKS(N_Vector yp, N_Vector res, void *user_data);

static int parseInput(UserData data, int argc, char *argv[]);

static int hunt(realtype x, realtype xx[], int n, int jguess);

static void polint(realtype *xdata, realtype *f, int n, realtype x,realtype *y, realtype *dy);

static void lookupDpdt(UserData data, realtype t, realtype *P, realtype *dPdt);

/* Subroutine to compute the Jacobian using numerical differencing: */
static int kinDlsBandDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac, UserData data, N_Vector scale, N_Vector tmp1, N_Vector tmp2);

/* Subroutine that prints the column titles in the output file: */
static void printInstructions();

/* Subroutine that prints the column titles in the output file: */
static void printHeader(UserData data);
static void printSensitivitiesHeader(UserData data);

/* Subroutine that prints the output of the IVP into the output file contained
 * in the UserData structure: */
static void printOutput(realtype t, N_Vector y, UserData data);

/* Subroutine that prints the output of the BVP into the output file contained
 * in the UserData structure: */
//static void printOutputKS(N_Vector yp, UserData data);

/* Subroutine that prints the residue of the BVP into the output file contained
 * in the UserData structure: */
//static void printResidueKS(N_Vector res, UserData data);

/* Subroutine that prints the species sensitivities into the output file
 * (speSensOutput) contained in the UserData structure: */
static void printSpeciesSensitivities(int index, N_Vector yp, N_Vector res, UserData data);

/* Subroutine that prints the sensitivities into the output file (ignSensOutput)
 * contained in the UserData structure: */
static void printIgnSensitivities(realtype sensCoeffs[], int indices[], UserData data);

static int parsedPdt(UserData data);

/* Print solver statistics for the IVP: */
static void PrintFinalStats(void *cvode_mem);

/* Print solver statistics for the BVP: */
static void PrintFinalStatsKS(void *kmem);

/* Subroutine that reports failure if memory hasn't been successfully allocated
 * to various objects: */
static int check_flag(void *flagvalue, const char *funcname, int opt);

int main(int argc, char *argv[])
{

	clock_t start, end;
     	double cpu_time_used;
     	start = clock();

  	int ier;			//error flag

  	UserData data;			//User defined data 
  	data = NULL;

  	/* Create and load problem data block. */
  	data = (UserData) malloc(sizeof *data);
	ier=parseInput(data,argc,argv);
	if(ier==-1)return(-1);

	/* Set the maximum number of time-steps that will be used in the BVP
	 * formulation of the problem: */
	int nout=50000;

	/* Create a temporary solution array to save the results of the IVP for
	 * use in the BVP:*/
  	N_Vector ytemp;
	data->npts=nout;
	data->KSneq=(data->nvar+1)*data->npts;
	ytemp=N_VNew_Serial(data->KSneq);

	{
		/* Solver Memory:*/
  		void *mem;	
		mem=NULL;

		/* Save the initial mass fractions: */
		data->Y0=N_VNew_Serial(data->nsp);	
		for(int k=1;k<=data->nsp;k++){
  			Y0(k)=data->gas->massFraction(k-1); 	
		}

		/*Assign variable indices:*/
		data->nt=0;
		data->ny=data->nt+1;
		data->ntau=data->ny+data->nsp;
		
		/*Get and save the no: of reactions:*/
		data->nreac=data->gas->nReactions();

		/* Solution vector of the IVP: */
  		N_Vector y;
  		y = NULL;
  		y = N_VNew_Serial(data->neq);
  		ier=check_flag((void *)y, "N_VNew_Serial", 0);

  		realtype t0, tret;			//times
		t0=tret=ZERO;

		/*Set the initial values:*/
		setInitialProfile(data, y);
		/*Print out the column names and the initial values:*/
		printHeader(data);
		printOutput(0.0e0,y,data);

  		/* Create a CVode solver object for solution of the IVP:  */
  		//mem = CVodeCreate(CV_BDF,CV_NEWTON);
  		mem = CVodeCreate(CV_BDF);
  		ier=check_flag((void *)mem, "CVodeCreate", 0);
		
		/*Associate the user data with the solver object: */
  		ier = CVodeSetUserData(mem, data);
  		ier=check_flag(&ier, "CVodeSetUserData", 1);

		/* Initialize the solver object by connecting it to the solution vector
		 * y: */
  		ier = CVodeInit(mem, residue, t0, y);
  		ier=check_flag(&ier, "CVodeInit", 1);
		
		/*Set the tolerances: */
  		ier = CVodeSStolerances(mem, data->rtol, data->atol);
  		ier=check_flag(&ier, "IDASStolerances", 1);

      		/* Create dense SUNMatrix for use in linear solves: */
  		SUNMatrix A;
      		A = SUNDenseMatrix(data->neq, data->neq);
      		ier=check_flag((void *)A, "SUNDenseMatrix", 0);

      		/* Create dense SUNLinearSolver object for use by CVode: */
  		SUNLinearSolver LS;
      		LS = SUNDenseLinearSolver(y, A);
      		ier=check_flag((void *)LS, "SUNDenseLinearSolver", 0);

		/* Call CVDlsSetLinearSolver to attach the matrix and linear solver to
		 * CVode: */
      		ier = CVDlsSetLinearSolver(mem, LS, A);
      		ier=check_flag(&ier, "CVDlsSetLinearSolver", 1);

      		/* Use a difference quotient Jacobian: */
      		ier = CVDlsSetJacFn(mem, NULL);
      		ier=check_flag(&ier, "CVDlsSetJacFn", 1);
		
		/*Save the initial values in the temporary solution vector:*/
		Ttemp(1)=T(1);
		for(int k=1;k<=data->nsp;k++){
			Ytemp(1,k)=Y(1,k);
		}
		tautemp(1)=tret;

		/*Begin IVP solution; Save solutions in the temporary solution vector;
		 * stop problem when the temperature reaches a certain value (like 400
		 * K) corresponding to the time for ignition :*/
		int i=1;
		bool BVPCutOff=false;
		int skip=10;
		int count=0;
		while (tret<=data->t1) {
			if(i<data->npts){
				if(!BVPCutOff){
					Ttemp(i+1)=T(1);
					for(int k=1;k<=data->nsp;k++){
						Ytemp(i+1,k)=Y(1,k);
					}
					tautemp(i+1)=tret;
				}
			}
			else{
				printf("Need more points!\n");
				printf("Hard-coded npts=%d not enough!\n",data->npts);
				data->IVPSuccess=false;
				break;
			}
			if(T(1)>=data->TIgn && !BVPCutOff){
				printf("Ignition Delay: %15.6es\n", tret);
				BVPCutOff=true;
			}
  			ier = CVode(mem, data->t1, y, &tret, CV_ONE_STEP);

  			if(check_flag(&ier, "CVode", 1)) {
				data->IVPSuccess=false;
				break;
			}
			else{
				data->IVPSuccess=true;
			}

			printOutput(tret,y,data);

			count++;

			if(tret>1e-06 && !BVPCutOff && count%skip==0){
				i++;
				count=0;
			}

		}

		data->npts=i;
		data->KSneq=(data->nvar+1)*data->npts;

  		/* Print remaining counters and free memory. */
		PrintFinalStats(mem);
  		CVodeFree(&mem);
  		N_VDestroy_Serial(y);
	}
	
	/*Create a solution vector yp, dimensioned such that the maximum no: of
	 * grid points associated with it is the no: of time-steps (i) used
	 * above. Fill in its values using the temporary solution vector. */
  	N_Vector yp;
	if(data->IVPSuccess){
		yp=N_VNew_Serial(data->KSneq);
		for(int j=1;j<=data->npts;j++){
			Tp(j)=Ttemp(j);
			for(int k=1;k<=data->nsp;k++){
				if(fabs(Ytemp(j,k))>1e-16){
					Yp(j,k)=Ytemp(j,k);
				}
				else{
					Yp(j,k)=ZERO;
				}
			}
			taup(j)=tautemp(j);
		}

		data->tauIgn=taup(data->npts);	//Save the ignition delay time
		data->TIgn=Tp(data->npts);	//Save the temperature corresponding to
						//the ignition delay
		/* Create a grid (in time) to be used in the solution of the BVP: */
  		data->x = N_VNew_Serial(data->npts);
  		if(check_flag((void *)data->x, "N_VNew_Serial", 0)) return(1);

		if(data->imposedPdt){
  			data->Pp    = N_VNew_Serial(data->npts);
  			if(check_flag((void *)data->Pp, "N_VNew_Serial", 0)) return(1);
  			data->dPdtp = N_VNew_Serial(data->npts);
  			if(check_flag((void *)data->dPdtp, "N_VNew_Serial", 0)) return(1);
		}

		realtype P,dPdt;
		for (int j = 1; j <=data->npts; j++) {
			if(data->imposedPdt){	
				//P=dPdt=ZERO;
				lookupDpdt(data, taup(j), &P, &dPdt);
				Pp(j)=P*Cantera::OneAtm;
				dPdtp(j)=dPdt*Cantera::OneAtm;
				//printf("%15.6e\t%15.6e\n",Pp(j),dPdtp(j));
			}
			x(j)=taup(j)/data->tauIgn;
			taup(j)=data->tauIgn;
		}
		if(data->tArray!=NULL){
  			N_VDestroy_Serial(data->tArray);
		}
		if(data->PArray!=NULL){
  			N_VDestroy_Serial(data->PArray);
		}
		if(data->dPdtArray!=NULL){
  			N_VDestroy_Serial(data->dPdtArray);
		}
		//printOutputKS(yp,data);
	}

  	N_VDestroy_Serial(ytemp);	//Destroy the temporary solution vector
	printf("No: of time-steps: %d\n",data->npts);
	
	/********************************************************/

	/*Begin solution of the BVP here:*/
	if(data->IVPSuccess && data->sens){
		/*Create a KINSOL solver object for the solution of the BVP:*/
		/* Solver Memory:*/
  		void *mem; mem = NULL;
  		mem = KINCreate();
  		//if (check_flag((void *)mem, "KINCreate", 0)) return(1);

		/*Associate the user data with the solver object: */
  		ier = KINSetUserData(mem, data);
  		ier=check_flag(&ier, "KINSetUserData", 1);

		/* Initialize the solver object by connecting it to the solution vector
		 * yp and the residue "residueKS": */
  		ier = KINInit(mem, residueKS, yp);
  		ier = check_flag(&ier, "KINInit", 1);

  		/* Specify stopping tolerance based on residual: */
  		ier = KINSetFuncNormTol(mem, data->ftol);
  		ier = check_flag(&ier, "KINSetFuncNormTol", 1);
		
		/* Create banded SUNMatrix for use in linear solves; bandwidths are
		 * dictated by the dependence of the solution at a given time on one
		 * time-step ahead and one time-step behind:*/
		SUNMatrix J;
		int mu,ml;
		mu = data->nvar+2; ml = data->nvar+2;
  		//J = SUNBandMatrix(data->KSneq, mu, ml, mu+ml);
  		J = SUNBandMatrix(data->KSneq, mu, ml);
  		ier = check_flag((void *)J, "SUNBandMatrix", 0);

      		/* Create dense SUNLinearSolver object for use by KINSOL: */
  		SUNLinearSolver LSK;
  		LSK = SUNBandLinearSolver(yp, J);
  		ier = check_flag((void *)LSK, "SUNBandLinearSolver", 0);

		/* Call KINDlsSetLinearSolver to attach the matrix and linear solver to
		 * KINSOL: */
  		ier = KINDlsSetLinearSolver(mem, LSK, J);
  		ier = check_flag(&ier, "KINDlsSetLinearSolver", 1);

  		/* No scaling used: */
		N_Vector scale;
		scale=NULL;
  		scale = N_VNew_Serial(data->KSneq);
  		//ier = check_flag((void *)scale, "N_VNew_Serial", 0);
  		N_VConst_Serial(ONE,scale);

  		/* Force a Jacobian re-evaluation every mset iterations: */
  		int mset = 100;
  		ier = KINSetMaxSetupCalls(mem, mset);
  		//ier = check_flag(&ier, "KINSetMaxSetupCalls", 1);

  		/* Every msubset iterations, test if a Jacobian evaluation
  		   is necessary: */
  		int msubset = 1;
  		ier = KINSetMaxSubSetupCalls(mem, msubset);
  		//ier = check_flag(&ier, "KINSetMaxSubSetupCalls", 1);

		ier=KINSetPrintLevel(mem,2);

  		/* Solve the BVP! */
		if(data->solveBVP){
  		  ier = KINSol(mem,     /* KINSol memory block */
  		               yp,      /* initial guess on input; solution vector */
  		               KIN_NONE,//KIN_LINESEARCH,/* global strategy choice */
  		               scale,   /* scaling vector, for the variable cc */
  		               scale);  /* scaling vector for function values fval */

  		  if (check_flag(&ier, "KINSol", 1)){
		  	data->BVPSuccess=false;
		  }else{
		  	data->BVPSuccess=true;
		  	/* Get scaled norm of the system function */
		  	ier = KINGetFuncNorm(mem, &data->ftol);
		  	//ier = check_flag(&ier, "KINGetfuncNorm", 1);
		  	printf("\nComputed solution (||F|| = %g):\n\n",data->ftol);
		  	printf("KinSOL Ignition Delay: %15.6es\n", taup(1));
		  }
		}
  		//ier=check_flag(&ier, "KINSol", 1);

		//printOutputKS(yp,data);
		fclose(data->output);
		PrintFinalStatsKS(mem);
  		KINFree(&mem);
  		N_VDestroy_Serial(scale);
	}

	/********************************************************/
	/*Begin sensitivity analysis here:*/
	if(data->BVPSuccess || !data->solveBVP){
		/* Create banded SUNMatrix for use in linear solves; bandwidths
		 * are dictated by the dependence of the solution at a given
		 * time on one time-step ahead and one time-step behind:*/
		SUNMatrix J;
		int mu,ml;
		mu = data->nvar+1; ml = data->nvar+1;
  		//J = SUNBandMatrix(data->KSneq, mu, ml, mu+ml);
  		J = SUNBandMatrix(data->KSneq, mu, ml);
  		ier = check_flag((void *)J, "SUNBandMatrix", 0);

		/*Create a residue vector and two temporary vectors:*/
  		N_Vector res,tmp1,tmp2, scale;
		res=N_VNew_Serial(data->KSneq);
		tmp1=N_VNew_Serial(data->KSneq);
		tmp2=N_VNew_Serial(data->KSneq);
  		scale = N_VNew_Serial(data->KSneq);
  		N_VConst_Serial(ONE,scale);

		/*Evaluate the residue using the solution computed by solving
		 * the BVP above:*/
		ier = residueKS(yp, res, data);
		ier = check_flag(&ier, "residueKS", 1);
		//printResidueKS(res,data);
       
       		/*Compute the Jacobian and store it in J:*/	
		kinDlsBandDQJac(yp, res, J, data, scale, tmp1, tmp2);

		/*Create a linear solver object for the banded system; in the
		 * problem Ax=b, A is J and x is tmp2:*/
  		SUNLinearSolver LS;
  		LS = SUNBandLinearSolver(res, J);

		/*Initialize the solver and perform an LU factorization by
		 * running SUNLinSolSetup:*/
		ier = SUNLinSolInitialize_Band(LS);
		ier = check_flag(&ier, "SUNLinSolInitialize", 1);
		ier = SUNLinSolSetup_Band(LS,J);
		ier = check_flag(&ier, "SUNLinSolSetup", 1);

		/*Create an array to store the logarithmic sensitivity
		 * coefficients:*/
		realtype sensCoeffs[data->nreac];

		/*Create an array to store the indices of the reactions (needed
		 * for sorting the sensitivity coefficients):*/
		int indices[data->nreac];

		/*The sensitivities are computed as follows:
		 * For a system equations
		 *    F(y;α)=0
		 *    where F is the residue of the BVP
		 *    	    y is the solution of the BVP
		 *    	    α is a parameter 
		 * => (∂F/∂y)(∂y/∂α)+(∂F/∂α) = 0
		 * => (∂F/∂y)(∂y/∂α)=-(∂F/∂α)
		 * Therefore, the sensitivity coefficient (∂y/∂α) can be
		 * computed by solving a system Ax=b, where A=(∂F/∂y),
		 * x=(∂y/∂α), and b=-(∂F/∂α)!*/

		/*Run a loop over all the reactions: */
		printf("\nRunning Sensitivity Analysis...\n");
		/*Enable sensitivity analysis:*/
		data->sensitivityAnalysis=true;
		realtype oneOverRel=ONE/(data->rel);
		for(int j=0;j<data->nreac;j++){

			/*Save the index of the reaction whose collision
			 * frequency (A in the Arrhenius law k=A*exp(-Eₐ/RT))
			 * is to be perturbed:*/ 
			data->pert_index=j;	

			/*Compute the perturbed residue and store it in tmp1:*/
			ier=residueKS(yp, tmp1, data);
			ier=check_flag(&ier, "residueKS", 1);
			//printResidueKS(tmp1,data);

			/*Compute the difference F(α+Δα)-F(α) and store it in
			 * tmp2:*/
			N_VLinearSum(ONE,tmp1,-ONE,res,tmp2);

			/*Divide the difference quotient by -Δα and store the
			 * result in tmp1. This is the numerical approximation
			 * to -(∂F/∂α):*/
			N_VScale(-oneOverRel,tmp2,tmp1);

			/*Solve the linear system (thankfully the LU
			 * factiorization has already been carried out once and
			 * for all above!) and store the solution in tmp2:*/
			ier=SUNLinSolSolve_Band(LS,J,tmp2,tmp1,ZERO);
			ier=check_flag(&ier, "SUNLinSolSolve", 1);

			/*Divide the result by the ignition delay time to get
			 * the logarithmic sensitivity coefficient and save it
			 * in sensCoeffs:*/
			sensCoeffs[j]=tautmp2(1)/taup(1);

			if(sensCoeffs[j]!=sensCoeffs[j]){
				printf("\nNaN! Quitting Program! Try different tolerances!\n");
				data->sensSuccess=false;
				break;
			}else{
				data->sensSuccess=true;
				printf("Reaction %5d: %15.6e\n",j,sensCoeffs[j]);
			}

			/*Print out the temperature and species mass-fraction
			 * sensitivities:*/
			if(data->writeSpeciesSensitivities){
				printSpeciesSensitivities(j, yp, tmp2, data);
			}

			/*Save the index of the reaction:*/
			indices[j]=j;
		}

		if(data->sensSuccess){

		  /*Sort the sensitivities in ascending order. Note the advancing
		 * of the beginning indices of the sensCoeffs and indices
		 * arrays. This is due to the Numerical recipes convention for
		 * array indexing used in sort subroutine:*/

		  if(data->sort){
			sort(data->nreac,sensCoeffs-1,indices-1);
		  }
		  /*Print out the sensitivities:*/
		  printIgnSensitivities(sensCoeffs,indices,data);
		}

  		N_VDestroy_Serial(scale);
  		N_VDestroy_Serial(res);
  		N_VDestroy_Serial(tmp1);
  		N_VDestroy_Serial(tmp2);
		fclose(data->ignSensOutput);
		fclose(data->speSensOutput);
	}
	
	/*Free memory and delete all the vectors and user data:*/
	if(yp!=NULL){
  		N_VDestroy_Serial(yp);
		printf("BVP Solution Vector Deleted!\n");
	}
	if(data->Y0!=NULL){
  		N_VDestroy_Serial(data->Y0);
		printf("Initial Mass fraction Vector Deleted!\n");
	}
	if(data->x!=NULL){
  		N_VDestroy_Serial(data->x);
		printf("Grid for BVP deleted!\n");
	}
	if(data->Pp!=NULL){
  		N_VDestroy_Serial(data->Pp);
		printf("P array deleted!\n");
	}
	if(data->dPdtp!=NULL){
  		N_VDestroy_Serial(data->dPdtp);
		printf("dPdt array deleted!\n");
	}
	if(data->gas!=NULL){
		delete data->gas;
		printf("Gas Deleted!\n");

	}
	if(data->imposedPdt){
		if(data->acc!=NULL){
			gsl_interp_accel_free(data->acc);
		}
		if(data->accdot!=NULL){
			gsl_interp_accel_free(data->accdot);
		}
		if(data->spline!=NULL){
			gsl_spline_free(data->spline);
		}
		if(data->splinedot!=NULL){
			gsl_spline_free(data->splinedot);
		}
	}
	if(data!=NULL){
		/* Free the user data */
  		free(data);             
		printf("User data structure deleted!\n");
	}


 	end = clock();
     	cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
	printf("Elapsed cpu time: %15.6e s\n", cpu_time_used);

	return(0);
}

static int parseInput(UserData data, int argc, char *argv[]){

	/*Set defaults here:*/
	/*****************************************************/
	/*Relative tolerance for IVP:*/
  	data->rtol = RCONST(1e-14);
	/*Absolute tolerance for IVP:*/
  	data->atol = RCONST(1e-14);
	/*Residue tolerance for BVP:*/
	data->ftol = RCONST(1e-13);
	/*Final Time:*/
	data->t1   = RCONST(10.0);
	/*Solve constant pressure problem:*/
	data->constantVolume=false;	
	/*Do not impose a P vs t curve:*/
	data->imposedPdt=false;		
	/*Disable writing species sensitivities:*/
	data->writeSpeciesSensitivities=false;
	/*Enable sensitivity analysis:*/
	data->sens=false;
	/*Index to start with for pressure lookup (see hunt):*/
	data->jguess=0;			
	/*Set rate of change of pressure to zero:*/
	data->dPdt=ZERO;
	/*Disable sensitivity analysis:*/
	data->sensitivityAnalysis=false;
	/*Find the relative perturbation constant:*/
	data->rel=SUNRsqrt(UNIT_ROUNDOFF);
	/*Set flags that indicate success of various stages:*/
	data->IVPSuccess=false;
	data->BVPSuccess=false;
	data->sensSuccess=false;
	/*TIgnFac:*/
	data->TIgnFac=0.20;
	/*Sorting*/
	data->sort=false;
	/* Solve two-point BVP: */
	data->solveBVP=true;
	/*****************************************************/

	int ier;
	int opt;
	char mech[BUFSIZE+1];
	char comp[BUFSIZE+1];
	bool enteredT0, enteredP, enteredMech, enteredComp;
	enteredT0=enteredP=enteredMech=enteredComp=false;
	while((opt=getopt(argc,argv,"a:r:f:T:P:m:c:t:q:vsodn")) != -1){
		switch(opt){
			case 'a':
				data->atol=RCONST(atof(optarg));
				break;
			case 'r':
				data->rtol=RCONST(atof(optarg));
				break;
			case 'f':
				data->ftol=RCONST(atof(optarg));
				break;
			case 'T':
				data->T0=RCONST(atof(optarg));
				enteredT0=true;
				break;
			case 'P':
				data->P=RCONST(atof(optarg))*Cantera::OneAtm; 
				enteredP=true;
				break;
			case 'm':
				snprintf(mech,BUFSIZE,"%s",optarg);
				enteredMech=true;
				break;
			case 'c':
				snprintf(comp,BUFSIZE,"%s",optarg);
				enteredComp=true;
				break;
			case 't':
				data->t1=RCONST(atof(optarg));
				enteredT0=true;
				break;
			case 'v':
				data->constantVolume=true;
				break;
			case 'S':
				data->writeSpeciesSensitivities=true;
				break;
			case 's':
				data->sens=true;
				break;
			case 'o':
				data->sort=true;
				break;
			case 'd':
				data->imposedPdt=true;
				break;
			case 'q':
				data->TIgnFac=RCONST(atof(optarg));
				break;
			case 'n':
				data->solveBVP=false;
				break;
			default:
				printInstructions();
				return(-1);
		}
	}
	if(!enteredT0){
		printf("Not specified initial Temperature! Exiting...\n");
		printInstructions();
		return(-1);
	}
	else if(!enteredP){
		printf("Not specified Pressure! Exiting...\n");
		printInstructions();
		return(-1);
	}
	else if(!enteredMech){
		printf("Not specified Mechanism! Exiting...\n");
		printInstructions();
		return(-1);
	}
	else if(!enteredComp){
		printf("Not specified Composition! Exiting...\n");
		printInstructions();
		return(-1);
	}
	else{
		if(data->imposedPdt){
			ier=parsedPdt(data);
			if(ier==-1){
				return(-1);
			}
		}

		printf("Required inputs provided:\n"); 
		printf("\natol: %15.6e\n", data->atol); 
		printf("\nrtol: %15.6e\n", data->rtol); 
		printf("\nftol: %15.6e\n", data->ftol); 
		printf("\nT0  : %15.6e K\n", data->T0); 
		printf("\nP   : %15.6e Pa\n", data->P); 
		printf("\nmech: %s\n", mech); 
		printf("\ncomp: %s\n", comp); 
		printf("\nt1  : %15.6e s\n", data->t1); 
		printf("\nconstantVolume  : %d\n", data->constantVolume); 
		printf("\nimposedPdt  : %d\n", data->imposedPdt); 
		printf("\nwriteSpeciesSensitivities  : %d\n", 
				data->writeSpeciesSensitivities); 
		printf("Proceeding...\n\n");
	}

	/*Define Gas here:*/
	data->gas = new Cantera::IdealGasMix(mech);

	/* Set the initial state of the gas: */
	data->gas->setState_TPX(data->T0,
				data->P,
				comp);

  	/* Create output file for the solution of the IVP: */
	data->output=fopen("output.dat","w");

	/* Create output file for the ignition sensitivities (post BVP): */
	data->ignSensOutput=fopen("ignitionSensitivities.dat","w");

	/* Create output file for the species sensitivities (post BVP): */
	data->speSensOutput=fopen("speciesSensitivities.dat","w");

	data->rho0=data->gas->density();
	if(data->constantVolume){
		data->gas->equilibrate("UV");
	}
	else{
		data->gas->equilibrate("HP");
	}
	data->TIgn=data->T0+(data->gas->temperature()
			    -data->T0)*data->TIgnFac;
	printf("Ignition Temperature: %15.6e K\n",data->TIgn);
	data->gas->setState_TPX(data->T0,
				data->P,
				comp);

	/*Get and save the no: of species:*/
	data->nsp=data->gas->nSpecies();	
	/*set the no: of variables for the IVP:*/
	data->nvar=data->nsp+1;			 
	/*set the no: of equations for the IVP:*/
	data->neq=data->nvar;			

	return(0);
}

static int setInitialProfile(UserData data, N_Vector y)
{
	/*This routine sets the initial temperature and mass fractions for the
	 * solution of the initial value problem:*/
  	N_VConst(ZERO, y);

	T(1)=data->gas->temperature();
	for (int k = 1; k <=data->nsp; k++) {
  		Y(1,k)=Y0(k);
	}

  	return(0);

}

static int residue(realtype t, N_Vector y, N_Vector ydot, void *user_data)
{

	/*This is the sub-routine that computes the right hand side F(y) in the
	 * problem ∂y/∂t = F(y). 
	 *
	 * The energy conservation equation is 
	 * 	  ∂T/∂t = (1/ρcₚ)(dP/dt -∑ ώᵢhᵢ)
	 *
	 * T has units of K
	 * 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.
	 * ρ has units of kg/m³
	 * cₚ has units J/kg/K
	 * ώᵢ has units of kmol/m³/s
	 * dP/dt has units of Pa/s
	 *
	 * In the case of a constant volume problem, the energy conservation
	 * equation can be re-written as
	 * 	  ∂T/∂t = -(1/ρcᵥ)∑ ώᵢeᵢ
	 * eᵢ has units J/kmol
	 * cᵥ has units J/kg/K
	 *
	 *
	 * The species conservation equation is
	 * 	∂Yᵢ/∂t = ώᵢWᵢ/ρ
	 * Wᵢ has units of kg/kmol
	 */

	/*Assign data structure type to user_data:*/
	UserData data;
	data=(UserData)user_data;

	/*Get no: of species:*/
	int nsp=data->nsp;

	/*Create temporary arrays to store enthalpy (non-dimensional) and ώᵢ:*/
	realtype wdot[nsp];
	realtype enthalpy[nsp];

    	try {
		/*Set the gas conditions:*/
		if(data->constantVolume){
  			data->gas->setMassFractions_NoNorm(&Y(1,1)); 	
  			data->gas->setTemperature(T(1)); 	        
  			data->gas->setDensity(data->rho0);
			data->P=data->gas->pressure();
		}
		else if(data->imposedPdt){
  			data->gas->setMassFractions_NoNorm(&Y(1,1)); 	
  			data->gas->setTemperature(T(1)); 	        
			realtype P,dPdt;
			P=dPdt=0.0e0;
			lookupDpdt(data, t, &P, &dPdt);
			data->P=P*Cantera::OneAtm;
			data->dPdt=dPdt*Cantera::OneAtm;
  			data->gas->setPressure(data->P);
		}
		else{
  			data->gas->setMassFractions_NoNorm(&Y(1,1)); 	
  			data->gas->setTemperature(T(1)); 	        
  			data->gas->setPressure(data->P);
		}

		realtype rho=data->gas->density();	//get ρ
		realtype cp=data->gas->cp_mass();	//get cₚ
		data->gas->getNetProductionRates(wdot);	//get ώᵢ
		data->gas->getEnthalpy_RT(enthalpy);	//get hᵢ/RT

		if(data->constantVolume){
			for(int k=1;k<=nsp;k++){
				enthalpy(k)=enthalpy(k)-ONE;
			}
			cp=data->gas->cv_mass();
		}

		realtype sum=ZERO;
		for(int k=1;k<=nsp;k++){
			Ydot(1,k)=wdot(k)*(data->gas->molecularWeight(k-1))/rho;
			sum=sum+wdot(k)*enthalpy(k)*T(1);
		}
		sum=sum*Cantera::GasConstant;
		Tdot(1)=(data->dPdt-sum)/(rho*cp);

    	} catch (Cantera::CanteraError& err) {
	    printf("Error:\n");
	    printf("%s\n",err.what());
	    return(-1);
    	}

	return(0);
}

static int residueKS(N_Vector yp, N_Vector res, void *user_data)
{

	/*This is the sub-routine that computes the right hand side F(y) in the
	 * problem F(y) = 0. Here the problem is specifically that of a
	 * homogeneous isobaric batch reactor formulated as a boundary value
	 * problem (BVP). The unknowns in this problem are T, Yᵢ, and τ.
	 *
	 * The energy conservation equation is 
	 * 	  ∂T/∂x = (τ/ρcₚ)(dP/dt - ∑ ώᵢhᵢ)
	 *
	 * T has units of K
	 * x is *unitless*
	 * τ has units of time
	 * 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.
	 * ρ has units of kg/m³
	 * cₚ has units J/kg/K
	 * ώᵢ has units of kmol/m³/s
	 * dP/dt has units of Pa/s
	 *
	 * In the case of a constant volume problem, the energy conservation
	 * equation can be re-written as
	 * 	  ∂T/∂x = -(τ/ρcᵥ)∑ ώᵢeᵢ
	 * eᵢ has units J/kmol
	 * cᵥ has units J/kg/K
	 *
	 * The species conservation equation is
	 * 	∂Yᵢ/∂x = τώᵢWᵢ/ρ
	 *
	 * Wᵢ has units of kg/kmol
	 *
	 * The equation (trivial) to maintain a banded structure in the linear
	 * solver is 
	 * 	∂τ/∂x = 0
	 */
	/*Assign data structure type to user_data:*/
	UserData data;
	data=(UserData)user_data;

	/*Get no: of species:*/
	int nsp=data->nsp;
	
	/*Create temporary arrays to store enthalpy (non-dimensional) and ώᵢ:*/
	realtype wdot[nsp];
	realtype enthalpy[nsp];
	
	/*Create temperorary variables to store grid spacing and summation
	 * terms:*/
	realtype dx;
	realtype sum=ZERO;
	data->dPdt=ZERO;

    	try {
		/*If sensitivity analysis has been enabled, perturb the collision
		 * frequency by using cantera's setMultiplier function:*/
		if(data->sensitivityAnalysis){
			data->gas->setMultiplier(data->pert_index,ONE+data->rel);
			//printf("%d\t%15.9e\n",data->pert_index,ONE+data->rel);
		}

		Tres(1)=Tp(1)-data->T0;
		for(int k=1;k<=nsp;k++){
			Yres(1,k)=Yp(1,k)-Y0(k);
		}
		taures(1)=taup(2)-taup(1);


		for(int i=2;i<=data->npts;i++){

			if(data->constantVolume){
  				data->gas->setMassFractions_NoNorm(&Yp(i,1)); 	
  				data->gas->setTemperature(Tp(i)); 	        
  				data->gas->setDensity(data->rho0);
			}
			else if(data->imposedPdt){
				//printf("%15.6e\t%15.6e\t%15.6e\n",x(i),Tp(i),Pp(i));
  				data->gas->setMassFractions_NoNorm(&Yp(i,1)); 	
  				data->gas->setTemperature(Tp(i)); 	        
  				data->gas->setPressure(Pp(i));

				data->dPdt=dPdtp(i);
			}
			else{
  				data->gas->setMassFractions_NoNorm(&Yp(i,1)); 	
  				data->gas->setTemperature(Tp(i)); 	        
  				data->gas->setPressure(data->P);
			}

			realtype rho=data->gas->density();
			realtype cp=data->gas->cp_mass();
			data->gas->getNetProductionRates(wdot);
			data->gas->getEnthalpy_RT(enthalpy);

			if(data->constantVolume){
				for(int k=1;k<=nsp;k++){
					enthalpy(k)=enthalpy(k)-ONE;
				}
				cp=data->gas->cv_mass();
			}

			dx=x(i)-x(i-1);
			sum=ZERO;

			for(int k=1;k<=nsp;k++){
				Yres(i,k)=(Yp(i,k)-Yp(i-1,k))/dx;
				Yres(i,k)+=-taup(i)*wdot(k)*(data->gas->molecularWeight(k-1))/rho;
				sum+=wdot(k)*enthalpy(k)*Tp(i);
			}
			sum=sum*Cantera::GasConstant;
			Tres(i)=(Tp(i)-Tp(i-1))/dx;
			Tres(i)+=taup(i)*(sum-data->dPdt)/(rho*cp);
			if(i==data->npts){
				taures(i)=Tp(i)-data->TIgn;
			}
			else{
				taures(i)=taup(i+1)-taup(i);
			}

		}

		/*If sensitivity analysis has been enabled, *un*perturb the collision
		 * frequency:*/
		if(data->sensitivityAnalysis){
			data->gas->setMultiplier(data->pert_index,ONE);
		}

    	} catch (Cantera::CanteraError& err) {
	    printf("Error:\n");
	    printf("%s\n",err.what());
	    return(-1);
    	}

	return(0);

}

static void lookupDpdt(UserData data, realtype t, realtype *P, realtype *dPdt)
{

	realtype *dPdtData;
	dPdtData = N_VGetArrayPointer_Serial(data->dPdtArray);

	int npts=data->tArraySize;

	if(t<=data->t1){
	
		*P=gsl_spline_eval(data->spline,t,data->acc);
		*dPdt=gsl_spline_eval(data->splinedot,t,data->accdot);
	}
	else{
		*P=P(npts-1);
		*dPdt=dPdt(npts-1);
	}
}

/*------------------------------------------------------------------
  kinDlsBandDQJac

  This routine generates a banded difference quotient approximation
  to the Jacobian of F(u).  It assumes a SUNBandMatrix input stored
  column-wise, and that elements within each column are contiguous.
  This makes it possible to get the address of a column of J via the
  function SUNBandMatrix_Column() and to write a simple for loop to
  set each of the elements of a column in succession.
 
  NOTE: Any type of failure of the system function her leads to an
        unrecoverable failure of the Jacobian function and thus of
        the linear solver setup function, stopping KINSOL.
  ------------------------------------------------------------------*/
static int kinDlsBandDQJac(N_Vector u, N_Vector fu, SUNMatrix Jac,
                    UserData data, N_Vector scale, N_Vector tmp1, N_Vector tmp2)
{
	realtype inc, inc_inv;
	N_Vector futemp, utemp;
	sunindextype group, i, j, width, ngroups, i1, i2;
	sunindextype N, mupper, mlower;
	realtype *col_j, *fu_data, *futemp_data, *u_data, *utemp_data, *uscale_data;
	int retval = 0;
	//KINDlsMem kindls_mem;
	
	//KINMem kin_mem;
	//kin_mem = (KINMem) mem;
	
	///* access DlsMem interface structure */
	//kindls_mem = (KINDlsMem) kin_mem->kin_lmem;
	
	/* access matrix dimensions */
	N = SUNBandMatrix_Columns(Jac);
	mupper = SUNBandMatrix_UpperBandwidth(Jac);
	mlower = SUNBandMatrix_LowerBandwidth(Jac);
	
	/* Rename work vectors for use as temporary values of u and fu */
	futemp = tmp1;
	utemp  = tmp2;
	
	/* Obtain pointers to the data for ewt, fy, futemp, y, ytemp */
	fu_data     = N_VGetArrayPointer(fu);
	futemp_data = N_VGetArrayPointer(futemp);
	u_data      = N_VGetArrayPointer(u);
	//uscale_data = N_VGetArrayPointer(kin_mem->kin_uscale);
	uscale_data = N_VGetArrayPointer(scale);
	utemp_data  = N_VGetArrayPointer(utemp);
	
	/* Load utemp with u */
	N_VScale(ONE, u, utemp);
	
	/* Set bandwidth and number of column groups for band differencing */
	width   = mlower + mupper + 1;
	ngroups = SUNMIN(width, N);
	
	//UserData data;
	//data=(UserData) kin_mem->kin_user_data;
	
	for (group=1; group <= ngroups; group++) {
	  
	  /* Increment all utemp components in group */
	  for(j=group-1; j < N; j+=width) {
	    //inc = kin_mem->kin_sqrt_relfunc*SUNMAX(SUNRabs(u_data[j]),
	    //                                       ONE/SUNRabs(uscale_data[j]));
	    //inc = data->rel*SUNRabs(u_data[j]);
	    inc = data->rel*SUNMAX(SUNRabs(u_data[j]),
	                           ONE/SUNRabs(uscale_data[j]));
	    utemp_data[j] += inc;
	  }
	
	  /* Evaluate f with incremented u */
	  //retval = kin_mem->kin_func(utemp, futemp, kin_mem->kin_user_data);
	  retval=residueKS(utemp, futemp, data);
	  if (retval != 0) return(retval); 
	
	  /* Restore utemp components, then form and load difference quotients */
	  for (j=group-1; j < N; j+=width) {
	    utemp_data[j] = u_data[j];
	    col_j = SUNBandMatrix_Column(Jac, j);
	    //inc = kin_mem->kin_sqrt_relfunc*SUNMAX(SUNRabs(u_data[j]),
	    //                                       ONE/SUNRabs(uscale_data[j]));
	    //inc = kin_mem->kin_sqrt_relfunc*SUNRabs(u_data[j]);
	    //inc = data->rel*SUNRabs(u_data[j]);
	    inc = data->rel*SUNMAX(SUNRabs(u_data[j]),
	                           ONE/SUNRabs(uscale_data[j]));
	    inc_inv = ONE/inc;
	    i1 = SUNMAX(0, j-mupper);
	    i2 = SUNMIN(j+mlower, N-1);
	    for (i=i1; i <= i2; i++)
	      SM_COLUMN_ELEMENT_B(col_j,i,j) = inc_inv * (futemp_data[i] - fu_data[i]);
	  }
	}
	
	///* Increment counter nfeDQ */
	//kindls_mem->nfeDQ += ngroups;
	
	
	return(0);
}

static int hunt(realtype x, realtype xx[], int n, int jguess)
{
	int jlo=jguess;
	int jm,jhi,inc;
	int ascnd;

	ascnd=(xx[n-1] >= xx[0]);
	if (jlo <= 0 || jlo > n-1) {
		jlo=0;
		jhi=n;
	} else {
		inc=1;
		if ((x >= xx[jlo]) == ascnd) {
			if (jlo == n-1) return(jlo);
			jhi=(jlo)+1;
			while ((x >= xx[jhi]) == ascnd) {
				jlo=jhi;
				inc += inc;
				jhi=(jlo)+inc;
				if (jhi > n-1) {
					jhi=n;
					break;
				}
			}
		} else {
			if (jlo == 1) {
				jlo=0;
				return (jlo);
			}
			jhi=(jlo)--;
			while ((x < xx[jlo]) == ascnd) {
				jhi=(jlo);
				inc <<= 1;
				if (inc >= jhi) {
					jlo=0;
					break;
				}
				else jlo=jhi-inc;
			}
		}
	}
	while (jhi-(jlo) != 1) {
		jm=(jhi+(jlo)) >> 1;
		if ((x >= xx[jm]) == ascnd)
			jlo=jm;
		else
			jhi=jm;
	}
	if (x == xx[n-1]) jlo=n-2;
	if (x == xx[0]) jlo=1;

	return(jlo);
}

static void polint(realtype *xdata, realtype *f, int n, realtype x, realtype *y, realtype *dy){
	int i,m,ns=1;
	realtype den,dif,dift,ho,hp,w;
	realtype c[n+1],d[n+1];
	dif=fabs(x-xdata[1]);
	for (i=1;i<=n;i++) {
		if ( (dift=fabs(x-xdata[i])) < dif) {
			ns=i;
			dif=dift;
		}
		c[i]=f[i];
		d[i]=f[i];
	}
	*y=f[ns--];
	for (m=1;m<n;m++) {
		for (i=1;i<=n-m;i++) {
			ho=xdata[i]-x;
			hp=xdata[i+m]-x;
			w=c[i+1]-d[i];
			if ( (den=ho-hp) == 0.0) printf("Error in routine polint!\n");
			den=w/den;
			d[i]=hp*den;
			c[i]=ho*den;
		}
		*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
	}
}

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50
#define NR_END 1
#define FREE_ARG char*

static unsigned long *lvector(long nl, long nh)
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
	unsigned long *v;

	v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
	if (!v) printf("allocation failure in lvector()");
	return v-nl+NR_END;
}

static void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

static void sort(unsigned long n, double arr[], int ind[])
{
	unsigned long i,ir=n,j,k,l=1,*istack;
	int jstack=0;
	double a,temp;
	int b;

	istack=lvector(1,NSTACK);
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				b=ind[j];
				for (i=j-1;i>=l;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
					ind[i+1]=ind[i];
				}
				arr[i+1]=a;
				ind[i+1]=b;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP(arr[k],arr[l+1])
			SWAP(ind[k],ind[l+1])
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
				SWAP(ind[l],ind[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
				SWAP(ind[l+1],ind[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP(arr[l],arr[l+1])
				SWAP(ind[l],ind[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			b=ind[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j]);
				SWAP(ind[i],ind[j]);
			}
			arr[l+1]=arr[j];
			ind[l+1]=ind[j];
			arr[j]=a;
			ind[j]=b;
			jstack += 2;
			if (jstack > NSTACK) printf("NSTACK too small in sort.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	free_lvector(istack,1,NSTACK);
}
#undef M
#undef NSTACK
#undef SWAP

static void printInstructions(){
	printf("\nInputs either incomplete or wrong!\n");
	printf("\n-a <absolute tolerance for IVP>\n");
	printf("\n-r <relative tolerance for IVP>\n");
	printf("\n-f <function tolerance for BVP>\n");
	printf("\n-T <Initial Temperature in Kelvin>\n");
	printf("\n-P <Initial Pressure in atm>\n");
	printf("\n-m <mechanism file (cti or xml)>\n");
	printf("\n-c <composition in mole fractions>\n");
	printf("\n-t <Total simulation time in seconds>\n");
	printf("\n-v :Enables constant volume simulation\n");
	printf("\n-s :Enables ignition delay sensitivity output\n");
	printf("\n-S :Enables species mass fraction sensitivity output (works only with -s)\n");
	printf("\n-d :Enables manual dPdt entryspecies sensitivity output\n");
	printf("\n-o :Enables sorting of ignition delay sensitivity output\n");
	printf("\n-q :factor that determines ignition temperature (default: 0.2)\n");
	printf("\n-n :Disables the boundary value problem solver\n");
	printf("\nexample: <executable> ");
	printf("-T 1200.0 -P 1.0 -m gri30.cti");
	printf(" -c H2:1.0,N2:3.76,O2:1.0");
	printf(" -t 1e-03\n");
}

static void printHeader(UserData data)
{
	fprintf((data->output), "%15s\t","#1");
	for (int k = 1; k <=data->nvar; k++) {
		fprintf((data->output), "%15d\t",k+1);
	}
	fprintf((data->output), "\n");

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

}

static void printSensitivitiesHeader(UserData data)
{
	fprintf((data->speSensOutput), "%15s\t","#1");
	for (int k = 1; k <=data->nvar; k++) {
		fprintf((data->speSensOutput), "%15d\t",k+1);
	}
	fprintf((data->speSensOutput), "\n");

	fprintf((data->speSensOutput), "%15s\t%15s\t",
			"#time(s)","Temp(K)");
	for (int k = 1; k <=data->nsp; k++) {
		fprintf((data->speSensOutput), "%15s\t",
				data->gas->speciesName(k-1).c_str());
	}
	fprintf((data->speSensOutput), "\n");

}

static void printOutput(realtype t, N_Vector y, UserData data)
{
	fprintf((data->output), "%15.6e\t%15.6e\t%15.6e\t",t,T(1),data->P);
	for (int k = 1; k <=data->nsp; k++) {
		fprintf((data->output), "%15.6e\t",Y(1,k));
	}
	fprintf((data->output), "\n");
}


//static void printOutputKS(N_Vector yp, UserData data)
//{
//	//printf("%d\n",data->npts);
//	fprintf(data->output, "\n\n");
//	for(int i=1;i<=data->npts;i++){
//		fprintf((data->output), "%15.6e\t%15.6e\t%15.6e\t",x(i)*data->tauIgn,Tp(i),data->P);
//		//printf("%15.6e\t%15.6e\n",x(i),Tp(i));
//		for (int k = 1; k <=data->nsp; k++) {
//			fprintf((data->output), "%15.6e\t",Yp(i,k));
//		}
//		fprintf(data->output, "\n");
//	}
//}

static void printSpeciesSensitivities(int index, N_Vector yp, N_Vector res, UserData data)
{
	char buf[50];
	std::string line;
	line=data->gas->reactionString(index);
	sprintf(buf,"%s",line.c_str());
	fprintf((data->speSensOutput), "#%d: %38s\t\n", index, buf);
	printSensitivitiesHeader(data);
	for(int i=1;i<=data->npts;i++){
		fprintf((data->speSensOutput), "%15.6e\t%15.6e\t",x(i)*taup(i),Tres(i)/Tp(i));
		for (int k = 1; k <=data->nsp; k++) {
			if(Yp(i,k)>=data->atol){
				fprintf((data->speSensOutput), "%15.6e\t",Yres(i,k)/(Yp(i,k)+1e-31));
			}
			else{
				fprintf((data->speSensOutput), "%15.6e\t",ZERO);
			}
		}
		fprintf((data->speSensOutput), "\n");
	}
	fprintf((data->speSensOutput), "\n\n");
}

//static void printResidueKS(N_Vector res, UserData data)
//{
//	fprintf((data->output), "\n\n");
//	for(int i=1;i<=data->npts;i++){
//		fprintf((data->output), "%15.6e\t%15.6e\t",x(i),Tres(i));
//		for (int k = 1; k <=data->nsp; k++) {
//			fprintf((data->output), "%15.6e\t",Yres(i,k));
//		}
//		fprintf((data->output), "\n");
//	}
//}

static void printIgnSensitivities(realtype sensCoeffs[], int indices[], UserData data)
{
	char buf[50];
	std::string line;
	for(int j=0;j<data->nreac;j++){
		line=data->gas->reactionString(indices[j]);
		sprintf(buf,"%s",line.c_str());
		fprintf((data->ignSensOutput), "%d\t\"%35s\"\t%15.6e\n",indices[j],buf,sensCoeffs[j]);
	}
}

static int parsedPdt(UserData data){
  	/* Open file containing P and dPdt as functions of time: */
	data->dPdtInput=fopen("dPdt.dat","r");
	if(data->dPdtInput==NULL){
		printf("The dPdt.dat file wasn't found!\n");
		return(-1);
	}
	char buf[1000];
	char comment[1];
	char *ret;
	int i=0;
	int j=0;
	while (fgets(buf,100, data->dPdtInput)!=NULL){
		comment[0]=buf[0];
		if(strncmp(comment,"#",1)!=0 &&
		   strncmp(comment,"\n",1)!=0){
			i++;
		}
	}
	printf("No: of rows with numbers: %d\n",i);
	rewind(data->dPdtInput);
	int nPts=i;
	data->tArraySize=nPts;
	i=0;

  	data->tArray 	= N_VNew_Serial(nPts);
  	data->PArray 	= N_VNew_Serial(nPts);
  	data->dPdtArray = N_VNew_Serial(nPts);

	while (fgets(buf,100, data->dPdtInput)!=NULL){
		comment[0]=buf[0];
		if(strncmp(comment,"#",1)==0 ||
		   strncmp(comment,"\n",1)==0){
			printf("Comment! Skip Line!\n");
		}
		else{
			j=0;
			ret=strtok(buf,", \t");
			t(i)=atof(ret);
			j++;
			while(ret!=NULL){
   		  		ret=strtok(NULL,", \t");
				if(j==1){
					P(i)=atof(ret);
				}
				else if(j==2){
					dPdt(i)=atof(ret);
				}
				j++;
			}
			i++;
   		}
	}
	fclose(data->dPdtInput);
	//for(int k=0;k<nPts;k++){
	//	printf("%15.6e\t%15.6e\t%15.6e\n",t(k),P(k),dPdt(k));
	//}
	data->t1=t(nPts-1);
	data->P=P(0)*Cantera::OneAtm;

	/*check the polynomial interpolation (testing only)*/
	//realtype t,P,dPdt;
	//t=1e-03;
	//P=dPdt=0.0e0;
	//lookupDpdt(data, t, &P, &dPdt);
	//printf("Interpolated P value %15.6e. \n",P);
	//printf("Interpolated dPdt value %15.6e. \n",dPdt);
	//
	
	//GSL additions here:
	data->acc=gsl_interp_accel_alloc();
	data->spline=gsl_spline_alloc(gsl_interp_steffen,data->tArraySize);

	data->accdot=gsl_interp_accel_alloc();
	data->splinedot=gsl_spline_alloc(gsl_interp_steffen,data->tArraySize);

	double* Pdata;
	double* dPdtdata;
	double* tdata;

	tdata=N_VGetArrayPointer_Serial(data->tArray);
	Pdata=N_VGetArrayPointer_Serial(data->PArray);
	dPdtdata=N_VGetArrayPointer_Serial(data->dPdtArray);

	gsl_spline_init(data->spline,tdata,Pdata,data->tArraySize);
	gsl_spline_init(data->splinedot,tdata,dPdtdata,data->tArraySize);

	return(0);
}

static void PrintFinalStats(void *cvode_mem)
{
	long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
	int flag;
	
	flag = CVodeGetNumSteps(cvode_mem, &nst);
	check_flag(&flag, "CVodeGetNumSteps", 1);
	flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
	check_flag(&flag, "CVodeGetNumRhsEvals", 1);
	flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
	check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
	flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
	check_flag(&flag, "CVodeGetNumErrTestFails", 1);
	flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
	check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
	flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
	check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
	
	flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
	check_flag(&flag, "CVDlsGetNumJacEvals", 1);
	flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
	check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
	
	flag = CVodeGetNumGEvals(cvode_mem, &nge);
	check_flag(&flag, "CVodeGetNumGEvals", 1);
	
	printf("\nFinal CVode Statistics:\n");
	printf("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
	       nst, nfe, nsetups, nfeLS, nje);
	printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n \n",
		 nni, ncfn, netf, nge);
}

static void PrintFinalStatsKS(void *kmem)
{
	long int nni, nfe, nje, nfeD;
	long int lenrw, leniw, lenrwB, leniwB;
	long int nbcfails, nbacktr;
	int flag;
	
	/* Main solver statistics */
	
	flag = KINGetNumNonlinSolvIters(kmem, &nni);
	check_flag(&flag, "KINGetNumNonlinSolvIters", 1);
	flag = KINGetNumFuncEvals(kmem, &nfe);
	check_flag(&flag, "KINGetNumFuncEvals", 1);
	
	/* Linesearch statistics */
	
	flag = KINGetNumBetaCondFails(kmem, &nbcfails);
	check_flag(&flag, "KINGetNumBetacondFails", 1);
	flag = KINGetNumBacktrackOps(kmem, &nbacktr);
	check_flag(&flag, "KINGetNumBacktrackOps", 1);
	
	/* Main solver workspace size */
	
	flag = KINGetWorkSpace(kmem, &lenrw, &leniw);
	check_flag(&flag, "KINGetWorkSpace", 1);
	
	/* Band linear solver statistics */
	
	flag = KINDlsGetNumJacEvals(kmem, &nje);
	check_flag(&flag, "KINDlsGetNumJacEvals", 1);
	flag = KINDlsGetNumFuncEvals(kmem, &nfeD);
	check_flag(&flag, "KINDlsGetNumFuncEvals", 1);
	
	/* Band linear solver workspace size */
	
	flag = KINDlsGetWorkSpace(kmem, &lenrwB, &leniwB);
	check_flag(&flag, "KINDlsGetWorkSpace", 1);
	
	printf("\nFinal KINSOL Statistics:\n");
	printf("nni      = %6ld    nfe     = %6ld \n", nni, nfe);
	printf("nbcfails = %6ld    nbacktr = %6ld \n", nbcfails, nbacktr);
	printf("nje      = %6ld    nfeB    = %6ld \n", nje, nfeD);
	printf("\n");
	printf("lenrw    = %6ld    leniw   = %6ld \n", lenrw, leniw);
	printf("lenrwB   = %6ld    leniwB  = %6ld \n", lenrwB, leniwB);

}

static int check_flag(void *flagvalue, const char *funcname, int opt)
{
	int *errflag;
	
	/* Check if SUNDIALS function returned NULL pointer - no memory allocated */
	if (opt == 0 && flagvalue == NULL) {
	  fprintf(stderr, 
	          "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
	          funcname);
	  return(1);
	}
	
	/* Check if flag < 0 */
	else if (opt == 1) {
	  errflag = (int *) flagvalue;
	  if (*errflag < 0) {
	    fprintf(stderr,
	            "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
	            funcname, *errflag);
	    return(1); 
	  }
	}
	
	/* Check if function returned NULL pointer - no memory allocated */
	else if (opt == 2 && flagvalue == NULL) {
	  fprintf(stderr,
	          "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
	          funcname);
	  return(1);
	}
	
	return(0);
}