/*=====================================================================\
|     function solving  o.d.e. using  Runge-Kutta-Fehlberg method      |
|                          of order 7 and 8                            |
|----------------------------------------------------------------------|
|                                                                      |
|   double rk78 ( *at, x[], *ah, tol, hmin, hmax, n, deriv )           |
|                                                                      |
|     at          pointerr ato the variable TIME                       |
|     x[]         vector POSITION                                      |
|     ah          pointer to the variable STEP                         |
|     tol         TOLERANCE                                            |
|     hmin, hmax  minimal and maximal step                             |
|     n           dimension                                            |
|     deriv       pointer to the function containing the VECTOR FIELD  |
|                                                                      |
|                 void  ( * deriv ) ( t, x[], n, dx[] )                |
|                                                                      |
|     it returns an estimation of the ERROR done                       |
|                 or -1 if the dimension is greater than NMAX          |
|                                                                      |
\=====================================================================*/



static double alfa[ 13 ] =
   {            0.,      2. / 27.,       1. / 9.,       1. / 6.,
          5. / 12.,           .5 ,       5. / 6.,       1. / 6.,
           2. / 3.,       1. / 3.,            1.,            0.,
                1.};

static double beta[ 79 ] =
   {            0.,      2. / 27.,      1. / 36.,      1. / 12.,
          1. / 24.,            0.,       1. / 8.,       5. /12.,
                0.,   - 25. / 16.,     25. / 16.,        .5e-1 ,
                0.,            0.,          .25 ,           .2 ,
      - 25. / 108.,            0.,            0.,   125. / 108.,
       - 65. / 27.,    125. / 54.,    31. / 300.,            0.,
                0.,            0.,    61. / 225.,     - 2. / 9.,
        13. / 900.,            2.,            0.,            0.,
        - 53. / 6.,    704. / 45.,   - 107. / 9.,     67. / 90.,
                3.,  - 91. / 108.,            0.,            0.,
        23. / 108., - 976. / 135.,    311. / 54.,   - 19. / 60.,
          17. / 6.,    - 1. / 12., 2383. / 4100.,            0.,
                0., - 341. / 164., 4496. / 1025.,  - 301. / 82.,
     2133. / 4100.,     45. / 82.,    45. / 164.,     18. / 41.,
         3. / 205.,            0.,            0.,            0.,
                0.,    - 6. / 41.,   - 3. / 205.,    - 3. / 41.,
          3. / 41.,      6. / 41.,            0.,-1777. / 4100.,
                0.,            0., - 341. / 164., 4496. / 1025.,
      - 289. / 82.,  2193. /4100.,     51. / 82.,    33. / 164.,
         12. / 41.,            0.,            1.};

static double c7[ 11 ] =
   {    41. / 840.,            0.,            0.,            0.,
                0.,    34. / 105.,      9. / 35.,      9. / 35.,
         9. / 280.,     9. / 280.,    41. / 840.};

static double c8[ 13 ] =
   {            0.,            0.,            0.,            0.,
                0.,    34. / 105.,      9. / 35.,      9. / 35.,
         9. / 280.,     9. / 280.,            0.,    41. / 840.,
        41. / 840.};


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


#define NMAX 4000    
/*  NMAX is the maximal number of particles multiply by 4 (dim=2) */

#define max(a,b) ((a)<(b) ? (b) : (a))
#define sgn(a)   ((a)<0 ? -1 : 1 )
#define abs(a)   ((a)<0 ? -(a) : (a))


/* r[i][j] contains the distance between particles i and j */


double rk78 ( double *at, double x[], double *ah, double tol,
              double hmin, double hmax, int n,
              void ( *deriv ) ( double, double *, int, double * ))
   {
   double tpon, tol1, err, nor, kh, beth, h1, k [NMAX][13] ,
          x7 [NMAX], x8 [NMAX], xpon [NMAX], dx [NMAX];
   register int j, l;
   int i, m;

   if ( n > NMAX ) return -1;
   do {
      /*---------------> Computing K(i,j) <-----------------------*/
      m = 0;
      for ( i = 0; i < 13; ++i )
         {
         tpon = *at + alfa[i] * *ah;
         for ( j = 0; j < n; xpon[j] = x[j], ++j );
         for ( l = 0; l < i; ++l )
            {
            m=m+1;
            beth = *ah * beta[m];
            for ( j = 0; j < n; xpon[j] += beth * k[j][l], j=j+1);
            }
         ( *deriv )( tpon, xpon, n, dx );
         for ( j = 0; j < n; k[j][i] = dx[j], j=j+1 );
         }
      /*------> Computing the 2 points and associated values <--------------*/
      err = nor = 0;
      for ( j = 0; j < n; j++ )
         {
         x7[j] = x8[j] = x[j];
         for ( l = 0; l < 11; ++l )
            {
            kh = *ah * k[j][l];
            x7[j] += kh * c7[l];
            x8[j] += kh * c8[l];
            }
         x8[j] += *ah * ( c8[11] * k[j][11] + c8[12] * k[j][12] );
         err += abs ( x8[j] - x7[j] );
         nor += abs ( x8[j] );
         }
      err /= n;
      /*---------------> Computing the new step size h <-------------------*/
      tol1 = tol * ( 1 + nor / 100 );
      if ( err < tol1 ) err = max ( err, tol1 / 256 );
      h1 = *ah;
      *ah *= 0.9 * pow ( tol1 / err, 0.125 );
      if ( abs ( *ah ) < hmin ) *ah = hmin * sgn ( *ah );
      if ( abs ( *ah ) > hmax ) *ah = hmax * sgn ( *ah );
		} while (( err >= tol1 ) && ( abs ( *ah ) > hmin ));
   *at += h1;
   for ( j = 0; j < n; x[j] = x8[j], j=j+1 );
   return err;
   }