#ifndef GSL_DEF
#define GSL_DEF

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

#endif

#include "residue.h"
#include "macros.h"
#include <cmath>
#include <stdio.h>
#include "timing.hpp"

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

    for (size_t i = 2; i <= nPts; i++) {
        TV1 = TV1 + fabs(T(i) - T(i - 1));
    }

    for (size_t i = 2; i <= nPts - 1; i++) {
        TV2 = TV2 + fabs((T(i + 1) - T(i)) / (R(i + 1) - R(i)) - (T(i) - T(i - 1)) / (R(i) - R(i - 1)));
    }

    XLEN = R(nPts) - R(1);
    B1 = R1 * XLEN / ((1.0 - P) * TV1);
    B2 = R2 * XLEN / ((1.0 - P) * TV2);

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

    DDX = R(nPts) - R(nPts - 1);
    WORK(nPts) = WORK(nPts - 1) + DDX;

    for (size_t i = 2; i <= nPts; i++) {
        WORK(i) = WORK(i) / WORK(nPts);
    }

    XNEW(1) = R(1);
    XNEW(nPts) = R(nPts);
    ISTART = 2;
    //double DETA;
    DETA = 1.0 / (nPts - 1);

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

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

        }
    }

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

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


void INTERPO(double *y, double *ydot, const size_t nvar, size_t nPts, const double XNEW[], const double XOLD[]) {
    double ytemp[nPts], ydottemp[nPts];

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

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

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

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

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

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

    }

    gsl_interp_accel_free(acc);
    gsl_spline_free(spline);

    gsl_interp_accel_free(accdot);
    gsl_spline_free(splinedot);

}

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

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

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

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

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

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

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

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

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

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

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

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

    double ytemp[nPts], ydottemp[nPts];

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

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

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

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

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

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

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

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

    gsl_interp_accel_free(acc);
    gsl_spline_free(spline);

    gsl_interp_accel_free(accdot);
    gsl_spline_free(splinedot);
}

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

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

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

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

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

/*a similar function will be used*/
/*therefore the original function will be commented out*/
// int setAlgebraicVariables(N_Vector* id, UserData data,const double* ydata){
// 	double *iddata;
//     //char* specPtr[2];

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

// 	for (size_t i = 0; i < data->dropMole.size(); i++){
// 		size_t index = specIndex(data,data->dropSpec[i]);
// 		data->k_drop.push_back(index);
// 	}

// 	printf("Oxidizer index: %lu\n",data->k_oxidizer);
// 	printf("Bath gas index:%lu\n",data->k_bath);
// 	printf("Droplet species index: ");
// 	for (size_t i = 0; i < data->dropMole.size(); i++)
// 	{
// 		printf("%lu\t",data->k_drop[i]);
// 	}
// 	printf("\n");

// 	/*set governing equations' id*/
// 	for (size_t i = 1; i <=data->npts; i++) {
// 		/*Algebraic variables: indicated by ZERO.*/
// 		Rid(i)=ZERO;
// 		Pid(i)=ZERO;
// 		Yid(i,data->k_bath)=ZERO;
// 		Mdotid(i)=ZERO;
// 	}
// 	Mdotid(1)=ZERO;
// 	Rid(1)=ONE;
// 	Yid(1,data->k_drop[0])=ZERO;
//     Yid(1,data->k_drop[1])=ZERO;
// 	Tid(data->npts)=ZERO;
// 	if(data->constantPressure){
// 		Pid(data->npts)=ONE;
// 	}else{
// 		Pid(data->npts)=ZERO;
// 		Rid(data->npts)=ONE;
// 	}
// 	if(data->dirichletInner){
// 		Tid(1)=ZERO;
// 	}else{
// 		//double T_boil = getLiquidMaxT(data->dropMole,P(1));
// 		//if(T(1) <= T_boil){
// 		//	Tid(1)=ONE;
// 		//}else{
// 		//	Tid(1)=ZERO;
// 		//}
// 		Tid(1)=ONE;
// 	}
// 	if(data->dirichletOuter){
// 		for (size_t k = 1; k <=data->nsp; k++) {
// 			if(k!=data->k_bath){
// 				Yid(data->npts,k)=ONE;
// 			}
// 			Yid(1,k)=ZERO;
// 			Tid(data->npts)=ONE;
// 		}
// 	}else{
// 		for (size_t k = 1; k <=data->nsp; k++){
// 			Yid(1,k)=ZERO;
// 			Yid(data->npts,k)=ZERO;
// 			Tid(data->npts)=ZERO;
// 		}
// 	}
// 	return(0);
// }

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

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

    for (size_t i = 0; i < data->dropMole.size(); i++) {
        size_t index = 0;
        for (size_t k = 1; k <= data->nsp; k++) {
            if (data->gas->speciesName(k - 1) == data->dropSpec[i]) {
                index = k;
            }
        }
//        size_t index = specIndex(data, data->dropSpec[i].c_str());
        data->k_drop.push_back(index);
    }

    printf("Oxidizer index: %lu\n", data->k_oxidizer);
    printf("Bath gas index:%lu\n", data->k_bath);
    printf("Droplet species index: ");

    for(auto & arg : data->k_drop){
        printf("%lu,\t",arg) ;
    }
    printf("\n");

    /*set governing equations' id*/
    for (size_t i = 1; i <= data->npts; i++) {
        /*Algebraic variables: indicated by ZERO.*/
        Rid(i) = ZERO;
        Pid(i) = ZERO;
        Yid(i, data->k_bath) = ZERO;
        Mdotid(i) = ZERO;
        Tid(i) = ONE;
    }
    Mdotid(1) = ZERO;

    /*set radius ids*/
    Rid(1) = ONE;
    Rid(data->l_npts) = ONE;

    /*set temperature ids*/
    Tid(1) = ZERO;
    Tid(data->l_npts) = ZERO;
    Tid(data->l_npts + 1) = ZERO;

    /*set liquid phase mass fraction (y) ids*/
    if (data->dropType == 0) {
        for (size_t k = 1; k <= data->nsp; k++) {
            for (size_t i = 1; i <= data->l_npts; i++) {
                if (i == 1) {
                    Yid(i, k) = ONE;
                } else {
                    Yid(i, k) = ZERO;
                }
            }
        }
    } else {
        for (size_t i = 1; i <= data->l_npts; i++) {
            for (size_t k = 1; k <= data->nsp; k++) {
                if (i == 1) {
                    if (k == data->k_drop[0] || k == data->k_drop[1]) {
                        Yid(i, k) = ZERO;
                    } else {
                        Yid(i, k) = ONE;
                    }
                } else if (i == data->l_npts) {
                    Yid(i, k) = ZERO;
                } else {
                    if (k == data->k_drop[0]) {
                        Yid(i, k) = ONE;
                    } else {
                        Yid(i, k) = ZERO;
                    }
                }
            }
        }

    }

//    for (size_t k = 1; k <= data->nsp; k++) {
//        Yid(data->l_npts + 1, k) = ZERO;
//    }

    /*set gas phase mass fraction (y) ids*/
    for (size_t i = (data->l_npts + 1); i <= (data->npts - 1); i++) {
        for (size_t k = 1; k <= data->nsp; k++) {
            if (i == (data->l_npts + 1)) {
                Yid(i, k) = ZERO;
            } else {
                if (k == data->k_bath) {
                    Yid(i, k) = ZERO;
                } else {
                    Yid(i, k) = ONE;
                }
            }
        }
    }

    //Yid(1,data->k_drop[0])=ZERO;
    //Yid(1,data->k_drop[1])=ZERO;
    //Tid(data->npts)=ZERO;

    // if(data->dirichletInner){
    // 	Tid(1)=ZERO;
    // }else{
    // 	//double T_boil = getLiquidMaxT(data->dropMole,P(1));
    // 	//if(T(1) <= T_boil){
    // 	//	Tid(1)=ONE;
    // 	//}else{
    // 	//	Tid(1)=ZERO;
    // 	//}
    // 	Tid(1)=ONE;
    // }
    if (data->constantPressure) {
        Pid(data->npts) = ONE;
    } else {
        Pid(data->npts) = ZERO;
        /*Rres(npts) seems to be algebric all the time?*/
        Rid(data->npts) = ONE;
    }
    if (data->dirichletOuter) {
        for (size_t k = 1; k <= data->nsp; k++) {
            if (k != data->k_bath) {
                Yid(data->npts, k) = ONE;
            }
            //Yid(1,k)=ZERO;
            Tid(data->npts) = ONE;
        }
    } else {
        for (size_t k = 1; k <= data->nsp; k++) {
            //Yid(1,k)=ZERO;
            Yid(data->npts, k) = ZERO;
            Tid(data->npts) = ZERO;
        }
    }

    printIddata(data,iddata);
    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);
    }
}
/*Similar to function:"setInitialCondition"*/
/*"readInitialCondition" will also be revised here for two-phase LTORC*/
// void readInitialCondition(FILE* input, double* ydata, const size_t nvar, const size_t nr, const size_t nPts, double Rg){

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

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

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

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

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

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

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

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

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

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

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

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

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

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

// 		gsl_spline_init(spline,xOld,ytemp,nRows);

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

// 	gsl_interp_accel_free(acc);
// 	gsl_spline_free(spline);

// }

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

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

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

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

    //get the # of rows in the input file:"initialCondition.dat, which equals to # of grid points"
    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];

    /*extract the initial data from file to y[] array*/
    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++;
        }
    }

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

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

    /*get the x(radial) coordinate from y[] array*/
    for (size_t j = 0; j < nRows; j++) {
        xOld[j] = y[j * nColumns + nr];
    }

    /*get the first index of droplet interface:adjacent r equals*/
    size_t interII;
    for (size_t i = 0; i < nRows; i++) {
        interII = i;
        if (xOld[i] == xOld[i + 1]) {
            break;
        }
    }

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

    double dX0, dXL;
    int NN = nPts - l_nPts - 2;

    /*re-calculate the spatial coordinate: uniform grid in liquid phase and non-uniform for gas(fixed ratio)*/
    dXL = (xOld[interII] - xOld[0]) / (l_nPts - 1);
    // dX0 = (xOld[nRows-1]-xOld[0])*(pow(Rg,1.0/NN) - 1.0)/(pow(Rg,(NN+1.0)/NN) - 1.0);
    xNew[0] = xOld[0];
    for (size_t i = 1; i < l_nPts; i++) {
        xNew[i] = xNew[i - 1] + dXL;
    }
    xNew[l_nPts - 1] = xOld[interII];
    xNew[l_nPts] = xNew[l_nPts-1];

    dX0 = (xOld[nRows - 1]- xOld[interII] ) * (pow(Rg, 1.0 / NN) - 1.0) / (pow(Rg, (NN + 1.0) / NN) - 1.0);

    double ratio = pow(Rg,1.0/NN);
    for (size_t ii = (l_nPts + 1); ii <= (nPts - 2); ii++) {
        xNew[ii] = xNew[ii - 1] + dX0 * pow(ratio,(ii-l_nPts-1)) ;
    }
    xNew[nPts - 1] = xOld[nRows - 1];

    /*pass variable data from ytemp[] to ydata[] array using spline interpolation*/
    /*interpolate the liquid phase state variables with old array of size (interII+1)*/
    double ytemp0[interII + 1];
    double xOld0[interII + 1];
    std::copy(xOld, xOld + interII + 1, xOld0);
    gsl_interp_accel *acc;
    gsl_spline *spline;
    acc = gsl_interp_accel_alloc();
    spline = gsl_spline_alloc(gsl_interp_steffen, interII + 1);
    for (size_t j = 0; j < nColumns; j++) {
        for (size_t k = 0; k < (interII + 1); k++) {
            ytemp0[k] = y[j + k * nColumns];
        }
        gsl_spline_init(spline, xOld0, ytemp0, interII + 1);

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

    /*interpolate the gas phase state variables with old array of size (nRows-interII-1)*/
    double ytemp1[nRows - interII - 1];
    double xOld1[nRows - interII - 1];
    std::copy(xOld + interII + 1, xOld + nRows, xOld1);
    gsl_interp_accel *acc1;
    gsl_spline *spline1;
    acc1 = gsl_interp_accel_alloc();
    spline1 = gsl_spline_alloc(gsl_interp_steffen, nRows - interII - 1);
    for (size_t j = 0; j < nColumns; j++) {

        for (size_t k = (interII + 1); k < nRows; k++) {
            // ytemp[k]=y[j+k*nColumns];
            ytemp1[k - interII - 1] = y[j + k * nColumns];
        }

        gsl_spline_init(spline1, xOld1, ytemp1, nRows - interII - 1);

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

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

/*a different function:initializePsiEtaGrid is implement*/
/*similar to the following function, it create a psi grid in the gas phase*/
/*and a eta grid in the liquid phase*/
int initializePsiGrid(double *ydata, double *psidata, UserData data) {

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

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

/*create the psi grid in the gas phase and eta grid in the liquid phase*/
/*note: psi has units of kg. while eta is dimensionless*/
int initializePsiEtaGrid(double *ydata, double *psidata, UserData data) {
    double rho, rhom, rhop;
    /*Create a psi grid that corresponds CONSISTENTLY to the spatial grid
	 * "R" created above. Note that the Lagrangian variable psi has  units
	 * of kg. */
    /*create the psi grid in gas phase*/
    psi(data->l_npts + 1) = ZERO;
    for (size_t i = data->l_npts + 2; i <= data->npts; i++) {
        data->gas->setState_TPY(T(i), P(i), &Y(i, 1));
        rho = data->gas->density();
        data->gas->setState_TPY(T(i - 1), P(i - 1), &Y(i - 1, 1));
        rhom = data->gas->density();
        //psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&data->metric);
        //psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
        psi(i) = psi(i - 1) + (R(i) - R(i - 1)) *
                              (rho * calc_area(R(i), &data->metric) + rhom * calc_area(R(i - 1), &data->metric)) / TWO;
    }

    /*create the \eta grid for liquid phase*/
    /*for the time being, liquid phase species are hard coded*/
    std::vector<std::string> composition = components(data->dropType);
    psi(data->l_npts) = ZERO;
    if (data->dropType == 0) {
        for (size_t i = data->l_npts - 1; i >= 1; i--) {
            rho = getLiquidDensity(T(i), P(i), composition);
            rhop = getLiquidDensity(T(i + 1), P(i + 1), composition);
            psi(i) = psi(i + 1) + (R(i + 1) - R(i)) *
                                  (rho * calc_area(R(i), &data->metric) + rhop * calc_area(R(i + 1), &data->metric)) /
                                  TWO;
        }
    } else if (data->dropType == 1) {
        for (size_t i = data->l_npts - 1; i >= 1; i--) {
            std::vector<double> moleFrac, moleFracp;
            moleFrac = getLiquidmolevec(data, ydata, i);
            moleFracp = getLiquidmolevec(data, ydata, i + 1);
            rho = getLiquidDensity(T(i), P(i), composition, moleFrac);
            rhop = getLiquidDensity(T(i + 1), P(i + 1), composition, moleFracp);
            psi(i) = psi(i + 1) + (R(i + 1) - R(i)) *
                                  (rho * calc_area(R(i), &data->metric) + rhop * calc_area(R(i + 1), &data->metric)) /
                                  TWO;
        }
    }
    /*according to the mehtod described by Stauch et al.*/
    /*eta should be dimensionless, so we normalize the eta with total mass of droplet*/
    double dropmass = psi(1);
    for (size_t i = data->l_npts; i >= 1; i--) {
        psi(i) = psi(i) / dropmass;
    }

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

/*A different setInitialCondition is implemented based on the problemType =2 */
/*As such, the following function is commented*/
// int setInitialCondition(N_Vector* y,
// 	    		N_Vector* ydot,
//             		UserData data){
// 	double* ydata;
// 	double* ydotdata;
// 	double* psidata;
// 	double* innerMassFractionsData, Rd, massDrop;
// 	double rhomhalf, lambdamhalf, YVmhalf[data->nsp];
// 	double f=ZERO;
// 	double g=ZERO;

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

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

// 	m=data->metric;

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

// 	//Define Grid:
// 	double dR=(data->domainLength)/((double)(data->npts)-1.0e0);
// 	double dv=(pow(data->domainLength,1+data->metric)-pow(data->firstRadius*data->domainLength,1+data->metric))/((double)(data->npts)-1.0e0);

//     //Initialize the R(i),T(i)(data->initialTemperature),Y(i,k),P(i)
//     for (size_t i = 1; i <=data->npts; i++) {
// 		if(data->metric==0){
// 			R(i)=Rd+(double)((i-1)*dR);
// 		}else{
// 			if(i==1){
// 				R(i)=ZERO;
// 			}else if(i==2){
// 				R(i)=data->firstRadius*data->domainLength;
// 			}else{
// 				R(i)=pow(pow(R(i-1),1+data->metric)+dv,1.0/((double)(1+data->metric)));
// 			}
// 		}
// 		T(i)=data->initialTemperature;
// 		for (size_t k = 1; k <=data->nsp; k++) {
//   	        	Y(i,k)=data->gas->massFraction(k-1);	//Indexing different in Cantera
// 		}
//   	        P(i)=data->initialPressure*Cantera::OneAtm;
// 	}
// 	R(data->npts)=data->domainLength+data->Rd;

// //    /********** test R(i) and the volumn between grids *****************/
// //    printf("Print the first 4 R(i)s and volumn between them: \n") ;
// //    printf("Grids: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e \n",R(1),R(2),R(3),R(4));
// //    double v1,v2,v3,v4,v5;
// //    v1 =pow(R(2),1+data->metric) - pow(R(1),1+data->metric);
// //    v2 =pow(R(3),1+data->metric) - pow(R(2),1+data->metric);
// //    v3 =pow(R(4),1+data->metric) - pow(R(3),1+data->metric);
// //    v4 =pow(R(5),1+data->metric) - pow(R(4),1+data->metric);
// //    v5 =pow(R(6),1+data->metric) - pow(R(5),1+data->metric);
// //    printf("Volumn: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e,5th:%2.6e\n",v1,v2,v3,v4,v5) ;


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

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

// 	initializePsiGrid(ydata,psidata,data);

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

//             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 Rg = data->Rg, dpsi0;
// 		int NN = data->npts-2;
// 		double psiNew[data->npts];

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

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

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

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

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

// 	//Set grid to location of maximum curvature calculated by r instead of x:
// 	if(data->adaptiveGrid){
// 		//data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
// 		//		  		data->grid->x, data->npts);
//         int maxCurvII = 0;
//         maxCurvII = maxCurvIndexR(ydata,data->nt,data->nvar,data->nr,data->npts);

//         data->grid->position = data->grid->x[maxCurvII];
// 		ier=reGrid(data->grid, data->grid->position);
// 		updateSolution(ydata, ydotdata, data->nvar,
// 		               data->grid->xOld,data->grid->x,data->npts);
// 		storeGrid(data->grid->x,data->grid->xOld,data->npts);

//         /******** Test the maxCurvPosition and related variables *******/
//         printf("The maxCurvPositition(based on r) is : %.6f,Temperature is : %.3f \n",data->grid->position,T(maxCurvII+1));
// 	}


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

// 	return(0);
// }

/*Revise the setInitialCondition function to be suitable for two-phase TORC*/
int setInitialCondition(N_Vector *y,
                        N_Vector *ydot,
                        UserData data) {
    double *ydata;
    double *ydotdata;
    double *psidata;
    double *innerMassFractionsData, Rd, massDrop;
    double rhomhalf, lambdamhalf, YVmhalf[data->nsp];
    double f = ZERO;
    double g = ZERO;

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

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

    m = data->metric;

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

    //Define Grid:
    double dR = (data->domainLength) / ((double) (data->g_npts) - 1.0e0);
    double dv = (pow(data->domainLength, 1 + data->metric) -
                 pow(data->firstRadius * data->domainLength, 1 + data->metric)) / ((double) (data->g_npts) - 1.0e0);


    // //Initialize the R(i),T(i)(data->initialTemperature),Y(i,k),P(i)
    // for (size_t i = 1; i <=data->npts; i++) {
    // 	if(data->metric==0){
    // 		R(i)=Rd+(double)((i-1)*dR);
    // 	}else{
    // 		if(i==1){
    // 			R(i)=ZERO;
    // 		}else if(i==2){
    // 			R(i)=data->firstRadius*data->domainLength;
    // 		}else{
    // 			R(i)=pow(pow(R(i-1),1+data->metric)+dv,1.0/((double)(1+data->metric)));
    // 		}
    // 	}
    // 	T(i)=data->initialTemperature;
    // 	for (size_t k = 1; k <=data->nsp; k++) {
    //         	Y(i,k)=data->gas->massFraction(k-1);	//Indexing different in Cantera
    // 	}
    //         P(i)=data->initialPressure*Cantera::OneAtm;
    // }
    // R(data->npts)=data->domainLength+data->Rd;


//    /********** test R(i) and the volumn between grids *****************/
//    printf("Print the first 4 R(i)s and volumn between them: \n") ;
//    printf("Grids: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e \n",R(1),R(2),R(3),R(4));
//    double v1,v2,v3,v4,v5;
//    v1 =pow(R(2),1+data->metric) - pow(R(1),1+data->metric);
//    v2 =pow(R(3),1+data->metric) - pow(R(2),1+data->metric);
//    v3 =pow(R(4),1+data->metric) - pow(R(3),1+data->metric);
//    v4 =pow(R(5),1+data->metric) - pow(R(4),1+data->metric);
//    v5 =pow(R(6),1+data->metric) - pow(R(5),1+data->metric);
//    printf("Volumn: 1st:%2.6e,2nd:%2.6e,3rd:%2.6e,4th:%2.6e,5th:%2.6e\n",v1,v2,v3,v4,v5) ;

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

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

    // initializePsiGrid(ydata,psidata,data);


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

    //         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 Rg = data->Rg, dpsi0;
    // 	int NN = data->npts-2;
    // 	double psiNew[data->npts];

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

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

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

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

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

    //Set grid to location of maximum curvature calculated by r instead of x:
    if (data->adaptiveGrid) {
        //data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
        //		  		data->grid->x, data->npts);
        int maxCurvII = 0;
        maxCurvII = maxCurvIndexR(ydata, data->nt, data->nvar, data->nr, data->npts);

        data->grid->position = data->grid->x[maxCurvII];
        ier = reGrid(data->grid, data->grid->position);
        updateSolution(ydata, ydotdata, data->nvar,
                       data->grid->xOld, data->grid->x, data->npts);
        storeGrid(data->grid->x, data->grid->xOld, data->npts);

        /******** Test the maxCurvPosition and related variables *******/
        printf("The maxCurvPositition(based on r) is : %.6f,Temperature is : %.3f \n", data->grid->position,
               T(maxCurvII + 1));
    }


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

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

void getGasMassFlux(UserData data,
                  double *ydata,
                  size_t gridPoint,
                  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));

    data->trmix->getSpeciesFluxes(1, &gradT, data->nsp,
                                  gradX, data->nsp, YV);
    //setGas(data,ydata,gridPoint);
}

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

    /*Declare and fetch nvectors and user data:*/
    double *ydata, *ydotdata, *resdata, *psidata, *innerMassFractionsData;

    UserData data;
    data = (UserData) user_data;
    size_t npts = data->npts;
    size_t nsp = data->nsp;
    size_t k_bath = data->k_bath;
    size_t l_npts = data->l_npts;
    int dropType = data->dropType;

    std::vector<size_t> k_drop = data->k_drop;
    std::vector<std::string> dropSpec = data->dropSpec;
    std::vector<std::string> composition = components(data->dropType);
    std::vector<double> dropMole = data->dropMole;

    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, Diffcoeffphalf, Diffcoeffmhalf,propYVmhalf,propYVphalf,
            rhophalf, lambdaphalf, YVphalf[nsp],
            Cpb, Cvb, Cp[nsp], Cpl[2], dHvl[2], rhol, wdot[nsp], enthalpy[nsp], energy[nsp], T_boil,
            tranTerm, diffTerm, srcTerm, advTerm,
            area, areamhalf, areaphalf, aream, areamhalfsq, areaphalfsq;

    /*Aliases for difference coefficients:*/
    double cendfm, cendfc, cendfp;
    /*Aliases for various grid spacings:*/
    double dpsip, dpsiav, dpsipm, dpsim, dpsimm;
    dpsip = dpsiav = dpsipm = dpsim = dpsimm = ONE;
    //double mass, mdotIn;
    double mass, massDrop;
    double sum, sum1, sum2, sum3;

    std::vector<double> vapheat, mole_,liquidCp,vapPres;

//    size_t j, k;
    int m;
    m = data->metric;                    //Unitless
    mass = data->mass;                //Units: kg

    /*evaluate and implement the governing equations in the liquid phase*/

    /*implement res @ left boundary i.e., first grid point*/
    Rres(1) = Rdot(1);
    if (dropType == 0) {
        // Yres(1,data->k_drop[0]) = Ydot(1,data->k_drop[0]);
        for (size_t k = 1; k <= nsp; k++) {
            Yres(1, k) = Ydot(1, k);
        }
    } else {
        for (size_t k = 1; k <= nsp; k++) {
            if (k == k_drop[0]) {
                Yres(1, k) = Y(2, k) - Y(1, k);
            } else if (k == k_drop[1]) {
                Yres(1, k) = ONE - Y(1, k_drop[0]) - Y(1, k);
            } else {
                Yres(1, k) = Ydot(1, k);
            }
        }
    }
    Pres(1) = P(2) - P(1);
    Mdotres(1) = Mdot(2) - Mdot(1);
    Tres(1) = T(2) - T(1);

    /*set governing equations at intermediate grid points */
    massDrop = dropletmass(data, ydata);
    if (dropType == 0) {
        for (size_t j = 2; j <= l_npts - 1; j++) {
            /*Mass:*/
            /* ∂r/∂ψ = 1/ρA */
            dpsim = - (psi(j) - psi(j - 1));
            dpsip = - (psi(j + 1) - psi(j));
            dpsiav = - (HALF * (psi(j + 1) - psi(j - 1)));
            rho = getLiquidDensity(T(j), P(j), composition);
            rhom = getLiquidDensity(T(j - 1), P(j - 1), composition);
            rhophalf = getLiquidDensity(HALF * (T(j) + T(j + 1)), HALF * (P(j) + P(j + 1)), composition);
            rhomhalf = getLiquidDensity(HALF * (T(j) + T(j - 1)), HALF * (P(j) + P(j - 1)), composition);

            Cpb = getLiquidCpb(T(j), P(j), composition);
            lambdaphalf = getLiquidCond(HALF * (T(j) + T(j + 1)), HALF * (P(j) + P(j + 1)), composition);
            lambdamhalf = getLiquidCond(HALF * (T(j) + T(j - 1)), HALF * (P(j) + P(j - 1)), composition);

            aream = calc_area(R(j - 1), &m);
            areamhalf = calc_area(HALF * (R(j) + R(j - 1)), &m);
            areaphalf = calc_area(HALF * (R(j + 1) + R(j)), &m);
            areamhalfsq = areamhalf * areamhalf;
            areaphalfsq = areaphalf * areaphalf;
            area = calc_area(R(j), &m);

            Rres(j) = ((R(j) - R(j - 1)) / dpsim) - massDrop * (TWO / (rhom * aream + rho * area));
            Mdotres(j) = Mdot(j + 1) - Mdot(j);
            Pres(j) = P(j + 1) - P(j);
            for (size_t k = 1; k <= nsp; k++) {
                Yres(j, k) = Y(j, k) - Y(j - 1, k);
            }

            tranTerm = Tdot(j);
            advTerm = psi(j) * (-Mdot(j)) / massDrop * (T(j) - T(j - 1)) / dpsim ;
            diffTerm = 1 / (Cpb * massDrop * massDrop) *
                       (rhophalf * areaphalfsq * lambdaphalf * (T(j + 1) - T(j)) / dpsip -
                        rhomhalf * areamhalfsq * lambdamhalf * (T(j) - T(j - 1)) / dpsim) / dpsiav;
            Tres(j) = tranTerm - advTerm - diffTerm;
        }

        /*Fill up the res at l_npts*/
        rho = getLiquidDensity(T(l_npts), P(l_npts), composition);
        area = calc_area(R(l_npts), &m);
        lambdamhalf = getLiquidCond(T(l_npts), P(l_npts), composition);
        lambdaphalf = getGasCond(data, ydata, l_npts+1);
        vapheat = getLiquidVH(P(l_npts+1), data->dropType);

        Rres(l_npts) = Mdot(l_npts) + rho * Rdot(l_npts) * area;
        Mdotres(l_npts) = Mdot(l_npts + 1) - Mdot(l_npts);
        Pres(l_npts) = P(l_npts + 1) - P(l_npts);
        for (size_t k = 1; k <= nsp; k++) {
            Yres(l_npts, k) = Y(l_npts, k) - Y(l_npts - 1, k);
        }
        Tres(l_npts) = lambdamhalf * (T(l_npts) - T(l_npts - 1)) / (R(l_npts) - R(l_npts - 1)) * area +
                       Mdot(l_npts) * vapheat[0] -
                       lambdaphalf * (T(l_npts + 2) - T(l_npts + 1)) / (R(l_npts + 2) - R(l_npts + 1)) * area;
    } else {
        /*binary component case*/
        /*implementation of liquid phase inner grid points*/
        for (size_t j = 2; j <= l_npts - 1; j++) {
            mole_ = getLiquidmolevec(data, ydata, j);

            dpsim = - (psi(j) - psi(j - 1));
            dpsip = - (psi(j + 1) - psi(j));
            dpsiav = - (HALF * (psi(j + 1) - psi(j - 1)));

            /*evaluate various central difference coefficients*/
            cendfm = -dpsip / (dpsim * dpsipm);
            cendfc = (dpsip - dpsim) / (dpsip * dpsim);
            cendfp = dpsim / (dpsip * dpsipm);
            /**************************************************/

            rho = getLiquidDensity(T(j), P(j), composition, mole_);
            rhom = getLiquidDensity(T(j - 1), P(j - 1), composition, mole_);
            rhophalf = getLiquidDensity(HALF * (T(j) + T(j + 1)), HALF * (P(j) + P(j + 1)), composition, mole_);
            rhomhalf = getLiquidDensity(HALF * (T(j) + T(j - 1)), HALF * (P(j) + P(j - 1)), composition, mole_);

            Cpb = getLiquidCpb(T(j), P(j), composition, mole_);
            liquidCp = getLiquidCp(T(j),P(j),composition);

            lambdaphalf = getLiquidCond(HALF * (T(j) + T(j + 1)), HALF * (P(j) + P(j + 1)), composition, mole_);
            lambdamhalf = getLiquidCond(HALF * (T(j) + T(j - 1)), HALF * (P(j) + P(j - 1)), composition, mole_);

            aream = calc_area(R(j - 1), &m);
            areamhalf = calc_area(HALF * (R(j) + R(j - 1)), &m);
            areaphalf = calc_area(HALF * (R(j + 1) + R(j)), &m);
            areamhalfsq = areamhalf * areamhalf;
            areaphalfsq = areaphalf * areaphalf;
            area = calc_area(R(j), &m);

            Diffcoeffphalf = getLiquidmassdiff(data, ydata, j, HALF * (T(j) + T(j + 1)));
            Diffcoeffmhalf = getLiquidmassdiff(data, ydata, j, HALF * (T(j) + T(j - 1)));

            Rres(j) = ((R(j) - R(j - 1)) / dpsim) - massDrop * (TWO / (rhom * aream + rho * area));
            Mdotres(j) = Mdot(j + 1) - Mdot(j);
            Pres(j) = P(j + 1) - P(j);

             propYVphalf = (-Diffcoeffphalf * (Y(j + 1, k_drop[0]) - Y(j, k_drop[0])) / (R(j + 1) - R(j)));
             propYVmhalf = (-Diffcoeffmhalf * (Y(j, k_drop[0]) - Y(j - 1, k_drop[0])) / (R(j) - R(j - 1)));

            liquidCp = getLiquidCp(T(j),P(j),composition);
            /*species equation*/
            for (size_t k = 1; k <= nsp; k++) {
                if (k != k_drop[0] && k != k_drop[1]) {
                    Yres(j, k) = Y(j, k) - Y(j - 1, k);
                } else if (k == k_drop[0]) {
                    tranTerm = Ydot(j, k);
                    advTerm = 1 / massDrop * (psi(j) * (-Mdot(j))) * (Y(j, k) - Y(j - 1, k)) / dpsim;
                    diffTerm = 1 / massDrop *
                               (areaphalf * rhophalf * propYVphalf
                                - areamhalf * rhomhalf * propYVmhalf )
                                / dpsiav;
                    Yres(j, k) = tranTerm - advTerm + diffTerm;
                } else{
                    Yres(j,k) = ONE - Y(j,k_drop[0]);
                }
            }
            /*energy equation*/
            tranTerm  = Tdot(j);
            advTerm = 1/massDrop *(psi(j)*(-Mdot(j)))* (T(j)-T(j-1))/(psi(j)-psi(j-1));
            double diffTerm1 = 1/(Cpb*massDrop*massDrop) * ( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/(psi(j+1)-psi(j)) )
                    - (rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/(psi(j)-psi(j-1)) ) ) /dpsiav;
            double diffTerm2 = area/(Cpb*massDrop)*rho*(cendfp*T(j+1)+cendfc*T(j)+cendfm*T(j-1)) *
                    (HALF*((propYVphalf+propYVmhalf)*liquidCp[0]+(-propYVphalf-propYVmhalf)*liquidCp[1])) ;
            diffTerm = diffTerm1 - diffTerm2 ;
            Tres(j) = tranTerm - advTerm - diffTerm ;
        }

        /*implementation of right boundary of liquid phase*/
        mole_ = getLiquidmolevec(data,ydata,l_npts);
        rho = getLiquidDensity(T(l_npts), P(l_npts), composition,mole_);
        area = calc_area(R(l_npts), &m);
        lambdamhalf = getLiquidCond(T(l_npts), P(l_npts), composition,mole_);
        lambdaphalf = getGasCond(data, ydata, l_npts+1);

        propYVmhalf = (-Diffcoeffmhalf * (Y(l_npts, k_drop[0]) - Y(l_npts - 1, k_drop[0])) / (R(l_npts) - R(l_npts - 1)));
        getGasMassFlux(data,ydata,l_npts+1,YVphalf);

        Diffcoeffmhalf = getLiquidmassdiff(data, ydata, l_npts, HALF * (T(l_npts) + T(l_npts - 1)));
        vapheat = getLiquidVH(P(l_npts), data->dropType);

        Rres(l_npts) = Mdot(l_npts) + rho * Rdot(l_npts) * area;
        Mdotres(l_npts) = Mdot(l_npts + 1) - Mdot(l_npts);
        Pres(l_npts) = P(l_npts + 1) - P(l_npts);

        /*species formualtion*/
        for (size_t k = 1; k <= nsp; k++) {
            if (k == k_drop[0]) {
                Yres(l_npts, k) = Mdot(l_npts) * Y(l_npts, k) + area * rho *
                                                                (-Diffcoeffmhalf * (Y(l_npts, k) - Y(l_npts - 1, k)) /
                                                                 (R(l_npts) - R(l_npts - 1)))
                                  - Mdot(l_npts) * Y(l_npts + 1, k) - area * YVphalf[k_drop[0] - 1];
            } else if (k == k_drop[1]) {
                Yres(l_npts, k) = ONE - Y(l_npts, k_drop[0]) - Y(l_npts, k);
            } else {
                Yres(l_npts, k) = Y(l_npts, k) - Y(l_npts - 1, k);
            }
        }

        /*temperature formulation*/
        double epsilon = Y(l_npts + 1, k_drop[0]) + YVphalf[k_drop[0] - 1] * area / Mdot(l_npts);
        Tres(l_npts) = lambdamhalf * (T(l_npts) - T(l_npts - 1)) / (R(l_npts) - R(l_npts - 1)) * area
                       + Mdot(l_npts) * (epsilon * vapheat[0] + (1 - epsilon) * vapheat[1])
                       - lambdaphalf * (T(l_npts + 2) - T(l_npts + 1)) / (R(l_npts + 2) - R(l_npts + 1)) * area;
    }

    //massDrop=data->massDrop;		//Units: kg
    //mdotIn=data->mdot*calc_area(R(npts),&m);	//Units: kg/s

//	/*evaluate properties at j=l_npts+1 *************************/
    setGas(data, ydata, l_npts+1);
    rhom = data->gas->density();
    Cpb = data->gas->cp_mass();        //J/kg/K
    Cvb = data->gas->cv_mass();        //J/kg/K
    //TODO: Create user input model for these. They should not be constant
    //Cpl=4.182e03;						//J/kg/K for liquid water
    //Cpl=  5.67508633e-07*pow(T(1),4) - 7.78060597e-04*pow(T(1),3) + 4.08310544e-01*pow(T(1),2)  //J/kg/K for liquid water
    //     - 9.62429538e+01*T(1)        + 1.27131046e+04;

    /*Update specific heat:Cpl for liquid proprane&n-heptane&dodecane*/
    /*Based on the existiong data,a linear expression is used*/
    // Cpl= 12.10476*T(1) - 746.60143; // Unit:J/(kg*K), Need to be improved later
    //Cpl[0] = 2.423*T(1)+1661.074; //Unit:J/(kg*K),range:1atm,propane
    //Cpl[0] = 3.336*T(1)+1289.5; //Unit:J/(kg*K),range:1atm,n-Heptane
    //Cpl[1] = 3.815*T(1)+1081.8; //Unit:J/(kg*K),range:1atm,n-Dodecane
//    double Cp_l = 0.0;
//    Cp_l = getLiquidCp(data->dropMole,T(1),P(1));
    //for(size_t i=0;i<=1;i++){
    //    Cp_l=Cp_l+Cpl[i]*data->dropMassFrac[i];
    //}

    //dHvl=2.260e6;						//J/kg	heat of vaporization of water
    //double Tr = T(1)/647.1;				//Reduced Temperature: Droplet Temperature divided by critical temperature of water
    //dHvl= 5.2053e07*pow(1-Tr,0.3199 - 0.212*Tr + 0.25795*Tr*Tr)/18.01;	//J/kg	latent heat of vaporization of water

    //double Tr1 = T(1)/369.8;            //Reduced Temperature;
    //dHvl[0] = 6.32716e5*exp(-0.0208*Tr1)*pow((1-Tr1),0.3766);    //Unit:J/kg,Latent Heat of Vaporizaiton of Liquid n-propane,NIST

    //double Tr2 = T(1)/540.2;            //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST
    //dHvl[0] = 5.366e5*exp(-0.2831*Tr2) * pow((1-Tr2),0.2831);    //Unit:J/kg, latent heat of vaporization of liquid n-heptane, Source: NIST

    //dHvl[1] = 358.118e3;  //Unit:J/kg,latent heat of vaporization of n-dodecane,Source:NIST

//    double dHv_l = 0.0;
//    dHv_l = getLiquidHv(data->dropMole,T(1),P(1));
    //for(size_t i=0;i<=1;i++){
    //    dHv_l=dHv_l+dHvl[i]*data->dropMassFrac[i];
    //}

    /*Following section is related to the property of water*/
    //rhol = 997.0;
    //massDrop=1.0/3.0*R(1)*R(1)*R(1)*997.0; //TODO: The density of the droplet should be a user input

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

//    aream = calc_area(R(1), &m);


    /*******************************************************************/
    /*Calculate values at j=l_npts+2 's m and mhalf*****************************/
//    getTransport(data, ydata, 1, &rhomhalf, &lambdamhalf, YVmhalf);
//    areamhalf = calc_area(HALF * (R(1) + R(2)), &m);
//    aream = calc_area(R(1), &m);
//    areamhalfsq = areamhalf * areamhalf;

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

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

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


    /*Mass:*/
    if (data->quasiSteady) {
//        Rres(1) = Rdot(1); // Should remain satisfied for quasi-steady droplet
        Rres(l_npts+1) = Rdot(l_npts) ;
    } else {
        Rres(l_npts+1) = R(l_npts+1) - R(l_npts) ;
//        Rres(1) = Rdot(1) + Mdot(1) / rhol / aream;
    }

    //TODO: Antoine Parameters should be part of user input, not hardcoded
    //Yres(1,k_drop)=Y(1,k_drop) - 1.0e5*pow(10,4.6543-(1435.264/(T(1)-64.848)))*data->gas->molecularWeight(k_drop-1)
    //			    / P(1) / data->gas->meanMolecularWeight();
    //Yres(1,k_drop)=Y(1,k_drop) - 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)))*data->gas->molecularWeight(k_drop-1)
    //			    / P(1) / data->gas->meanMolecularWeight();
    //Yres(1,k_drop)=Y(1,k_drop)-2.566785e-02;

    /*Following using the Antoine Parameters to calculate the pressure of volatile component(n-heptane)*/
    /*Antoine parameter from NIST website*/
    /*Raoult's Law is applied to calculate the partial vapor pressure*/
//    double p_i[2] = {0.0, 0.0};
//    p_i[0] = 1.0e5 * pow(10, 4.53678 - (1149.36 / (T(1) + 24.906))) * data->dropMole[0] *
//             data->gamma[0];    //unit:Pa,Helgeson and Sage,1967,Propane
//    //p_i[0] = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199)))*data->dropMole[0] ;
//    p_i[1] = 1.0e5 * pow(10, 4.02832 - (1268.636 / (T(1) - 56.199))) * data->dropMole[1] *
//             data->gamma[1];    //unit:Pa,n-Heptane
    //p_i[1] = 1.0e5*pow(10,4.10549-(1625.928/(T(1)-92.839)))*data->dropMole[1] ;    //Unit:Pa.Williamham and Taylor,n-dodecane

    /*for both dropletTypes(0,1),implementations of Rres,Tres,Pres are the same
     * only Yres and Mdotres are different */
    /*Species:*/
    sum = ZERO;
    /*single component case*/
    if (dropType == 0) {
        double hepPres = 1.0e5 * pow(10, 4.02832 - (1268.636 / (T(l_npts+1) - 56.199)));
        Yres(l_npts + 1, k_drop[0]) = Y(l_npts + 1, k_drop[0]) - hepPres * data->gas->molecularWeight(k_drop[0] - 1) /
                                                                 (P(l_npts + 1) * data->gas->meanMolecularWeight());
        sum = sum + Y(l_npts + 1, k_drop[0]);
        Mdotres(l_npts + 1) = Mdot(l_npts + 1) - (YVmhalf(k_drop[0]) * areamhalf) / (1 - Y(l_npts + 1, k_drop[0]));
        for (size_t k = 1; k <= nsp; k++) {
            if (k != k_bath && k != k_drop[0]) {
                Yres(l_npts + 1, k) = Y(l_npts + 1, k) * Mdot(l_npts + 1) + YVmhalf(k) * areamhalf;
                sum = sum + Y(l_npts + 1, k);
            }
        }
        Yres(l_npts + 1, k_bath) = ONE - sum - Y(l_npts + 1, k_bath);
    }

    /*binary componet case for gas phase LHS boundary*/
    if(dropType == 1){
        mole_ = getLiquidmolevec(data,ydata,l_npts);
        vapPres = getVapPressure(data,ydata,l_npts+1,mole_);

        for (size_t i = 0; i < k_drop.size(); i++) {
            Yres(l_npts+1, k_drop[i]) = Y(l_npts+1, k_drop[i]) - vapPres[i] * data->gas->molecularWeight(k_drop[i] - 1) /
                                                                 (P(l_npts+1) * data->gas->meanMolecularWeight());
            sum = sum + Y(l_npts+1, k_drop[i]);
        }
        Mdotres(l_npts+1) = Mdot(l_npts+1) - ((YVmhalf(k_drop[0]) + YVmhalf(k_drop[1])) * areamhalf) /
                (1 - Y(l_npts+1, k_drop[0]) - Y(l_npts+1, k_drop[1]));
        for (size_t k = 1; k <= nsp; k++) {
            if (k != k_bath && k != k_drop[0] && k != k_drop[1]) {
                //TODO: May need to include chemical source term
                Yres(l_npts+1, k) = Y(l_npts+1, k) * Mdot(l_npts+1) + YVmhalf(k) * areamhalf;
                //Yres(1,k)=YVmhalf(k)*areamhalf;
                sum = sum + Y(l_npts+1, k);
                //if(fabs(Mdot(1))>1e-14){
                ///Yres(1,k)=innerMassFractionsData[k-1]-
                ///  Y(1,k)-
                ///  (YVmhalf(k)*areamhalf)/Mdot(1);

                //}
                //else{
                //  		//Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1];
                //  		Yres(1,k)=Y(2,k)-Y(1,k);
                //	/*The below flux boundary condition makes the
                //	 * problem more prone to diverge. How does one
                //	 * fix this?*/
                //  		//Yres(1,k)=YVmhalf(k);
                //}
            }
        }
        Yres(l_npts+1, k_bath) = ONE - sum - Y(l_npts+1, k_bath);
    }



    //p_i = 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)));     //Unit:Pa,FOR WATER
    //Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1)
    //			    / P(1) / data->gas->meanMolecularWeight();




    //Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass



    /*Energy:*/
    if (data->dirichletInner) {
        //Tres(1)=Tdot(1);
        //Tres(1)=T(1)-data->innerTemperature;
        /*Following code should be revised in the future*/
        //Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHv_l - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop));
    } else {
        //Tres(1)=Tdot(1) -
        //	(rhomhalf*areamhalfsq*lambdamhalf*(T(2)-T(1))/((psi(2)-psi(1))*mass) - Mdot(1)*dHvl)
        //	/ (massDrop * Cpl);

        //Tres(1)=Tdot(1) -
        //	(areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl)
        //	/ (massDrop * Cpl);

        /**** Boundary Condition for Temperature @ 1st grid point,
		 *** Temperature should not exceed boiling temperature of each liquid component *****/
//        T_boil = getLiquidMaxT(data->dropMole, P(1));
//        if (T(1) <= T_boil) {
//            Tres(1) = Tdot(1) -
//                      (areamhalf * lambdamhalf * (T(2) - T(1)) / (R(2) - R(1)) - Mdot(1) * dHv_l)
//                      / (massDrop * Cp_l);
//        } else {
//            //Tres(1)=T(1)-T_boil;
//            Tres(1) = Tdot(1);
//        }
        Tres(l_npts+1) = T(l_npts+1) - T(l_npts) ;

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

    /*Pressure:*/
    Pres(l_npts+1) = P(l_npts+2) - P(l_npts+1);
    /*Fill up res with governing equations at inner points:*************/
    for (size_t j = l_npts+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;
        dpsip = (psi(j + 1) - psi(j)) ;
        dpsim = (psi(j) - psi(j - 1)) ;
        dpsiav = HALF * (psi(j + 1) - psi(j - 1)) ;
        dpsipm = (psi(j + 1) - psi(j - 1)) ;
        /***********************************///

        /*evaluate various central difference coefficients*/
        cendfm = -dpsip / (dpsim * dpsipm);
        cendfc = (dpsip - dpsim) / (dpsip * dpsim);
        cendfp = dpsim / (dpsip * dpsipm);
        /**************************************************/


        /*evaluate properties at j*************************/
        setGas(data, ydata, j);
        rho = data->gas->density();        //kg/m^3
        Cpb = data->gas->cp_mass();        //J/kg/K
        Cvb = data->gas->cv_mass();        //J/kg/K
        data->gas->getNetProductionRates(wdot); //kmol/m^3
        data->gas->getEnthalpy_RT(enthalpy);    //unitless
        data->gas->getCp_R(Cp);            //unitless
        area = calc_area(R(j), &m);            //m^2

        /*evaluate properties at p*************************/
        getTransport(data, ydata, j, &rhophalf, &lambdaphalf, YVphalf);
        areaphalf = calc_area(HALF * (R(j) + R(j + 1)), &m);
        areaphalfsq = areaphalf * areaphalf;
        /**************************************************///

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

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

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

        //enthalpy formulation:
        //tranTerm = Tdot(j) - (Pdot(j)/(rho*Cpb));
        tranTerm = Tdot(j);
        sum = ZERO;
        sum1 = ZERO;
        for (size_t 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); // Qdot is forced heating to cause ignition, sum is the conversion of chemical enthalpy to sensible enthalpy
        advTerm = (Mdot(j) * (T(j) - T(j - 1)) / dpsim);
        Tres(j) = tranTerm - (Pdot(j) / (rho * Cpb))
                  + advTerm
                  - diffTerm
                  + srcTerm;

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

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

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

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

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

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

        /*Assign values evaluated at p and phalf to m
		 * and mhalf to save some cpu cost:****************/
        areamhalf = areaphalf;
        areamhalfsq = areaphalfsq;
        aream = area;
        rhom = rho;
        rhomhalf = rhophalf;
        lambdamhalf = lambdaphalf;
        for (size_t 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;
    dpsim = (psi(npts) - psi(npts - 1)) ;
    Rres(npts) = ((R(npts) - R(npts - 1)) / dpsim) - (TWO / (rhom * aream + rho * area));

    /*Energy:*/
    if (data->dirichletOuter) {
        //Tres(npts)=T(npts)-data->wallTemperature;
        Tres(npts) = Tdot(npts);
    } else {
        Tres(npts) = T(npts) - T(npts - 1);
    }

    /*Species:*/
    sum = ZERO;
    if (data->dirichletOuter) {
        for (size_t k = 1; k <= nsp; k++) {
            if (k != k_bath) {
                Yres(npts, k) = Ydot(npts, k);
                sum = sum + Y(npts, k);
            }
        }
    } else {
        for (size_t k = 1; k <= nsp; k++) {
            if (k != k_bath) {
                Yres(npts, k) = Y(npts, k) - Y(npts - 1, k);
                //Yres(npts,k)=YVmhalf(k);
                sum = sum + Y(npts, k);
            }
        }
    }
    Yres(npts, k_bath) = ONE - sum - Y(npts, k_bath);


    /*Pressure:*/
    if (data->constantPressure) {
        Pres(npts) = Pdot(npts) - data->dPdt;
    } else {
        Pres(npts) = R(npts) - data->domainLength;
        //Pres(npts)=Rdot(npts);
    }

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

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

    return (0);

}

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

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

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

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

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

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

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

void printGlobalHeader(UserData data) {
    fprintf((data->globalOutput), "%8s\t", "#Time(s)");
    //fprintf((data->globalOutput), "%15s","S_u(exp*)(m/s)");
    fprintf((data->globalOutput), "%15s", "bFlux(kg/m^2/s)");
    fprintf((data->globalOutput), "%16s", "  IsothermPos(m)");
    fprintf((data->globalOutput), "%15s\t", "Pressure(Pa)");
    fprintf((data->globalOutput), "%15s\t", "Pdot(Pa/s)");
    fprintf((data->globalOutput), "%15s\t", "gamma");
    fprintf((data->globalOutput), "%15s\t", "S_u(m/s)");
    fprintf((data->globalOutput), "%15s\t", "Tu(K)");
    fprintf((data->globalOutput), "\n");
}

//
//void printSpaceTimeRates(double t, N_Vector ydot, UserData data)
//{
//	double *ydotdata,*psidata;
//  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);
//  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);
//	for (int i = 0; i < data->npts; i++) {
//		fprintf((data->ratesOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
//		for (int j = 0; j < data->nvar; j++) {
//			fprintf((data->ratesOutput), "%15.6e\t",ydotdata[j+i*data->nvar]);
//		}
//		fprintf((data->ratesOutput), "\n");
//	}
//	fprintf((data->ratesOutput), "\n\n");
//}
//
void printGlobalVariables(double t, N_Vector *y, N_Vector *ydot, UserData data) {
    double *ydata, *ydotdata, *innerMassFractionsData, *psidata;
    innerMassFractionsData = data->innerMassFractions;

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

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


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

        setGas(data, ydata, jj);
        Tu = T(jj);
        rho = data->gas->density();
        Cpb = data->gas->cp_mass();        //J/kg/K
        Cvb = data->gas->cv_mass();        //J/kg/K
        gamma = Cpb / Cvb;

    } else {
        while (T(j) > TAvg) {
            j = j + 1;
        }
        YAvg = innerMassFractionsData[data->k_oxidizer - 1] - Y(1, data->k_oxidizer);
        while (fabs((T(data->npts - jj - 1) - T(data->npts - jj)) / T(data->npts - jj)) > 1e-08) {
            jj = jj + 1;
        }

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

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


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

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

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

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

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


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

    double *ydata;
    ydata = N_VGetArrayPointer_OpenMP(*y); //access the data stored in N_Vector* y

    for (i = 1; i <= npts; i++) {
        setGas(data, ydata, i); //set the gas state at each grid point
        rho = data->gas->density(); //get the averaged density at each grid point, Unit:kg/m^3
        data->gas->getNetProductionRates(wdot_mole); //Unit:kmol/(m^3 * s)

        for (k = 1; k <= nsp; k++) {
            MW(k) = data->gas->molecularWeight(k - 1); //Unit:kg/kmol
        }

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

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

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

    }
}


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


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

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

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

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

void floorSmallValue(UserData data, N_Vector *y) {
    double *ydata;
    ydata = N_VGetArrayPointer_OpenMP(*y);
    //double sum = 0.00;
    size_t k_bath = data->k_bath;

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

    /*Dump the error to the bath gas*/
    for (size_t i = 1; i <= data->npts; i++) {
        double sum = 0.00;
        for (size_t k = 1; k <= data->nsp; k++) {
            if (k != k_bath) {
                sum = sum + Y(i, k);
            }
        }
        Y(i, k_bath) = ONE - sum;
    }


}

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


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

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

        for (size_t i = 1; i <= data->npts; i++) {
            atolT(i) = tempTol;
            atolR(i) = radTol;
            atolP(i) = presTol;
            atolMdot(i) = mdotTol;

            for (size_t k = 1; k <= data->nsp; k++) {
                if (k != data->k_bath) {
                    atolY(i, k) = massFracTol;
                } else {
                    atolY(i, k) = bathGasTol;
                }
            }
        }
    }

}

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

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

//
////rho : density of liquid phase
//double getLiquidRho(double dropMole[],double temp,double pres){
//     std::vector<std::string> fluids;
//     fluids.push_back("Propane");
//     fluids.push_back("n-Heptane");
//
//     std::vector<double> Temperature(1,temp),Pressure(1,pres);
//     std::vector<std::string> outputs;
//     outputs.push_back("D");
//
//	 std::vector<double> moles;
//	 moles.push_back(dropMole[0]);
//	 moles.push_back(dropMole[1]);
//     double density = CoolProp::PropsSImulti(outputs,"T",Temperature,"P",Pressure,"",fluids,moles)[0][0]; //Unit:kg/m^3
//
//     return density;
//}

////Cp,l : specific heat capacity at constant pressure of liquid phase
//double getLiquidCp(double dropMole[],double temp,double pres){
//     std::vector<std::string> fluids;
//     fluids.push_back("Propane");
//     fluids.push_back("n-Heptane");
//
//     std::vector<double> Temperature(1,temp),Pressure(1,pres);
//     std::vector<std::string> outputs;
//     outputs.push_back("CPMASS");
//
//	 std::vector<double> moles;
//	 moles.push_back(dropMole[0]);
//	 moles.push_back(dropMole[1]);
//     double cp = CoolProp::PropsSImulti(outputs,"T",Temperature,"P",Pressure,"",fluids,moles)[0][0]; //Unit:J/(kg*K)
//
//     return cp;
//}

////Hv : heat of vaporization of liquid phase at constant pressure
//double getLiquidHv(double dropMole[],double temp,double pres){
//     std::vector<std::string> fluids;
//     fluids.push_back("Propane");
//     fluids.push_back("n-Heptane");
//
//     std::vector<double> Temperature(1,temp),Pressure(1,pres);
//     std::vector<std::string> outputs;
//     outputs.push_back("H");
//
//	 std::vector<double> moles;
//	 moles.push_back(dropMole[0]);
//	 moles.push_back(dropMole[1]);
//
//	 std::vector<double> state1(1,1.0);
//	 std::vector<double> state2(1,0.0);
//
//     double H_V = CoolProp::PropsSImulti(outputs,"P",Pressure,"Q",state1,"",fluids,moles)[0][0]; //Unit:J/kg
//     double H_L = CoolProp::PropsSImulti(outputs,"P",Pressure,"Q",state2,"",fluids,moles)[0][0]; //Unit:J/kg
//     double delta_H = H_V - H_L;
//
//     return delta_H;
//}

////MaxT : maximum allowed tempereture of liquid phase
////Also the boiling point of more volitile component
//double getLiquidMaxT(double dropMole[],double pres){
//     std::vector<std::string> fluids;
//     fluids.push_back("Propane");
//     //fluids.push_back("n-Heptane");
//
//     std::vector<double> Pressure(1,pres);
//     std::vector<std::string> outputs;
//     outputs.push_back("T");
//
//	double temperature = CoolProp::PropsSI(outputs[0], "P", Pressure[0],"Q",1.0,fluids[0]);
//
//    return temperature;
//}

/*function overloading for single and bi-component*/
/*calculate the liquid phase density based on droplet type and T,P,X(Y)*/
/*single component case*/
double getLiquidDensity(const double temp, const double pres, const std::vector<std::string> &composition) {
    std::vector<double> Temperature(1, temp), Pressure(1, pres);

    double density = CoolProp::PropsSI("D", "T", Temperature[0], "P", Pressure[0], composition[0]); //Unit:kg/m^3
    return density;
}

/*binary component case*/
double getLiquidDensity(const double temp, const double pres, std::vector<std::string> &composition,
                        const std::vector<double> &mole) {
    std::vector<double> Temperature(1, temp), Pressure(1, pres);
    std::vector<std::string> outputs;
    outputs.push_back("D");

    double density = CoolProp::PropsSImulti(outputs, "T", Temperature, "P", Pressure, "", composition,
                                            mole)[0][0]; //Unit:kg/m^3
    return density;
}

double getLiquidCond(const double temp, const double pres, const std::vector<std::string> &composition) {
    std::vector<double> Temperature(1, temp), Pressure(1, pres);
    double k = CoolProp::PropsSI("conductivity", "T", Temperature[0], "P", Pressure[0], composition[0]);
    return k;  //Unit:W/m/K (kg*m/s^3/K)
}

/*binary component case*/
double getLiquidCond(const double temp, const double pres, std::vector<std::string> &composition,
                     const std::vector<double> &mole) {
    std::vector<double> Temperature(1, temp), Pressure(1, pres);
    std::vector<std::string> outputs;
    outputs.push_back("conductivity");

    double k = CoolProp::PropsSImulti(outputs, "T", Temperature, "P", Pressure, "", composition,
                                      mole)[0][0];
    return k; //Unit:W/m/K (kg*m/s^3/K)
}

/*Cpb:total liquid phase specific heat capacity at constant pressure*/
double getLiquidCpb(const double temp,const double pres, const std::vector<std::string>& composition){
    std::vector<double> Temperature(1, temp), Pressure(1, pres);
    double k = CoolProp::PropsSI("Cpmass", "T", Temperature[0], "P", Pressure[0], composition[0]); //Unit:kg/m^3
    return k;  //unit:J/kg/K
}

/*binary component case*/
double getLiquidCpb(const double temp, const double pres, const std::vector<std::string> &composition,
                    const std::vector<double> &mole) {
    std::vector<double> Temperature(1, temp), Pressure(1, pres);
    std::vector<std::string> outputs;
    outputs.push_back("Cpmass");

    double k = CoolProp::PropsSImulti(outputs, "T", Temperature, "P", Pressure, "", composition,
                                      mole)[0][0]; //Unit:kg/m^3
    return k; //unit:J/kg/K
}

/*get Cp vector for binary componet droplet*/
std::vector<double> getLiquidCp(const double temp, const double pres, const std::vector<std::string> &composition){
    std::vector<double> k ;
    double k_temp;
    std::vector<double> Temperature(1, temp), Pressure(1, pres);
    for(auto & species : composition){
        k_temp = CoolProp::PropsSI("Cpmass", "T", Temperature[0], "P", Pressure[0], species);
        k.push_back(k_temp) ;
    }
    return k;
}


/*get gas phase thermal conductivity*/
double getGasCond(UserData data, double *ydata, size_t gridPoint) {
    double k;
    setGas(data, ydata, gridPoint);
    k = data->trmix->thermalConductivity();
    return k; //unit:W/(m*K)
}

/*latent heat of vaporizaiton*/
/*unit: J/kg*/
std::vector<double> getLiquidVH(const double pres, const int dropType) {
    double H, H_V, H_L;
    std::vector<double> vapheat;
    std::vector<std::string> fluids;
    std::vector<double> Pressure(1, pres);

    fluids = components(dropType);

    if (dropType == 0) {
//        fluids.push_back("nHeptane");
        H_V = CoolProp::PropsSI("H", "P", Pressure[0], "Q", 1, fluids[0]);
        H_L = CoolProp::PropsSI("H", "P", Pressure[0], "Q", 0, fluids[0]);
        H = H_V - H_L;
        vapheat.push_back(H);
    } else {
        for (size_t i = 0; i < fluids.size(); i++) {
            H_V = CoolProp::PropsSI("H", "P", Pressure[0], "Q", 1, fluids[i]);
            H_L = CoolProp::PropsSI("H", "P", Pressure[0], "Q", 0, fluids[i]);
            H = H_V - H_L;
            vapheat.push_back(H);
        }
    }
    return vapheat;  //unit:J/kg
}

/*convert mass fraction to mole fraction*/
/*size of both arrays should be identical*/
void mass2mole(const std::vector<double> &mass, std::vector<double> &mole, UserData data) {
    mole.reserve(mass.size());
    double sum = 0.0;
    double x_temp;
    for (size_t i = 0; i < data->dropMole.size(); i++) {
        sum = sum + mass[i] / data->MW[i];
    }

    x_temp = mass[0] / (data->MW[0] * sum);
    mole.push_back(x_temp);
    mole.push_back(ONE - x_temp);
}

/*liquid mass diffusivity D for binary component case*/
double getLiquidmassdiff(UserData data, double *ydata, size_t gridPoint, const double temp) {
    std::vector<double> V_ = {75.86, 162.34}; //molar volume at normal boiling point, unit: cm^3/mol
    std::vector<double> visc_, mole_, mass_, D_inf_;
    double visc_temp, D, Dinf_temp;
    std::vector<std::string> composition = components(data->dropType);

    for (size_t i = 0; i < data->dropMole.size(); i++) {
        visc_temp = 1000.0 * CoolProp::PropsSI("V", "T", temp, "Q", 0, composition[i]);  //unit: mPa*s
        visc_.push_back(visc_temp);
    }
    mole_ = getLiquidmolevec(data, ydata, gridPoint);

    for (size_t i = 0; i < data->dropMole.size(); i++) {
        if (i == 0) {
            Dinf_temp = 9.89e-8 * temp * pow(visc_[1], -0.907) * pow(V_[0], -0.45) * pow(V_[1], 0.265);
            D_inf_.push_back(Dinf_temp);
        } else {
            Dinf_temp = 9.89e-8 * temp * pow(visc_[0], -0.907) * pow(V_[1], -0.45) * pow(V_[0], 0.265);
            D_inf_.push_back(Dinf_temp);
        }
    }

    D = mole_[1] * D_inf_[0] + mole_[0] * D_inf_[1]; //unit:cm^2/s
    return (D / 10000); //unit:m^2/s
}

//double dropletmass(UserData data,double* ydata){
//    double mass = 0.0 ;
//    double rho,rhop;
//    if (data->dropType==0){
//            for (size_t i = data->l_npts-1; i>=1 ; i--) {
//                std::vector<std::string> str;
//                str.push_back("nHeptane");
//                rho = getLiquidDensity(T(i), P(i), str);
//                rhop = getLiquidDensity(T(i + 1), P(i + 1), str);
//                mass = mass + (R(i + 1) - R(i)) * (rho * calc_area(R(i), &data->metric) +
//                                                           rhop * calc_area(R(i + 1), &data->metric)) / TWO;
//            }
//    }else if(data->dropType==1){
//            for (size_t i = data->l_npts-1; i>=1 ; i--){
//                std::vector<std::string> str;
//                std::vector<double> massFrac,massFracp,moleFrac,moleFracp;
//                str.push_back("propane");
//                str.push_back("nHeptane");
//                for (size_t j = 0; j < data->dropMole.size(); j++){
//                    massFrac.push_back(Y(i,data->k_drop[j]));
//                    massFracp.push_back(Y(i+1,data->k_drop[j]));
//                }
//                /*convert mass fraction to mole fraction*/
//                mass2mole(massFrac,moleFrac,data);
//                mass2mole(massFracp,moleFracp,data);
//                rho = getLiquidDensity(T(i),P(i),str,moleFrac);
//                rhop = getLiquidDensity(T(i+1),P(i+1),str,moleFracp);
//                mass =  mass +(R(i+1)-R(i))*(rho*calc_area(R(i),&data->metric)+rhop*calc_area(R(i+1),&data->metric))/TWO;
//        }
//    }
//    return mass ;
//}


// List of fluid components per dropType
std::vector<std::string> components(int dropType) {
    if (dropType == 0) {
        return {"nHeptane"};
    } else if (dropType == 1) {
        return {"propane", "nHeptane"};
    }
    return {};
}

/*calculat droplet mass*/
double dropletmass(UserData data, double *ydata) {
    double mass = 0.0;
    double rho, rhop;

    for (size_t i = data->l_npts - 1; i >= 1; i--) {
        std::vector<std::string> composition = components(data->dropType);

        if (data->dropType == 0) {
            rho = getLiquidDensity(T(i), P(i), composition);
            rhop = getLiquidDensity(T(i + 1), P(i + 1), composition);

        } else if (data->dropType == 1) {
            std::vector<double> moleFrac, moleFracp;
            moleFrac = getLiquidmolevec(data, ydata, i);
            moleFracp = getLiquidmolevec(data, ydata, i);

            rho = getLiquidDensity(T(i), P(i), composition, moleFrac);
            rhop = getLiquidDensity(T(i + 1), P(i + 1), composition, moleFracp);
        }

        double area = calc_area(R(i), &data->metric);
        double areap = calc_area(R(i + 1), &data->metric);
        mass += (R(i + 1) - R(i)) * (rho * area + rhop * areap) / TWO;
    }
    return mass;
}

std::vector<double> getLiquidmassvec(UserData data, double *ydata, int gridPoint) {
    std::vector<double> mass_;
    mass_.push_back(Y(gridPoint, data->k_drop[0]));
    mass_.push_back(Y(gridPoint, data->k_drop[1]));
    return mass_;
}

std::vector<double> getLiquidmolevec(UserData data, double *ydata, int gridPoint) {
    std::vector<double> mass_, mole_;
    mass_ = getLiquidmassvec(data, ydata, gridPoint);
    mass2mole(mass_, mole_, data);
    return mole_;
}

/*return droplet species vapor pressure at interface, unit:Pa*/
std::vector<double> getVapPressure(UserData data, double* ydata,int gridPoint,const std::vector<double> mole_){
    std::vector<double> vapPres;
    vapPres.reserve(mole_.size()) ;
    std::vector<double> gamma_;
    getGamma(mole_,gamma_) ;
    //propane
    double p1 = 1.0e5 * pow(10, 4.53678 - (1149.36 / (T(gridPoint) + 24.906))) * mole_[0] *
                gamma_[0];
    //heptane
    double p2 = 1.0e5 * pow(10, 4.02832 - (1268.636 / (T(gridPoint) - 56.199))) * mole_[1] *
                gamma_[1];
    vapPres.push_back(p1);
    vapPres.push_back(p2);
    return vapPres;
}

void printIddata(UserData data, double *iddata) {
    FILE *file;
    file = fopen("id.dat", "w");
    fprintf(file, "%15s\t%15s\t", "Rid", "Tid");
    for (size_t k = 1; k <= data->nsp; k++) {
        fprintf(file, "%15s\t", data->gas->speciesName(k - 1).c_str());
    }
    fprintf(file, "%15s\t", "Pid");
    fprintf(file, "%15s\n", "Mdotid");

    for (size_t i = 1; i <= data->npts; i++) {
        fprintf(file, "%15lu\t", (size_t) Rid(i));
        fprintf(file, "%15lu\t", (size_t) Tid(i));
        for (size_t k = 1; k <= data->nsp; k++) {
            fprintf(file, "%15lu\t", (size_t) Yid(i, k));
        }
        fprintf(file, "%15lu\t", (size_t) Pid(i));
        fprintf(file, "%15lu\n", (size_t) Mdotid(i));
    }

    fclose(file);
}

void printPsidata(UserData data, double* psidata) {
    FILE *file;
    file = fopen("psi.dat", "w");
//    fprintf(file, "%15s\t%15s\t", "Rid", "Tid");
    for (size_t j = 1; j <= data->npts; j++) {
        fprintf(file, "%lu\t", j);
    }
    fprintf(file, "\n");

    for (size_t j = 1; j <= data->npts; j++) {
//        fprintf(file, "%15lu\t", (size_t) Rid(i));
//        fprintf(file, "%15lu\t", (size_t) Tid(i));
//        for (size_t k = 1; k <= data->nsp; k++) {
//            fprintf(file, "%15lu\t", (size_t) Yid(i, k));
//        }
//        fprintf(file, "%15lu\t", (size_t) Pid(i));
//        fprintf(file, "%15lu\n", (size_t) Mdotid(i));
        fprintf(file,"%15.6e\t", psi(j));
    }
    fprintf(file, "\n");
    fclose(file);
}