/* Note: all floats have been changed to doubles in this function, and
 * all memory allocations using the function vector() have been changed
 * to dvector() which allocates doubles instead of floats. 
 *
 * Portions Copyright 1993 Massachusetts Institute of Technology
 */

#include <math.h>
#include <stdio.h>

#define MAXSTP 1000000000
#define TINY 1.0e-8

extern void ((*hd_energy[]) (double));  /* Pointer to hd energy functions           */
extern void calculate_config_energy();  /* Second energy function                   */
extern int model;                       /* Number of hd model in use                */
extern int calc_energy;                 /* Flag to calculate energy                 */
extern int ideal_flag;                  /* Flag indicating ideal hd model           */
extern double v[];                      /* vector containing all present data       */
extern double **output_data;            /* matrix to store output data              */
extern double *time;                    /* vector to store output time values       */
extern int index[];                     /* vector listing the indexes to store      */
extern int num_selected_output_var;     /* the number of output variables to store  */
extern double step;                     /* stepsize interval to store data          */
extern int debug;                       /* debugging variable                       */
int kount=0;                            /* global counter variable                  */
int points_to_store;                    /* maximum number of steps to store         */

void odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,ptr_derivs,ptr_rkqc)
double ystart[],x1,x2,eps,h1,hmin;
int nvar,nok[],nbad[];
void (*ptr_derivs)(double,double *,double *);
void (*ptr_rkqc)(double *,double *,int,double *,double,
		 double,double *,double *,double *,void (*)());
{
	int nstp,i,j;
	double xsav,x,hnext,hdid=0.0,h;
	double *yscal,*y,*dydx,*dvector(),**dmatrix();
	void nrerror(),free_dvector();

	/* Calculate the number of points_to_store from the step size
	 * interval for data storage and the initial and final times. */
	points_to_store = (int) ((x2 - x1)/step);

	/* Dynamically allocate space for the output data based on
	 * the points_to_store variable. */
	time = dvector(1,points_to_store);
	output_data = dmatrix(1,num_selected_output_var,1,points_to_store);

	/* Print out simulation time header. */
	printf("\n Destination Time: %lf\n", x2);
	printf("\n Current Time: ");

	yscal=dvector(1,nvar);
	y=dvector(1,nvar);
	dydx=dvector(1,nvar);
	x=x1;
	h=(x2 > x1) ? fabs(h1) : -fabs(h1);
	nok[0] = nbad[0] = kount = 0;
	for (i=1;i<=nvar;i++) y[i]=ystart[i];
	if (points_to_store > 0) xsav=x-step*2.0;
	for (nstp=1;nstp<=MAXSTP;nstp++) {

	  /* After a successful step returns accurate values for the state
	   * vector, call the function which calculates the derivatives in order
	   * to: (1) determine the quality-control scaling vector, and 
	   * (2) calculate new values for the first point of the next time
	   * step.  With these new values stored, then calculate energy if 
	   * desired, and store the results in the output vector.  */
	        (*ptr_derivs)(x,y,dydx);
	        if(calc_energy)
		  {
		    if(!ideal_flag)
		      (*(hd_energy[model]))(hdid);
		    calculate_config_energy(hdid);
		  }
			
		for (i=1;i<=nvar;i++)
		  yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
		if (points_to_store > 0) {
			if (fabs(x-xsav) > fabs(step)) {
				if (kount < points_to_store-1) {
					time[(++kount)]=x;
					for (i=1;i<=num_selected_output_var;i++) 
					  output_data[i][kount] = v[(index[(i-1)])];

					/* Now set the xsav variable to the
					 * next time step. */
					xsav = x1 + (kount-1.0)*step;
				      }
			}
		}
		if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;

                (*ptr_rkqc)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,ptr_derivs);

		if (hdid == h) ++(nok[0]); else ++(nbad[0]);
		if ((x-x2)*(x2-x1) >= 0.0) {
			for (i=1;i<=nvar;i++) ystart[i]=y[i];
			if (points_to_store) {
				time[(++kount)]=x;
				for (i=1;i<=num_selected_output_var;i++) 
					  output_data[i][kount] = v[(index[(i-1)])];
			      }
			free_dvector(dydx,1,nvar);
			free_dvector(y,1,nvar);
			free_dvector(yscal,1,nvar);
			return;
		}
		if (fabs(hnext) <= hmin) nrerror("Step size too small in ODEINT");

		/* If the next step size is larger than the data storage step
		 * size, step, then set hnext to much less than the data storage 
		 * step size so that no output data points are skipped and so that 
		 * the output data points are at roughly equally spaced time 
		 * intervals for an FFT. */
	        if (fabs(hnext) > (step/2.0)) hnext = (step/2.0);

		h=hnext;
	}
	nrerror("Too many steps in routine ODEINT");
}

#undef MAXSTP
#undef TINY