#include <math.h>
#include "kutta.h"
#define MAXSTP 10000
#define TINY 1.0e-30
#define MAXRUNG 13
#define NEQ 3
static double rwork[(MAXRUNG+2)*NEQ];


int kmax, kount;
double *xp, **yp, dxsav;

 void f(double x, double *y, double *dydx)
{

     dydx[0] = y[1]*y[2];
     dydx[1] = -y[0]*y[2];
     dydx[2] = -0.51*y[0]*y[1];
     return;
}  

/* void f(double x, double *y, double *dydx)
{
     dydx[0] = -10.0 * y[0];
     return;
}  */
     

/* user storage for intermediate results.              */

/* Runge-kutta driver with adaptive stepsize control.  Integration start 
   with ystart[1..nvar] form x1 to x2 with accuracy  eps.  h1 is set as
   a guessed first stepsize.  hmin as the minimum allowed stepsize. 
   On output nok and nbad are the number of good anf bad steps taken
   and ystart is replaced by values at the end of the integration interval.
   f is the user defined routine to calculate the right-hand side of 
   the differential equation.  while rkqs is the name of the stepper 
   routine to be used.             */

void odeint(double *ystart, int nvar, double x1, double x2, double eps,
            double h1, double hmin, int *nok, int *nbad,
            void (*f)(double, double *, double *),
            void (*rkqs)(double *, double *, int, double *, double,
                double, double *, double *, double *, 
                void (*)(double, double *, double *)))
{
    int nstp, i;
    double xsav, x, hnext, hdid, h;
    double *yscal, *y, *dydx;
   printf("inside odeint: eps:%lf, h1:%lf, hmin:%lf, x1:%f, x2:%f, nvar:%d\n",
         eps, h1, hmin, x1, x2, nvar);
  for(i=0;i<nvar;i++) printf("ystart:%lf\n", ystart[i]);
    yscal = VectorAlloc(nvar);
    y = VectorAlloc(nvar);
    dydx = VectorAlloc(nvar);
    x = x1;
    h = h1*(x2-x1)/fabs(x2-x1);  
  printf("inside odeint: h:%lf\n", h);

    *nok = (*nbad) = kount = 0;
    for(i=0; i<nvar; i++) y[i] = ystart[i];
    if(kmax > 0) xsav = x-dxsav*2.0;
    for(nstp = 1; nstp <= MAXSTP; nstp++) {
         (*f)(x,y,dydx);
         for(i=0; i<nvar; i++) {
               yscal[i] = fabs(y[i]) + fabs(dydx[i]*h) + TINY;   
  /*             yscal[i] = 1.0;                     */
            }
         if (kmax>0 && kount<kmax-1 && fabs(x-xsav)>fabs(dxsav)) {
               xp[kount] = x;
               for(i=0;i<nvar;i++) yp[kount][i]=y[i];
               kount++;
               xsav = x;
           }
         if((x+h-x2)*(x+h-x1)>0.0) h=x2-x;
         (*rkqs)(y, dydx, nvar, &x, h, eps, yscal, &hdid, &hnext, f);
         if(hdid == h) ++(*nok); else ++(*nbad);
         if((x-x2)*(x2-x1) >= 0.0) {
             for(i=0;i<nvar;i++) ystart[i] = y[i];
             if(kmax) {
               xp[++kount] = x;
               for(i=0; i<nvar; i++) yp[kount][i] = y[i];
              }
             VectorFree(nvar, dydx);
             VectorFree(nvar, y);
             VectorFree(nvar, yscal);
             return;
          }
          if(fabs(hnext) <= hmin) rkerror("Step size too small in odeint");
          h=hnext;
    }
    rkerror("Too many steps in routine odeint");
}


void main() {

   double ystart[1], x1, x2, eps, epsmin, h1, hmin, inc;
   int i, j, n, nbad, nok;


   x1 = 0.0;
   ystart[0] = 0.0;
   ystart[1] = 1.0;
   ystart[2] = 1.0;

   epsmin = FindMachEps();
   eps = epsmin * 10000000;
   h1 = eps * 1000000000;
   hmin = eps;
   n = NEQ;
  
   kmax = 15;
   dxsav = 1.8626;
   xp = VectorAlloc(kmax);
   yp = MatrixAlloc(kmax, kmax);

   
      x2 = x1 + 50*dxsav;
      rkstart(NEQ, rwork);
      odeint(ystart, n, x1, x2, eps, h1, hmin, &nok, &nbad,
             &f, &rkqs);

      printf("Integration result:   \n");
   for(i=0;i<kount;i++) {
      printf("x      y1        y2        y3\n");
      printf("%f    ", xp[i]);
       for(j=0;j<n; j++)
           printf("%25.22lf     ", yp[i][j]);
      printf("\n");
      printf("ok operation:%d, bad operation:%d\n", nok, nbad);
     }

}