/*			       03.MAY.2003			  rungestep.cc

*/

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

#include "corso.h"

void RungeStep(double *x,double *y,int nVar,double step,int FirstMidLast,
	       void (*derivs)(double *x,double *y, double *dydx))
{
  int i,k,nSav;
  static double h,h6,hh,x0,xh;
  static double *dydx_0,*dydx_1,*dydx_2,*dydx_3,*yt;

  // .......................................... if first step allocate arrays
  if (FirstMidLast==0)
    { dydx_0=new double [nVar];
      dydx_1=new double [nVar];
      dydx_2=new double [nVar];
      dydx_3=new double [nVar];
      yt=new double [nVar];
      // ........................................................ step values
      h=step;
      hh=h/2.0;
      h6=h/6.0;
    }
  // .......................................... else if last stap free memory
  else if (FirstMidLast==2)
    { delete [] yt;
      delete [] dydx_3;
      delete [] dydx_2;
      delete [] dydx_1;
      delete [] dydx_0;
      return;
    }
  // .................................................... initial derivatives
  (*derivs)(x,y,dydx_0);
  // ........................................... first estimate of midpoint y
  for (i=0;i<nVar;i++) yt[i]=y[i]+hh*dydx_0[i];
  // ............................................................. midpoint x
  x[0]+=hh;
  // ................................................... midpoint derivatives
  (*derivs)(x,yt,dydx_1);
  // .......................................... second estimate of midpoint y
  for (i=0;i<nVar;i++) yt[i]=y[i]+hh*dydx_1[i];
  // ................................ second estimate of midpoint derivatives
  (*derivs)(x,yt,dydx_2);
  // ......................................... first estimate of last point y
  for (i=0;i<nVar;i++) yt[i]=y[i]+h*dydx_2[i];
  // ................................................. last point derivatives
  x[0]+=hh;
  // ..................................... estimate of last point derivatives
  (*derivs)(x,yt,dydx_3);
  // ......................................... final estimate of last point y
  for (i=0;i<nVar;i++) 
    y[i]+=h6*(dydx_0[i]+2.0*(dydx_1[i]+dydx_2[i])+dydx_3[i]);
}