#include <stdio.h>

#include "integrate.h"
#include "nrutil.h"
#include "input.h"

/* --- Global variables ---------------------------------------------------- */

double *xp, **yp;
double dxsav;
int kmax, kount;
struct Input input;

/* --- Handles the numerical integration. Many of the parameters are set here
       instead of input parameters.
       Note that the routine used is Stifbs. --      ----------------------- */


void integrate(double zstart, double ystart[], int nv, int *nout, double *zout, double **yout)

{
  double z1, z2;
  double htry, hmin;
  double eps;
  int nok, nbad, i;

  void (*p)(double z, double y[], double fvec[]);

  /* z1, z2 start and end values of integration
   * htry = stepsize to try
   * hmin = minimum stepsize
   * eps = accuracy
   * nok, nbad = number of ok/bad steps
   * p = pointer to function type derivs
   * kmax = maximum integer to store steps (not max number of steps)
   * kount = actual max integer of steps stored
   * dxsav = approximate interval to save y[i] values
   * xp, yp are the x and y values saved at intervals dxsav
   */


  p = derivs;

  z1 = input.zstart; z2 = input.zend; 
  htry = 1.0e-05; hmin = 1.0e-06; eps = 1.0e-07;
  kmax = 100000; dxsav = 50.0;  

  xp = dvector(1,kmax);
  yp = dmatrix(1,nv,1,kmax);

  odeint(ystart,nv,z1,z2,eps,htry,hmin,&nok,&nbad,p); 

  *nout = kount;
  for (i=1; i<=kount; i++) {
    zout[i] = xp[i];
    yout[1][i] = yp[1][i]; yout[2][i] = yp[2][i]; yout[3][i] = yp[3][i];
  }

  free_dvector(xp,1,kmax);
  free_dmatrix(yp,1,nv,1,kmax);

}