#ifndef MACRO_DEF
#define MACRO_DEF
//#define SUNDIALS_DOUBLE_PRECISION 1
//#define SUNDIALS_SINGLE_PRECISION 1

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

/* In order to keep begin the index numbers from 1 instead of 0, we define
 * macros here. Also, we define macros to ease referencing various variables in
 * the sundials nvector.
 */
#define psi(i)   psidata[i-1]

#define T(i)   ydata[((i-1)*data->nvar)+data->nt]
#define Y(i,k) ydata[((i-1)*data->nvar)+data->ny+k-1]
#define R(i)   ydata[((i-1)*data->nvar)+data->nr]
#define P(i)   ydata[((i-1)*data->nvar)+data->np]

#define Tdot(i)   ydotdata[((i-1)*data->nvar)+data->nt]
#define Ydot(i,k) ydotdata[((i-1)*data->nvar)+data->ny+k-1]
#define Rdot(i)   ydotdata[((i-1)*data->nvar)+data->nr]
#define Pdot(i)   ydotdata[((i-1)*data->nvar)+data->np]

#define Tres(i)   resdata[((i-1)*data->nvar)+data->nt]
#define Yres(i,k) resdata[((i-1)*data->nvar)+data->ny+k-1]
#define Rres(i)   resdata[((i-1)*data->nvar)+data->nr]
#define Pres(i)   resdata[((i-1)*data->nvar)+data->np]

#define Tid(i)   iddata[((i-1)*data->nvar)+data->nt]
#define Yid(i,k)   iddata[((i-1)*data->nvar)+data->ny+k-1]
#define Rid(i)   iddata[((i-1)*data->nvar)+data->nr]
#define Pid(i)   iddata[((i-1)*data->nvar)+data->np]

#define Yav(i) Yav[i-1]
#define YAvg(i) YAvg[i-1]
#define YVmhalf(i) YVmhalf[i-1]
#define YVphalf(i) YVphalf[i-1]
#define X(i) X[i-1]
#define Xp(i) Xp[i-1]
#define Xgradhalf(i) Xgradhalf[i-1]
#define XLeft(i) XLeft[i-1]
#define XRight(i) XRight[i-1]
#define gradX(i) gradX[i-1]
#define wdot(i) wdot[i-1]
#define enthalpy(i) enthalpy[i-1]
#define energy(i) energy[i-1]
#define Cp(i) Cp[i-1]

#define atolT(i)   atolvdata[((i-1)*data->nvar)+data->nt]
#define atolY(i,k)   atolvdata[((i-1)*data->nvar)+data->ny+k-1]
#define atolR(i)   atolvdata[((i-1)*data->nvar)+data->nr]
#define atolP(i)   atolvdata[((i-1)*data->nvar)+data->np]


#define constraintsY(i,k)   constraintsdata[((i-1)*data->nvar)+data->ny+k-1]

#endif