#include "rkf45.h"

real_8
      k;              /* parameter for f() */

main()

   /* Test n-dimensional Runge Kutta Fehlberg integrator. */

{
   real_8
      x_init,         /* initial value of independent variable */
      x_fin,          /* final   value of independent variable */
      x_low,          /* low  intermediate value of indep. variable */
      x_high,         /* high intermediate value of indep. variable */
      y_init[NEQN],   /* initial value of dep. variable */
      y_low[NEQN],    /* initial value of dep. variable for step */
      y_high[NEQN],   /* value of dep. var. after step */
      y_fin[NEQN],    /* final answer, if integration was successful */
      y_prime[NEQN],  /* storage space for derivatives */
      abs_err[NEQN],  /* absolute error tolerated on integration */
      rel_err[NEQN],  /* relative error tolerated on integration */
      epsilon[NEQN],  /* total tolerance on integration */
      step,           /* step size given to rkf45 */
      new_step,       /* new step recommended by rkf45 */
      min_step;       /* smallest tolerable step size */
   boolean
      finished;       /* status of stepping */
   integer
      n,              /* number of dimensions */
      fail,           /* status of rkf45 */
      m,              /* used for computing starting step size */
      i;              /* dummy index for array indexing */
   boolean
      f2();           /* function which calculates derivatives */
      f3();           /* function which calculates derivatives */

   /* Initial conditions */
/* n = 1;                single first  order equation for f2 */
   n = 2;             /* single second order equation for f3 */

   k = 6.283185308;
   x_init = 0.0;
   y_init[0] = 1.0;   /* y(0) */
   y_init[1] = 0.0;   /* generate cosine */
   x_fin  = 10.0;

   /* Initialization for integration itself */
   x_high = x_init;
   for (i = 0; i < n; ++i)
      y_high[i] = y_init[i];
   m = 100;
   new_step = (x_fin - x_init) / m;
   rel_err[0] = 1.0e-3;
   abs_err[0] = 1.0e-2;
   rel_err[1] = 5.0 * rel_err[0];
   abs_err[1] = 5.0 * abs_err[0];
   min_step = 1.0e-5;
   fail = 0;
   finished = false;

   printf("Start tstrkf45\n");
   /* Stepping loop with automatic step size selection. */
   while ((fail == 0)  && ! finished)
   {
      step = new_step;
      x_low = x_high;
      x_high = x_low + step;
      for (i= 0; i < n; ++i)
      {
         y_low[i] = y_high[i];
         epsilon[i] = abs_err[i] + fabs(y_low[i]) * rel_err[i];
      }

      if (((x_high - x_fin) / step) > 0)
      {
         x_high = x_fin;
         finished = true;
      }
      rkf45(n, f3, y_prime, x_low, y_low, x_high, epsilon, step,
             min_step, y_high, &new_step, &fail);

      /* record progress of integration  
      printf("%14.4e %8.3f %14.4e %14.4e %14.4e\n", 
             (x_high - x_low), x_high, y_high[0], y_high[1], new_step);
      */     
   } /* end stepping */ ;

   if (fail == 0)        /* integration was successful */
   {

      printf(" Solution = \n");
      for (i = 0; i < n; ++i)
      {
         y_fin[i] = y_high[i];
         printf("%14.4e", y_fin[i]);
      }
      printf("\n");
   }
   else
      printf(" RKF45 failed to find a solution; fail=%d\n", fail);

} /* end main */


boolean f2(n, x, y, y_prime)
  /* Generate exponential function from first order eqn.

     y' = k  y

  */

   integer
      n;               /* number of dependent variables */
   real_8
      x,               /* independent variable */
      y[NEQN],         /* dependent variables */
      y_prime[NEQN];   /* derivatives to be calculated */
{
   if (k * x < 230.0)  /* within good range of x */
   {
      y_prime[0] = k * y[0];
      return(true);
   }
   else
      return (false);  /* say solution blows up */
} /* end f2 */

boolean f3(n, x, y, y_prime)
  /* Generate trig. function from second order eqn.

             2
     y'' = -k   y

   by transforming to 2 first order equations,

     y'  = y
      1     2

            2
     y'  = k  y
      2        1
  */

   integer
      n;               /* number of dependent variables */
   real_8
      x,               /* independent variable */
      y[NEQN],         /* dependent variables */
      y_prime[NEQN];   /* derivatives to be calculated */
{
   if (k * x < 230.0)  /* within good range of x for demo */
   {
      y_prime[0] = y[1];
      y_prime[1] = -k * k * y[0];
      return(true);
   }
   else
      return (false);  /* say solution blows up  for demo. */
} /* end f3 */