/*			       24.APR.2003		      rungkutcs.cc

*/

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

#include "corso.h"

void RungKutCS(double *x,double *vstart,int nVar,double xMin, double xMax, 
	       int nStep,double *xp,double **yp,int DeltaSave,int *nSaved,
	       void (*derivs)(double *x,double *, double *))
{
  int i,k,nSav;
  double h,h6,hh,x0,xh;
  double *dv,*dym,*dyt,*v,*vout,*yt;

  // .......................................................... allocate arrays
  dv=new double [nVar];
  dym=new double [nVar];
  dyt=new double [nVar];
  v=new double [nVar];
  vout=new double [nVar];
  yt=new double [nVar];
  // ........................................................... initial values
  for (i=0;i<nVar;i++) v[i]=yp[i][0]=vstart[i];
  x[0]=xp[0]=xMin;
  nSav=1;
  // ..................................................................... step
  h=(xMax-xMin)/nStep;
  hh=h/2.0;
  h6=h/6.0;
  // ......................................................... Runge-Kutta loop
  for (k=0;k<nStep;k++) 
    { (*derivs)(x,v,dv);
      // ............................................................. midpoint
      x[0]+=hh;
      // ........................................................... first step
      for (i=0;i<nVar;i++) yt[i]=v[i]+hh*dv[i];
      (*derivs)(x,yt,dyt);
      // .......................................................... second step
      for (i=0;i<nVar;i++) yt[i]=v[i]+hh*dyt[i];
      (*derivs)(x,yt,dym);
      // ........................................................... third step
      for (i=0;i<nVar;i++) 
	{ yt[i]=v[i]+h*dym[i];
	  dym[i] += dyt[i];
	}
      // .......................................................... fourth step
      x[0]+=hh;
      (*derivs)(x,yt,dyt);
      // .......................................................... fourth step
      for (i=0;i<nVar;i++) vout[i]=v[i]+h6*(dv[i]+dyt[i]+2.0*dym[i]);
      // ......................................................................
      for (i=0;i<nVar;i++) v[i]=vout[i];
      if (((k+1)%DeltaSave)==0)
	{ xp[nSav]=x[0];
          for (i=0;i<nVar;i++) yp[i][nSav]=v[i];
	  nSav++;
	}
    }
  // ......................................... number of points saved in array
  *nSaved=nSav;
  // ............................................................. free memory
  delete [] yt;
  delete [] vout;
  delete [] v;
  delete [] dyt;
  delete [] dym;
  delete [] dv;
}