#ifndef SUNDIALS_DEF
#define SUNDIALS_DEF
#include <sundials/sundials_types.h>
#include <nvector/nvector_openmp.h>
#endif

#ifndef PRINT_DEF
#define PRINT_DEF
#include <string.h>	//for strings
#include <stdio.h>	//for printf,scanf
#include <stdlib.h>	//for atoi, atof
#endif

#ifndef CANTERA_DEF
#define CANTERA_DEF
#include <cantera/IdealGasMix.h>
#include <cantera/transport.h>
#endif

#include "UserData.h"


double maxGradPosition(const double* y, const size_t nt, 
	               const size_t nvar, const double* x, size_t nPts);

int maxGradIndex(const double* y, const size_t nt, 
	         const size_t nvar, const double* x, size_t nPts);

double maxCurvPosition(const double* y, const size_t nt, 
	               const size_t nvar, const double* x, size_t nPts);

int maxCurvIndex(const double* y, const size_t nt, 
	               const size_t nvar, const double* x, size_t nPts);

double isothermPosition(const double* y, const double T, const size_t nt, 
	                const size_t nvar, const double* x, const size_t nPts);

int setAlgebraicVariables(N_Vector *id,UserData data);

inline double calc_area(double x,int* i);

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

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

double systemMass(double* ydata, UserData data);

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

int setInitialCondition(N_Vector* y, 
	    		N_Vector* ydot, 
            		UserData data);


inline void setGas(UserData data, double *ydata, size_t gridPoint);

void getTransport(UserData data, 
		  double *ydata, 
		  size_t gridPoint,
		  double *rho,
		  double *lambda,
		  double *YV);

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

void trackFlameOH(N_Vector y,UserData data);
void trackFlame(N_Vector y,UserData data);
size_t BathGasIndex(UserData data);
size_t oxidizerIndex(UserData data);

inline double Qdot(double* t, 
	           double* x, 
	           double* ignTime, 
	           double* kernelSize, 
	           double* maxQdot);

void printSpaceTimeHeader(UserData data);
void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data);
void printSpaceTimeRates(double t, N_Vector ydot, UserData data);
void printGlobalHeader(UserData data);
void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data);
void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data);

void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data);
void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data);