/* ici tu crees ta condition initiale */
  
  
  double reltol; /* relative tolerance */
  double ca; /* current crank angle degree */
  double t=t0; /* current time (in s) */
  double tout=t0+dt; /* next time step (in s) */
  
  UserData data=NULL;
  N_Vector abstol=NULL;
  void *cvode_mem=NULL;    
  int flag;
    
  /* Create serial vector of length neq for abstol */
  abstol = N_VNew_Serial(neq); 
  if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(EXIT_FAILURE);

  data = MALLOC(sizeof *data);  /* Allocate data memory */
  if(check_flag((void *)data, "malloc", 2)) return(EXIT_FAILURE);

/* cette structure contient tout plein de donnees intervenant dans la fonction a integrer */
  data=CreateUserData(e,mec,temp,m,wosch);  

  /* Set the scalar relative tolerance */
  reltol = rtol;

   /* 
     Call CVodeCreate to create the solver memory :
     
     CV_BDF     specifies the Backward Differentiation Formula
     CV_NEWTON  specifies a Newton iteration

     A pointer to the integrator problem memory is returned and stored in
     cvode_mem.
    */

  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(EXIT_FAILURE);

/* j'alloue toute la memoire qu'il me faut :
Chemistry est la fonction que j'integre (c'est la fonction f de xdot = f(t,x))
atol_, abstol, reltol sont des niveaux d'erreurs absolues et relatives que je me donne
t0 est le temps initial (en general t0 = 0)
y est ma condition initiale
data est une structure qui contient tout ce que je veux (par exemple des constantes intervenant dans la fonction f
neq est le nombre d'equations (donc y est de longueur neq)
dt est le pas de temps que je me fixe
*/
  int sma=SolverMemoryAllocation(cvode_mem,Chemistry,atol_,abstol,reltol,t0,y,data,neq,dt);
  if(sma==1)
  {
    fprintf(stderr,"%s %d : Error : Memory Solver Allocation failed\n",__FILE__,__LINE__);
    exit(EXIT_FAILURE);
  }
   

  /****************************************************************************
   *                                                                          *
   *                             Integration                                  *
   *                                                                          *
   ****************************************************************************/
     
  fprintf(stdout,"\nDifferential System Integration\n");
  
  
  /* writing the initial condition on the hard disk */
  PrintOutput(fileOutput,t0,y,e,mec,m);
     
  const double tfin=60./e->rpm; /* final time of integration*/  

  /* In loop, call CVode, test for error and g results on the hard disk */
  while(tout<=tfin)
  {    
    /* ici on integre entre t et tout. A la fin de l'integration, t vaut tout */
	flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
    if (check_flag(&flag, "CVode", 1))
    {
      fprintf(stderr,"%s:%d : Error : CVode failed at t = %g\nExit program\n",__FILE__,__LINE__,t);
      return(EXIT_FAILURE);     
    }
      
	/* par securite j'impose t = tout (probleme d'arrondi, de propagation d'erreur etc) */  
    t=tout;
    
    /* incrementation time */
    ca=e->rpm6*t-180.;
	
	/* j'incremente le pas de temps. Comme le problem est raide, je prends un dt different. Je m'explique :
	je sais que la raideur a lieu entre deux instants t1 et t2 avec t1 < t2. Pour gagner en nombre d'iteration (et
	donc en temps de calcul) je dis que pour 0 <= t <= t1 et pour t2 < t < tfin dt = dtmax et pour 
	t1 <= t <= t2 dt = dtmin ou dtmin < dtmax sont 2 valeurs que je me donne
	*/
    dt=(ca<ss->cadmin || ca>ss->cadmax) ? dtmax : dtmin;
   
    /* set maximum step size */
    flag = CVodeSetMaxStep(cvode_mem,dt);
    if (check_flag(&flag, "CVodeSetMaxStep", 1)) return 1;
    
	/* incrementation du temps : j'ai du plus haut que t valait tout. On increment tout de dt puis on va 
	integrer entre t et tout
	*/
    tout += dt;
    
    /* write outputs on the hard disk ==> ainsi j'ai ma solution a chaque pas de temps */
    PrintOutput(fileOutput,t,y,e,mec,m);
      
  } /* end of while(tout<=tfin) */
  
  /* -------------------------------------------------------- */
  

/* Cette fonction est un copier - coller de CVode (car on a acces au code source !) */
static int check_flag(void *flagvalue, char *funcname, int opt)
{
    int *errflag;

    /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
    if (opt == 0 && flagvalue == NULL) {
	fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
		funcname);
	return(1); }

    /* Check if flag < 0 */
    else if (opt == 1) {
	errflag = (int *) flagvalue;
	if (*errflag < 0) {
	    fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
		    funcname, *errflag);
	    return(1); }}

    /* Check if function returned NULL pointer - no memory allocated */
    else if (opt == 2 && flagvalue == NULL) {
	fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
		funcname);
	return(1); }

    return(0);
}

static int SolverMemoryAllocation(void * cvode_mem,CVRhsFn fun,double atol_,N_Vector abstol,double
reltol,double t0,N_Vector y,UserData data,unsigned int neq,double dt)
{
/*  Allocates memory for the solver

    Returns 1 if allocation failed, 0 else.
  
    cvode_mem : the solver
    fun : name of the function f  y'=f(t,y)
    atol_ : absolute tolerance
    abstol : vector of the absolute tolerance for all the unknows
    reltol : relative tolerance
    t0 : intial time
    y : initial condition
    data : user data
    neq : number of equations
    dt : time step
 */
 
 unsigned int i;
 int flag;
 
 /* Set the vector absolute tolerance */
  for(i=0;i<neq;++i) NV_Ith_S(abstol,i)=atol_;
  
   /* 
    Call CVodeMalloc to initialize the integrator memory : 
       
    cvode_mem is the pointer to the integrator memory returned by
              CVodeCreate
    f         is the user's right hand side function in y'=f(t,y)
    t0        is the initial time
    y         is the initial dependent variable vector
    CV_SV     specifies scalar relative and vector absolute tolerances
    &reltol   is a pointer to the scalar relative tolerance
    abstol    is the absolute tolerance vector
   */

  flag = CVodeMalloc(cvode_mem, fun, t0, y, CV_SV, reltol, abstol);
  if (check_flag(&flag, "CVodeMalloc", 1)) return 1;

  /* Set the pointer to user-defined data */
  flag = CVodeSetFdata(cvode_mem, data);
  if(check_flag(&flag, "CVodeSetFdata", 1)) return 1;

  /* Call CVDense to specify the CVDENSE dense linear solver */
  flag = CVDense(cvode_mem, neq);
  if (check_flag(&flag, "CVDense", 1)) return 1;

  /* set maximum step size */
  flag = CVodeSetMaxStep(cvode_mem,dt);
  if (check_flag(&flag, "CVodeSetMaxStep", 1)) return 1;

  /* Set maximum order of the BDF method */
  flag = CVodeSetMaxOrd(cvode_mem,2);
  if (check_flag(&flag, "CVodeSetMaxOrd", 1)) return 1;
    
  /* Set maximum number of steps to be taken by the solver to reach the next
     output time */
  //flag=CVodeSetMaxNumSteps(cvode_mem,3000);
  //if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1;

  /* Set the Jacobian routine to Jac (user-supplied) and
     the pointer to the user-defined block data. */
  flag = CVDenseSetJacFn(cvode_mem, NULL, data);
  if (check_flag(&flag, "CVDenseSetJacFn", 1)) return 1;
  
  return 0;
}