#include <math.h>
#include <stdlib.h>
#include "rk4e.h"

#define	max(a, b)	((a >= b)? a : b)

/*------------------------------------------------------------------------------
  Take one error-controlled step of integration of a set of equations
     dy
     -- = f[x, y]
     dx
  starting from xin, yin.

  The stepsize dx is adjusted by factors of 2 from the input stepsize
  so that all y[i] remain within a desired relative or absolute error:
     error in y[i] <= max(relerr * |y[i]|, abserr) * erry[i].
  Note that a y[i] passes the error test if it passes
  EITHER the relative OR the absolute test.
  The error test is skipped if erry[i] is negative.

  The actual stepsize used is never smaller than the input stepsize dx.
  If the input stepsize fails to achieve the desired error,
  then the step is rejected, the stepsize is halved,
  and the step is repeated, iteratively until the error test passes.

  If the error test passes with the input stepsize dx,
  then the step is accepted.
  If the error test passes by a sufficient margin,
  then the output stepsize dx is set to double the input stepsize dx.
  However, the actual stepsize taken is still the input stepsize dx,
  not the doubled output stepsize.  The output stepsize dx is
  to be considered as the recommended stepsize for the next step.

  Typically, one might set erry[i] = 1 for all y[i],
  so that all y[i] are treated equally.
  To suppress error checking of a y[i], set its erry[i] = -1.

  As an example of relative and absolute error, one might set, say,
  relerr = 1e-14, abserr = 1e-30.
  Then a y[i] passes the error test if it passes:
     the relative error test if y[i] is large, |y[i]| > 10^-16,
  or the absolute error test if y[i] is small, |y[i]| < 10^-16.

  Normally the relative error relerr (or more generally, relerr * erry[i])
  should be larger than the machine precision
  (about 10^-16 for standard 64-bit doubles).
  If the relative error relerr is set to machine precision or smaller,
  then the relative error test will pass only if the estimated error
  is zero to machine precision, which might be problematic
  since the estimation of errors is uncertain at machine precision.
  If the relative error is zero, then the absolute error test
  will always take precedence.

  The absolute error should be set equal to the value below which
  a y[i] is to be considered effectively equal to zero.
  Is y[i] = 10^-300 significant?
  If so, then it may be fine to set abserr = 0.
  In the above example where abserr = 1e-30, y[i] is considered
  effectively zero if |y[i]| <= 1e-30.

  The routine returns an error if the stepsize dx is zero.

   Input: deriv = derivative function.
	  xin = initial value of x.
	  ny = length of y array.
	  yin = initial values of y array.
	  relerr = desired relative error.
	  abserr = desired absolute error.
	  erry = desired error in each y relative to the overall error;
		 if erry[i] < 0, then the error check for y[i] is skipped.
  Output: xout = final value of x.
	  yout = final values of y array.
  Input/Output: dx = suggested stepsize on input;
		     recommended stepsize on output,
		     always a power of 2 times the suggested stepsize.

  Return value: 0 = ok;
		1 = stepsize went to zero.
*/
int rk4e(derivf deriv, double xin, double *xout, int ny, double *yin, double *yout, double *dx, double relerr, double abserr, double *erry, void *pass)
{
/* factor used in estimating error;
   15 = 2^N - 1 where N is the integration order */
#define	ERR_ESTIMATE_FACTOR	15.

/* double stepsize if err is less than ERR_DOUBLE_STEPSIZE;
   1/32 = 2^[-(N+1)] where N is the integration order */
#define	ERR_DOUBLE_STEPSIZE	(1. / 32.)

    /* allocate memory for temporary arrays */
    <your code>

    /* initialize to no error condition */
    ier = 0;

    /* whether integration step succeeded */
    done = 0;

    /* number of times the stepsize has been halved */
    nhalved = 0;

    /* keep going until the error is within desired limits */
    while (!done) {
	/* check for zero stepsize */
	if (*dx == 0.) {
	    ier = 1;
	    break;
	}

	/* one double-sized step */
	<your code>

	/* two normal-sized steps */
	<your code>

	/* estimate of error in each y */
	<your code>

	/* worst error */
	<your code>

	/* worst error exceeds desired */
	if (<your code>) {
	    <your code>

	/* worst error within desired */
	} else {
	    <your code>

	}
    }

    /* improved estimate of y from combination of the two estimates */
    <your code>

    /* free memory for allocated arrays */
    <your code>

    return(ier);
}