source: proiecte/NBody/NBI/NBIserial.c @ 162

#define CODENAME "NBI"
#define VERSION    "1"

/* dimension of physical space */
#define DIM                     3
/* max number of bodies, obviously */
#define MAXNBODIES              2100

/* TOLERANCES */

/* the extrapolation scheme stops when
        Tj,k - T(j-1),k becomes smaller than this */
#define GBSTOL                  5.0e-32

/* The solution of Kepler's equation is accepted when
   when iterates change by an amount smaller than this  */
#define KEPLERSOLVERTOLERANCE   5.0e-16

/* what method and order to use to initialize the CSN integrator, GBS or HWH  */
#define CSNINITGBS              1
#define CSNINITGBSORDER         4       
#define CSNINITHWH              2
#define CSNINITHWHORDER         4
int csninit = CSNINITGBS;


/* which extrapolation sequence to use in the GBS method */
#define EXSEQHARM int extrapolseq[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
#define EXSEQBUL  int extrapolseq[] = { 0, 1, 2, 3, 4, 6, 8, 12, 16, 24, 32};

/* cange this line accordingly */
EXSEQBUL

#define NOTTESTPARTICLE         1
#define TESTPARTICLE            0

/* the header files we need, only standard stuff */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/* Macro definitions */
#define NUMMETHODS      4
#define SKP 1
#define HWH 2
#define CSN 3
#define GBS 4

/* it is really the number of stages for GBS */
char *methodnames[] = { "   ", "SKP", "HWH", "CSN", "GBS" };
#define MAXCSN  15
#define MAXGBS  10
int MAXORDERS[] = { 0, 4, 4, MAXCSN, MAXGBS };

/* 
   MACROS WHICH DEFINE WHAT IS STORED FOR THE PARTICLES
        from 1 to 6: x y z vx vy vz
                  7: G*mass
        from 8 to 10:  accelerations in x, y and z
   e.g., particle i's y velocity is
         parts[i][5]
*/

#define TMPSIZE      (2*DIM+1) 
/* Note: we need only the G*mass for the actual integrations */ 
#define GMASSCO      (7)
#define ACC          (7)
/* you can add more data for particles here
   but keep NPARTDATA to specify the total number of
   coordinates and other components  */
#define NPARTDATA    (ACC+3)
/* this is used from index 1 in both indeces */
#define DECLARR [MAXNBODIES+1][NPARTDATA+1]

/* Cowell-Stormer-Newman array */
#define DECLSTORMER [MAXNBODIES+1][DIM+1][MAXCSN+2]
/* Gragg-Bulirsch-Stoer array */
#define DECLBS [MAXNBODIES+1][2*DIM+1][2][MAXGBS+2]

#define DIM2         (2*DIM)

/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                GLOBAL DATA
        ARRAYS THAT HOLD THE PARTICLES AND AUXILIARY DATA 
*/

/* ARRAYS THAT HOLD THE PARTICLES' DATA.
        This can be moved into main() without
        any adverse effect.  */
double parts DECLARR;

/* This holds test particle indices.
   It seems easier to keep this global. It is used, e.g.,
   by PlainAccel to avoid the computation of force between
   test particles.  */
int tpstore[MAXNBODIES+1];

/* this is needed for error messages from deep within the integrator */
FILE *FLOG;

/* Local error monitoring for the Cowell-Stormer integrator */
double CSNerror = 0;


/* FUNCTION AND STRUCT DECLARATIONS */

/* the tree to store Jacobi coordinates */
struct tree {
        double tmass; /* total G*mass */
        double rm1;   /* total mass / mass1 */
        double rm2;   /* total mass / mass2 */ 
        double bc[DIM2+1]; /* barycenter: position and velocity
                                for actual particle these  are pos. and vel. */
        double acc[DIM+1];
/* nonzero particle index indicates that this is one with finite mass
 or it is an array of test particles which can be processed in a simple loop,
   zero means that this is a node in the tree
*/
        int particleindex;  /* nonzero indicates actual particle(s) */
        struct tree *p1;
        struct tree *p2;
        int *testparticles; /* nonzero indicates array of test particles */
        };


int CStep(int ord, double h, int nump, double parts DECLARR, 
                struct tree *treein);

int GBSStep(int ord, double h, int nump, double parts DECLARR);

/* makes a simple chain, usual planetary motion case */
struct tree *MakeSimpleTree(int nump);
/* makes an arbitrary tree from a textual description, i.e.,
  parses the text and builds the tree */
struct tree *MakeTree(FILE *f);

/* Fill tree with barycenters, used to initialize the tree */
int FillCMinTree(struct tree *tree, double p DECLARR);
/* copy particle coordinates from tree;
   since the tree is the primary storage for the coordinates
   of particles with finite mass, one needs to call this
   before accessing data in the "parts" array
 */
int ParticlesFromTree(struct tree *tree, double p DECLARR);

/* The actual HWH integrator */
int HWHStep(int order, struct tree *tree, double stepsize, double p DECLARR);

/* Kinetic-potential energy split */
int SKPStep(int order, double stepsize, int nump, double p DECLARR);


/* access integer arrays in trees that hold the indeces
        of test particles attached to that point,
        returns pointer to testparticles,
        included for the sake of completeness, not actually used */
int **GetTestParticles(struct tree *tree, int ind);


/*    VERY SIMPLE IO, you can write your own if you need something more */
int readdata(FILE *f, int num, double p DECLARR) {
        int i, j; double tmp;
        for ( i = 1; i <= num; ++i ) {
                        fscanf(f, "%le", &tmp); p[i][GMASSCO] = tmp; }
        for ( i = 1; i <= num; ++i ) {
                for ( j = 1; j <= 6; ++j ) {
                        fscanf(f, "%le", &tmp); p[i][j] = tmp; } }
        return num; }

int printdataA(int num, FILE *f, double p DECLARR) {
        int i, j;
        for ( i = 1; i <= num; ++i ) {
                for ( j = 1; j <= DIM; ++j ) { fprintf(f, "%.14le  ", p[i][j]); } 
                fprintf(f, "%c", '\n');
                for ( j = DIM+1; j <= 2*DIM; ++j ) { fprintf(f, "%.14le  ", p[i][j]); } 
                fprintf(f, "%c", '\n');
                }
        return num; }

int printdata(int num, FILE *f, double p DECLARR, int c1, int c2) {
        int i, j;
        for ( i = 1; i <= num; ++i ) {
                for ( j = c1; j <= c2; ++j ) {
                        fprintf(f, "%.14le  ", p[i][j]); } 
                fprintf(f, "%c", '\n');
                }
        return num; }

/* total energy for particles with finite mass */
double energy(int nump, int dim, double p DECLARR, int *tps) {
        double h = 0, tmp, tmp2;
        int i,j,k;
        for ( i = 1; i <= nump; ++i ) {
                if ( tps[i] == 0 ) continue;
                tmp = 0.0;
                for ( k = 1; k <= dim; ++k ) {
                        tmp = tmp + p[i][k+dim]*p[i][k+dim]; }
                tmp = tmp*p[i][GMASSCO]/2.0;
                h = h + tmp;
                for ( j = i+1; j <= nump; ++j ) {
                        tmp = 0.0;
                        for ( k = 1; k <= dim; ++k ) {
                                tmp2 = p[i][k]-p[j][k];
                                tmp = tmp+tmp2*tmp2; }
                        tmp = 1/sqrt(tmp);
                        h = h - tmp*p[i][GMASSCO]*p[j][GMASSCO]; }
                }
        return h; }

/* z component of angular momentum */
double angmomz(int nump, int dim, double p DECLARR, int *tps) {
        double az = 0.0;
        int i;
        for ( i = 1; i <= nump; ++i ) {
                if ( tps[i] == 0 ) continue;
                az += (p[i][1]*p[i][2+dim]-p[i][2]*p[i][1+dim])*p[i][GMASSCO]; }
        return az; }

/* transforms all particles to coordinates relative to barycenter */
int tobarycenter(int nump, int dim, double p DECLARR) {
        double bc[10];
        int i, j;
        for (j = 0; j <= 2*dim; ++j ) bc[j] = 0.0;
        for ( i = nump; i >= 1; --i ) {
                bc[0] += p[i][GMASSCO];
                for (j = 1; j <= 2*dim; ++j ) bc[j] += p[i][j]*p[i][GMASSCO]; }
        for (j = 1; j <= 2*dim; ++j ) bc[j] = bc[j]/bc[0];
        for ( i = nump; i >= 1; --i ) {
                for (j = 1; j <= 2*dim; ++j ) p[i][j] -= bc[j]; }
        return 0; }




/* """"""""""""""""""""""""""""""""""""""""""""""""""""""""""
        MAIN 
*/

int main(int argc, char **argv) {
FILE *f, *fin, *flog, *flocerr;
char sicode[100];
int icode; /* SKP = 1 HWH = 2 CSN = 3 GBS = 4 */
struct tree *tree; int num, num2;
int ii, i, order = 1, nin, nps;
double stepsize;
double time = 0.0, ftime, ptime;
char fnamin[1000]; char id[1000];
double az, az0, toten, toten0;

/* get command-line argument */
if ( argc < 2 ) { exit(1); }
sscanf(argv[1], "%s", fnamin); 

/* read input file */
fin = fopen(fnamin, "r");
fscanf(fin, "%s", id); strcpy(fnamin, id);
/* create output and log file */
f = fopen(id, "w"); strcat(fnamin, ".log"); flog = fopen(fnamin, "w"); 
FLOG = flog;

#define FAILEDIN { fclose(fin); fclose(f); fclose(flog); exit(1); }

fprintf(flog, "ID = %s\n", id);
fscanf(fin, "%s", sicode); 
fscanf(fin, "%d", &order); fprintf(flog, "ORDER = %2d\n", order);
icode = 0;
for ( icode = 1; icode <= NUMMETHODS; ++icode ) {
        if ( strcmp(sicode, methodnames[icode] ) == 0 ) { 
                if ( order < 1 || order > MAXORDERS[icode] ) {
fprintf(flog, " Order %2d for method %s is not implemented \n", order, sicode);
        FAILEDIN }
                break; } }
if ( icode > NUMMETHODS ) {
fprintf(flog, " Invalid method code %s\n", sicode);
        FAILEDIN }


fscanf(fin, "%le", &stepsize); fprintf(flog, "STEPSIZE = %.12f\n", stepsize);
fscanf(fin, "%le", &ftime); fprintf(flog, "FINAL TIME = %.12f\n", ftime);
fscanf(fin, "%le", &ptime); fprintf(flog, "PRINTING TIME = %.12f\n", ptime);
fscanf(fin, "%d", &num); fprintf(flog, "NUMBER OF PARTS = %2d\n", num);
if ( num > MAXNBODIES ) {
fprintf(flog, " Too many particles, change MAXNBODIES macro in code ");
        FAILEDIN }


/* initialize tpstore */
for ( i = 1; i <= num+2; ++i ) { tpstore[i] = NOTTESTPARTICLE; }

/* this will TESTPARTICLE for test particles in tpstore */
tree = MakeTree(fin);

/* For CSN create error file */
if ( icode == CSN ) {
strcpy(fnamin, id); strcat(fnamin, ".lerr"); flocerr = fopen(fnamin, "w"); }

num2 = readdata(fin, num, parts);
/* DONE READING ALL INPUTS */
fclose(fin);

/* zero out test particles' tpstore, again, just in case */
for ( i = 1; i <= num; ++i ) { 
        if ( parts[i][GMASSCO] == 0.0 ) { tpstore[i] = TESTPARTICLE; } }

if ( num != num2 ) exit(1);
fprintf(flog, "GMs\n");
 printdata(num, flog, parts, GMASSCO, GMASSCO); 
fprintf(flog, "Inits \n");
         printdataA(num, flog, parts);
fflush(flog);

/* WATCH OUT FOR THIS 
   The transformation to barycenter
   tends to improve things. 
*/
tobarycenter(num, DIM, parts);
fprintf(flog, "Actual Barycentric Inits \n"); printdataA(num, flog, parts);

fprintf(flog, "Actual integrator is %s \n", sicode); 
fprintf(flog, "CODE: %s \n", CODENAME); 
fprintf(flog, "VERSION: %s \n", VERSION); 
fflush(flog);


/* z component of total angular momentum,
  it is used to monitor accumulation of roundoff */
az0 = angmomz(num, DIM, parts, tpstore);
toten0 = energy(num, DIM, parts, tpstore);



/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
        ACTUAL INTEGRATION
*/

nin = (int) (fabs(ptime)/fabs(stepsize));  /* the length of the inner loop, i.e., 
                                        between printing */
if ( nin == 0 ) {
        fprintf(flog, "%s\n", 
        " Stepsize is larger than printing time, aborted");
        exit(1); }
nps = (int) (fabs(ftime)/fabs(ptime));  /* the lenght of the outer loop, 
                        i.e., printings from beginning to end */

/* WARNING: VERY LONG INTEGRATIONS CAN PRODUCE VERY LARGE
        INTEGERS AND ONE HAS TO MAKE SURE THAT THEY ARE STILL
        SMALLER THAN THE LARGEST REPRESENTABLE INTEGER 
  ANOTHER WARNING: ROUNDOFF EFFECTS MAY CAUSE nin AND/OR nps TO BE
  DIFFERENT FROM WHAT THE USER EXPECTS
*/


/* print the first data point, i.e., the initial data */ 
#define PRINTDATA fprintf(f, "%.12f   ", time);\
        az = angmomz(num, DIM, parts, tpstore); \
        toten = energy(num, DIM, parts, tpstore); \
        fprintf(f, "%.16e  ", (toten-toten0)/toten0); \
        fprintf(f, "%.16e\n", (az-az0)/az0); \
         printdataA(num, f, parts);

#define PRINTERROR \
        if ( icode == CSN ) { \
        fprintf(flocerr, "%.12f    %.12e\n", time, CSNerror); }

PRINTDATA
PRINTERROR


/* MAIN LOOP */
for ( ii = 1; ii <= nps; ++ii ) {
for ( i = 1; i <= nin; ++i ) {
switch(icode) {
        case SKP:
          SKPStep(order, stepsize, num, parts); break;
        case HWH:
          HWHStep(order, tree, stepsize, parts); break;
        case CSN:
          CStep(order, stepsize, num, parts, tree); break;
        case GBS:
          GBSStep(order, stepsize, num, parts); break;
        default: 
                fprintf(flog, "Not implemented\n");
                fclose(flog); fclose(f); exit(1); break;
        } 
/* DONE ONE TIME STEP
        ADD MORE CODE HERE IF NECESSARY 
        variables you can use:
        double time;
        double parts DECLARR;
        Note that actual time should be computed as
        actualtime = time+stepsize*i
        */

}

 time += stepsize*nin;

 if ( icode == HWH ) { ParticlesFromTree(tree, parts); }

/* CHANGE THIS PART OR ADD YOUR CODE IF NECESSARY */
        
        PRINTDATA
        PRINTERROR

 if ( stepsize > 0.0 && time > ftime ) break;
 if ( stepsize < 0.0 && time < ftime ) break;

        }


fclose(f);
fprintf(flog, "\n Finished \n");
fclose(flog);
if ( icode == CSN ) fclose(flocerr);
        
return 0; }


/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        FOR TEST: THE SIMPLE SPLIT, I.E., Kinetic - Potential
        It does not provide special treatment for test particles.
*/
int PlainSplitKineticStep(int num, double stepsize, double p DECLARR) {
        int i, k;
        for ( i = 1; i <= num; ++i ) {
                for ( k = 1; k <= DIM; ++k ) {
                        p[i][k] += stepsize*p[i][k+DIM]; } }
        return 0; }

int PlainSplitPotentialStep(int num, double stepsize, double p DECLARR) {
        int i, j, k;
        double tmp, tmp2, tmp3;
        for ( i = 1; i <= num; ++i ) {
                tmp = 0.0;
                for ( j = i+1; j <= num; ++j ) {
                        tmp2 = 0.0;
                        for ( k = 1; k <= DIM; ++k ) {
                                tmp3 = p[i][k] - p[j][k];
                                tmp2 = tmp2 + tmp3*tmp3; }
                        tmp3 = 1.0/sqrt(tmp2);
                        tmp2 = tmp3/tmp2;
                        for ( k = 1; k <= DIM; ++k ) {
                                tmp3 = p[i][k] - p[j][k];
                                tmp = -tmp3*tmp2*p[j][GMASSCO]*stepsize;
                                p[i][k+DIM] += tmp;
                                tmp = tmp3*tmp2*p[i][GMASSCO]*stepsize;
                                p[j][k+DIM] += tmp; }
                        }
        }
        return 0; }



/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        COWELL-STORMER INTEGRATOR */


/* this computes accelerations, used both for CSN and GBS,
        not to be confused with Accel used for HWH */

int PlainAccel(int num, double p DECLARR) {
        int i, j, k; int tpi, tpj;
        double tmp, tmp2, tmp3;
        double massi, massj;
        for ( i = 1; i <= num; ++i ) {
                for ( k = 1; k <= DIM; ++k ) { p[i][ACC+k] = 0.0; } }
/* for testing things
for ( i = 1; i <= num; ++i ) {
        for ( k = 1; k <= DIM; ++k ) { p[i][ACC+k] = 1.0; } }
return 0;
*/
        for ( i = num; i >= 1; --i ) { 
                tpi = tpstore[i];
                massi = p[i][GMASSCO];
                for ( j = i+1; j <= num; ++j ) {
                        tpj = tpstore[j];
                        if ( tpi == TESTPARTICLE && tpj == TESTPARTICLE ) continue;
                        tmp2 = 0.0;
                        for ( k = 1; k <= DIM; ++k ) {
                                tmp3 = p[i][k] - p[j][k];
                                tmp2 = tmp2 + tmp3*tmp3; }
                        tmp3 = 1.0/sqrt(tmp2);
                        tmp2 = tmp3/tmp2;
                        massj = p[j][GMASSCO];
                        for ( k = 1; k <= DIM; ++k ) {
                                tmp3 = p[i][k] - p[j][k];
                                tmp = tmp3*tmp2;
                        if ( tpj == NOTTESTPARTICLE ) {
                                p[i][ACC+k] -= (tmp*massj); }
                        if ( tpi == NOTTESTPARTICLE ) {
                                p[j][ACC+k] += (tmp*massi); }
                                }
                        }
        }
        return 0; }

/* Cowell-Stormer-Newman 14(6)th order */


/* Newman's 14(6)th order Stormer's coeffs */

double CSalpha[] =  {   0.0,
                        1.000000000000000000000,
                        0.000000000000000000000,
                        0.083333333333333333333,
                        0.083333333333333333333,
                        0.079166666666666666666,
                        0.075000000000000000000,
                        0.071345899470899470899,
                        0.068204365079365079365,
                        0.065495756172839506172,
                        0.063140432098765432098,
                        0.061072649861712361712,
                        0.059240564123376623376,
                        0.057603625837453135733,
                        0.056129980884507174190,
                        0.054794379107071147139
                                                 };
 
 
double CSbeta[] =  {    0.0,
                        1.50000000000000000000,
                        0.33333333333333333333,
                        0.37500000000000000000,
                        0.35277777777777777777,
                        0.33402777777777777777,
                        0.31924603174603174603,
                        0.30736607142857142857,
                        0.29757660934744268077,
                        0.28933077050264550264,
                        0.28225737868098979210,
                        0.27609762576993479771,
                        0.27066578505957226195,
                        0.26582499331955247133,
                        0.26147199790503706399,
                        0.25752728193649391773
                                               };
 
        

double  sumfrombottom(double *a, int nr, double *coef) {
        int i;
        double sum = 0.0;
        for ( i = nr; i >= 1;  --i ) {
                sum += a[i]*coef[i]; }
        return sum; }

int differencing(double *a, int nr, double in) {
        int i, j, p;
        double tmp1, tmp2;
        tmp1 = a[1]; a[1] = in;
        for ( i = 2; i <= nr; ++i ) {
                tmp2 = a[i]; a[i] = a[i-1]-tmp1; tmp1 = tmp2; } 
        return 0; }



static int CStepCount = 0;
static int CStepInit = 0;
static struct tree *initCStree;

int CStep(int ord, double h, int nump, 
                double parts DECLARR, struct tree *treein) {
        double *a, asav, vx, sum, tmp; int p,j, i;
static double df DECLSTORMER;
/* actual local order to determine how terms to sum */
if ( CStepInit == 0 ) {
        initCStree = treein;
        CStepInit = 1;
        for ( i = 1; i <= ord+1; ++i ) {
        for ( p = 1; p <= nump; ++p ) {
        for ( j = 1; j <= DIM; ++j ) {
                df[p][j][i] = 0.0; } } }        }
if ( CStepInit == 1 ) {
/* initialization stage */
                for ( p = 1; p <= nump; ++p ) {
                for ( j = 1; j <= DIM; ++j ) {
                        df[p][j][0] = -parts[p][j]; } }

                if ( csninit == CSNINITHWH ) {
                        HWHStep(CSNINITHWHORDER, initCStree, h, parts); 
                        ParticlesFromTree(initCStree, parts);  }
                else { GBSStep(CSNINITGBSORDER, h, nump, parts); }
                if ( CStepCount > ord ) CStepInit = 2;
                for ( p = 1; p <= nump; ++p ) {
                for ( j = 1; j <= DIM; ++j ) {
                        df[p][j][0] += parts[p][j]; } }
                }
else {
/* actual run */
CSNerror = 0;
for ( p = 1; p <= nump; ++p ) {
for ( j = 1; j <= DIM; ++j ) {
/* formula: x(n+1) - x(n) = (x(n) - x(n-1)) + h^2*sum
        with X(n) = x(n)-x(n-1)
        X(n+1) = X(n) + h^2*sum
        x(n+1) = X(N+1) + x(n)
X(n) is in a[0] 
*/
        a = df[p][j];
        asav = a[0];
        sum = sumfrombottom(a, ord, CSalpha);
        sum *= h*h;
        a[0] += sum;
        parts[p][j] += a[0];
/* vx = (1.0/h)*(x(n)-x(n-1)) + h*sum */
        sum = h*sumfrombottom(a, ord, CSbeta);
        vx = asav/h + sum;
        parts[p][DIM+j] = vx;
/* record error */
        tmp = fabs(CSalpha[ord]*a[ord]*h*h);
        if ( CSNerror < tmp ) CSNerror = tmp;
        } } }
PlainAccel(nump, parts);
for ( p = 1; p <= nump; ++p ) {
for ( j = 1; j <= DIM; ++j ) {
        differencing(df[p][j], ord, parts[p][ACC+j]);
        parts[p][ACC+j] = 0.0;
        } }
++CStepCount;
return 0; }



/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        BULIRSCH-STOER
        actually, GRAGG-BULIRSCH-STOER, following Hairer et al., 1993
*/

/* Gragg's symmetric integrator, this is used to generate Tj,1  */
int GStep(int nump, double parts DECLARR, double s DECLARR, double h, int len) {
int k, j, p, d, knext, kprev, m;
double dif, difsum, tmp1, astep;
static double y0 DECLARR;
static double y1 DECLARR;
/* save originial in y0 */
for ( p = 1; p <= nump; ++p ) {
                y0[p][GMASSCO] = parts[p][GMASSCO];
for ( d = 1; d <= DIM; ++d ) { 
                y0[p][d] = parts[p][d]; 
                y0[p][d+DIM] = parts[p][d+DIM]; } }
/* consult Hairer et al. why we need this
        in short: we need a method which has only even powers
        of h in the expansion of its error, i.e., a symmetric method
 */
h = h/2.0;
PlainAccel(nump, y0);
for ( p = 1; p <= nump; ++p ) {
/* load masses into y1 */
                y1[p][GMASSCO] = y0[p][GMASSCO];
for ( d = 1; d <= DIM; ++d ) { 
                y1[p][d] = y0[p][d]+h*y0[p][DIM+d]; 
                y1[p][d+DIM] = y0[p][d+DIM]+h*y0[p][ACC+d]; } }
/* compute y2 (m=1), y3' (m=2), y4 (m=3), y5' (m=4) y(2n+1)' (m=2*n) 
        n = len
        */
for ( m = 1; m <= 2*len; ++m ) {
        if ( (m%2) ) {
        for ( p = 1; p <= nump; ++p ) {
        for ( d = 1; d <= DIM; ++d ) { 
                y1[p][d] = y0[p][d]+2.0*h*y1[p][d+DIM]; 
                s[p][d] = y1[p][d]; y0[p][d] = y1[p][d]; 
                }  } }
        else {
        PlainAccel(nump, y1);
        for ( p = 1; p <= nump; ++p ) {
        for ( d = 1; d <= DIM; ++d ) { 
                y0[p][d+DIM] = y1[p][d+DIM]; 
                y1[p][d+DIM] = y0[p][d+DIM]+2.0*h*y1[p][ACC+d]; 
                s[p][d+DIM] = (y0[p][d+DIM]+y1[p][d+DIM])/2.0; 
                } } 
                }
        }
return 0; }


int GBSStep(int bsord, double h, int nump, double parts DECLARR) {
int k, j, jfinal, p, d, jnext, jprev, len;
static double gbsarr DECLBS;
double dif, difsum, tmp1, astep;
int *n = extrapolseq;
static double s DECLARR;
/* for testing GStep
        astep = h; len = 1;
        GStep(nump, parts, s, astep, len);
for ( p = 1; p <= nump; ++p ) {
for ( d = 1; d <= 2*DIM; ++d ) { parts[p][d] = s[p][d]; } } 
return 0;
*/
jfinal = bsord;
for ( j = 1; j <= bsord; ++j ) {
        difsum = 0.0;
/* map j and j-1 to storage index */
        jnext = ((j)%2); 
        jprev = ((j-1)%2);
/* compute Tj,1 */
        len = n[j];
        astep = h/len;
        GStep(nump, parts, s, astep, len);
        for ( p = 1; p <= nump; ++p ) {
        for ( d = 1; d <= 2*DIM; ++d ) { gbsarr[p][d][jnext][1] = s[p][d]; } }
/* compute Tj,k */
for ( k = 1; k <= j-1; ++k ) {
/* Aitken-Neville for nonsymmetric underlying method */
        tmp1 = ((double) n[j])/((double) n[j-k]); 
/* Aitken-Neville for symmetric underlying method (square of previous) */
        tmp1 = tmp1*tmp1;
for ( p = 1; p <= nump; ++p ) {
for ( d = 1; d <= 2*DIM; ++d ) {
        dif = (gbsarr[p][d][jnext][k] - gbsarr[p][d][jprev][k]);
        if ( k == j-1 ) difsum += fabs(dif);
        gbsarr[p][d][jnext][k+1] = gbsarr[p][d][jnext][k] + dif/(tmp1-1.0);
        } }
        } 
/* stop computing more is reached roundoff level */
        if ( j > 1 && difsum <= GBSTOL ) { jfinal = j; break; }
        }
j = jfinal; k = j; jnext = ((j)%2); 
for ( p = 1; p <= nump; ++p ) {
for ( d = 1; d <= 2*DIM; ++d ) { parts[p][d] = gbsarr[p][d][jnext][k]; } } 
return 0; }
        





/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        Actual  HWH integration */
int MoveCMandKepler(struct tree *tree, double *cmove, double stepsize, 
                double p DECLARR);
int KeplerOnTree(struct tree *tree, double stepsize, double p DECLARR);
int RecursivePerturbations(struct tree *tree, double stepsize, double p DECLARR);
int RecPertAcc(struct tree *tree, double stepsize, double p DECLARR);

/* fourth order coeffs */
static double Coeff4x0 = 0.0;
static double Coeff4x1 = 0.0;

static int SympCalled = 0;

void initSympCoefs() {
        double cbrt2 = pow(2.0, 1.0/3); //cbrt(2.0);
        Coeff4x0 = -cbrt2/(2.0-cbrt2);
        Coeff4x1 = 1.0/(2.0-cbrt2);
        SympCalled = 1;
}
        

static int HWHCalled = 0;

int HWHStep(int order, struct tree *tree, double stepsize, double p DECLARR) {
        if ( SympCalled == 0 ) { initSympCoefs(); }
        if ( HWHCalled == 0) {
/* initialize tree */
                FillCMinTree(tree, p); 
                HWHCalled = 1; }
        if ( order == 1 ) {
                KeplerOnTree(tree, stepsize, p);
                RecursivePerturbations(tree, stepsize, p);
                return 0; }
        if ( order == 2 ) {
                KeplerOnTree(tree, stepsize/2.0, p);
                RecursivePerturbations(tree, stepsize, p);
                KeplerOnTree(tree, stepsize/2.0, p); 
                return 0; }
        if ( order == 4 ) {
                HWHStep(2, tree, stepsize*Coeff4x1, p);
                HWHStep(2, tree, stepsize*Coeff4x0, p);
                HWHStep(2, tree, stepsize*Coeff4x1, p);
                return 0; }
        printf(" Order %2d  is not implemented, Aborted", order);
        exit(1); }


int SKPStep(int order, double stepsize, int nump, double p DECLARR) {
        if ( SympCalled == 0 ) { initSympCoefs(); }
        if ( order == 1 ) {
                PlainSplitKineticStep(nump, stepsize, p);
                PlainSplitPotentialStep(nump, stepsize, p);
                return 0; }
        if ( order == 2 ) {
                PlainSplitKineticStep(nump, stepsize/2.0, p);
                PlainSplitPotentialStep(nump, stepsize, p);
                PlainSplitKineticStep(nump, stepsize/2.0, p);
                return 0; }
        if ( order == 4 ) {
                SKPStep(2, stepsize*Coeff4x1, nump, p);
                SKPStep(2, stepsize*Coeff4x0, nump, p);
                SKPStep(2, stepsize*Coeff4x1, nump, p);
                return 0; }
        printf(" Order %2d  is not implemented, Aborted", order);
        exit(1); }



/* +++++++++++++++++++++++++++++++++++++++++++++++++++++ 
/       TREE
/
/
*/


int inittree(struct tree *tmp) {
        int i; tmp->tmass = 0.0; tmp->particleindex = 0;
        tmp->p1 = 0; tmp->p2 = 0;
        tmp->testparticles = 0; tmp->rm1 = 0; tmp->rm2 = 0;
        for ( i = 1; i <= DIM2; ++i ) { tmp->bc[i] = 0; }
        for ( i = 1; i <= DIM; ++i ) { tmp->acc[i] = 0; }
        return 0; }


struct tree *newtree() {
        struct tree *tmp = (struct tree*) malloc(sizeof(struct tree)); 
        inittree(tmp); return tmp; }

void deletetree(struct tree *p) { 
        if ( p->p1 ) deletetree(p->p1);
        if ( p->p2 ) deletetree(p->p2);
        free((void *) p); }


struct tree *MakeSimpleTree(int nump) {
        int i;
        struct tree *root;
        struct tree *t1, *p1, *p2;
        root = newtree();
        t1 = root;
        for ( i = nump; i > 1; --i ) {
                p1 = newtree();
                p2 = newtree();
                p1->particleindex = 0;
                p2->particleindex = i;
                t1->p1 = p1; t1->p2 = p2;
                t1 = p1; }
        t1->particleindex = 1;
        return root; }

int eatwhite(FILE *f) {
        int c; 
        for ( ; ; ) {
                c = getc(f);
                if ( !isspace(c) ) break; }
        ungetc(c, f); 
        return c; }


/* just an array to hold indeces read in,
        portions of this are assigned to
        tpa's below */
int tmpReadList[MAXNBODIES+1000];
/* bookkeeping */
int freelist = 0;

int *ReadList(FILE *f) {
        int i1, i2; int c; int i; int *arr;
        c = eatwhite(f);
        if ( c != '[' ) {
                printf(" Error: bad list, aborted \n"); exit(1); }
        c = getc(f);
        fscanf(f, "%d - %d", &i1, &i2);
        c = eatwhite(f);
        if ( c != ']' ) {
                printf(" Error: bad list, aborted \n"); exit(1); }
        c = getc(f); 
        arr = tmpReadList+freelist;
/* put zero in tpstore */
        for ( i = i1; i <= i2; ++i ) tpstore[i] = TESTPARTICLE;
        for ( i = i1; i <= i2; ++i ) arr[i-i1+1] = i;
        arr[0] = i2-i1+1;
        freelist += i2-i1+2;
        return arr; }
        
                

struct tree *MakeTree(FILE *f) {
        int i; int c; int pind; int *tpa;
        struct tree *root;
        struct tree *p1, *p2;
        c = eatwhite(f);
        if ( c != '(' ) {
                if ( isdigit(c) ) {
/* single particle */
                        fscanf(f,"%d",&pind);
                        root = newtree();
                        root->particleindex = pind;
                        return root; } 
/* set of test particles */
                tpa = ReadList(f);
                root = newtree();
                root->particleindex = tpa[1];
                root->testparticles = tpa;
                return root; }
        c = getc(f);
        root = newtree();
        p1 = MakeTree(f);
        p2 = MakeTree(f);
        root->p1 = p1; root->p2 = p2;
        root->particleindex = 0;
        c = eatwhite(f);
        if ( c != ')' ) {
                printf(" Error: bad tree, aborted \n"); exit(1); }
        c = getc(f);
        return root; }



/* not used now but might come handy */
struct tree *getparticleintree(struct tree *tree, int ind) {
        struct tree *tmp;
        if ( tree->particleindex ) { 
                if ( tree->particleindex == ind ) { return tree; }
                return 0; }
        tmp = getparticleintree(tree->p1, ind);
        if ( tmp != 0 ) return tmp;
        tmp = getparticleintree(tree->p2, ind);
        if ( tmp != 0 ) return tmp;
        return 0; }
        

/* not used now but might come handy */
int **GetTestParticles(struct tree *tree, int ind) {
        struct tree *tmp = getparticleintree(tree, ind);
        if ( tmp == 0 ) return 0;
        return &(tmp->testparticles); }
        

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ACTUAL NUMERICAL STUFF ON TREE
*/

/* fills in initial data: masses, mass ratios, barycenters */
int FillCMinTree(struct tree *tree, double p DECLARR) {
        int j; double *bc = tree->bc;
        int i, ind, tp, num, *tpa = tree->testparticles;
        for ( j = 1; j <= DIM; ++j ) tree->acc[j] = 0.0;
        tree->tmass = 0.0; tree->rm1 = 0.0; tree->rm2 = 0.0;
        for ( j = 1; j <= DIM2; ++j ) bc[j] = 0.0;
        ind = tree->particleindex;
           if ( ind ) {
              if ( tpa != 0 ) {  num = tpa[0];
                for ( tp = 1; tp <= num; ++tp ) { 
                        ind = tpa[tp]; if ( ind < 1 ) continue;
                        for ( i = 1; i <= DIM; ++i ) p[ind][ACC+i] = 0.0; } }  
               else { tree->tmass = p[ind][GMASSCO];
                for ( j = 1; j <= 2*DIM; ++j ) bc[j] = p[ind][j]; }
           return 0; } 
        FillCMinTree(tree->p1, p);
        FillCMinTree(tree->p2, p);
        tree->tmass = tree->p1->tmass + tree->p2->tmass;
        tree->rm1 = (tree->p1->tmass)/(tree->tmass); 
        tree->rm2 = (tree->p2->tmass)/(tree->tmass);
        for ( j = 1; j <= DIM2; ++j ) {
                bc[j] = tree->p1->bc[j] * tree->rm1
                + tree->p2->bc[j] * tree->rm2; } 
        return 0; }

/* copy from tree into p */
int ParticlesFromTree(struct tree *tree, double p DECLARR) {
        int j; double *bc = tree->bc;
        int i, ind, *tpa = tree->testparticles;
        ind = tree->particleindex;
        if ( ind ) {
           if ( tpa == 0 ) {
                tree->tmass = p[ind][GMASSCO];
                for ( j = 1; j <= DIM2; ++j ) p[ind][j] = bc[j]; }
           return 0; } 
        ParticlesFromTree(tree->p1, p);
        ParticlesFromTree(tree->p2, p);
        return 0; }


/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        KEPLER STEP ON TREE
*/

/* computes Stumpff functions several ways, depending on x */
int StumpffFunctions(double *vals, double x);
/* advances the Kepler motion 1/2*v^2 - mu/|x| for time deltat */
double KeplerStumpffStep(int dim, double mu, double *xv, double deltat);
/* return only the increment in xv */
double IncrementKeplerStumpffStep(int dim, double mu, double *xv, double deltat);

/* NOTE: this does not CHANGE centers of masses */
int KeplerOnTree(struct tree *tree, double stepsize, double p DECLARR) { 
/* move the center of mass */
        double cmove[10]; int i;
        for ( i = 1; i <= 2*DIM; ++i ) cmove[i] = 0.0;
/* move by 1/2*(Tmass)*v^2 */
        for ( i = 1; i <= DIM; ++i ) cmove[i] = (tree->bc[DIM+i])*stepsize; 
        MoveCMandKepler(tree, cmove, stepsize, p);
return 0; }


/* Assume that test particles are attached to 2nd */
int DoTestPartKepler(struct tree *tree, double *cmove, 
                        double stepsize, double p DECLARR) {
        int i, tp, num, *tpa, ind;
        double mu, tmpkep[TMPSIZE], *bc;
        tpa = tree->p2->testparticles; 
        num = tpa[0];
        if ( num > 0 ) {
                bc = tree->p1->bc; mu = tree->tmass;
        for ( tp = 1; tp <= num; ++tp ) { ind = tpa[tp]; if ( ind < 1 ) continue;
                for ( i = 1; i <= DIM2; ++i )   tmpkep[i] = p[ind][i] - bc[i];
                IncrementKeplerStumpffStep(DIM, mu, tmpkep, stepsize); 
                for ( i = 1; i <= DIM2; ++i ) p[ind][i] += (cmove[i]+tmpkep[i]);
                } }
        return 0; }

/* move the center of mass and do a Kepler step on relative
  coordinates, transmit the new "moves" down */
int MoveCMandKepler(struct tree *tree, double *cmove, double stepsize, 
                double p DECLARR) {
        int i, *tpa, ind;
        double tmp[TMPSIZE], tmpkep[TMPSIZE], *xv1, *xv2, mu;
        double *bc = tree->bc;
        for ( i = 1; i <= DIM2; ++i ) { bc[i] += cmove[i]; }
/* do nothing else for actual particles */
        if ( tree->particleindex ) { return 0; }
        if ( tree->p2->testparticles ) {        
/* this is a pair where p2 is an array of test particles */
                DoTestPartKepler(tree, cmove, stepsize, p); 
                MoveCMandKepler(tree->p1, cmove, stepsize, p); 
                return 0; }
        mu = tree->tmass;
        xv1 = tree->p1->bc; xv2 = tree->p2->bc;
        for ( i = 1; i <= DIM2; ++i ) { tmpkep[i] = xv2[i]-xv1[i]; }
/* do Kepler step on tmpkep, we get back the increment */
        IncrementKeplerStumpffStep(DIM, mu, tmpkep, stepsize); 
/* use only the relative change */
        for ( i = 1; i <= DIM2; ++i ) { tmp[i] = cmove[i] - tmpkep[i]*(tree->rm2); }
        MoveCMandKepler(tree->p1, tmp, stepsize, p);
        for ( i = 1; i <= DIM2; ++i ) { tmp[i] = cmove[i] + tmpkep[i]*(tree->rm1); }
        MoveCMandKepler(tree->p2, tmp, stepsize, p);
        return 0; }
        
/*  ++++++++++++++++++++++++++++++++++++++++++++++++++++++
        PERTURBATION ACCELERATIONS
*/

/* simple acceleration */
int Accel(double *x1, double *x2, double *acc) {
        double tmp[10], r = 0.0; int i; 
        for ( i = 1; i <= DIM; ++i ) {
                tmp[i] = x1[i]-x2[i];
                r = r + tmp[i]*tmp[i]; }
        r = sqrt(r);
        r = r*r*r;
        for ( i = 1; i <= DIM; ++i ) { acc[i] = -tmp[i]/r; }
        return 0; }

int AdvanceTestParticleVel(struct tree *tree, double stepsize, double p DECLARR) {
                int i, ind, tp, num, *tpa = tree->testparticles;
                if ( tpa == 0 ) return 0; num = tpa[0];
        for ( tp = 1; tp <= num; ++tp ) {
                ind = tpa[tp]; if ( ind < 1 ) continue;
                for ( i = 1; i <= DIM; ++i ) {
                        p[ind][DIM+i] += stepsize*p[ind][ACC+i]; 
                        p[ind][ACC+i] = 0.0; }
                        } 
        return 0; }
                        
        
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        PERTURBATIONS ON TREE
*/

int RecursivePerturbations(struct tree *tree, double stepsize, double p DECLARR) {
        int j, ind; double *bc = tree->bc;
/* advance attached test particles */
        if ( tree->particleindex ) {
/* if particle we are ready to advance v */
                if ( tree->testparticles ) {
                        AdvanceTestParticleVel(tree, stepsize, p); 
                        return 0; }
/* this is a nonzero mass actual particle, advance v */
                for ( j = 1; j <= DIM; ++j ) {
                        bc[DIM+j] += stepsize*tree->acc[j]; tree->acc[j] = 0.0; }
                return 0; }
/* if not, first we take care of accs. for pairs with the two parts  
        IMPORTANT TO DO THIS FIRST
*/
        RecPertAcc(tree, stepsize, p);
/* now descend and do the same for both parts */
RecursivePerturbations(tree->p2, stepsize, p);
RecursivePerturbations(tree->p1, stepsize, p);
/* now regenerate center of mass etc. */
if ( tree->p2->testparticles ) {
        for ( j = 1; j <= DIM2; ++j ) { bc[j] = tree->p1->bc[j]; }
        return 0; }
for ( j = 1; j <= DIM2; ++j ) {
        bc[j] = tree->p1->bc[j] * tree->rm1
        + tree->p2->bc[j] * tree->rm2; }
return 0; }

/* assume test particles in tree2 for dir > 1, reverse signs for dir < 0 */
int TestParticleAcc(int dir, double *cacc, struct tree *tree1, struct tree *tree2, 
                                double p DECLARR) {
                int i, ind, tp, num, *tpa = tree2->testparticles;
                double acc[TMPSIZE];
                double *bc = tree1->bc;
                double tmp, gm1 = tree1->tmass;
                num = tpa[0];
        for ( tp = 1; tp <= num; ++tp ) {
                ind = tpa[tp]; if ( ind < 1 ) continue;
/* rewrite this to avoid bad numerical attitude */
                        Accel(bc, &p[ind][0], acc);
                for ( i = 1; i <= DIM; ++i ) {
/* since test particles are in 2nd */
                        tmp = acc[i]-cacc[i]*dir;
                        p[ind][ACC+i] += -gm1*tmp; }
                        } 
        return 0; }
                        

/* c1 and c2 are for possible future use */
int RecPertAccSub(struct tree *tree1, 
        struct tree *tree2, double *cacc, double *c1, double *c2, double p DECLARR) {
/* do the work if both are particles */
        int i, ind1, ind2; double acc[TMPSIZE];
        double tmp, gm1 = tree1->tmass, gm2 = tree2->tmass;
        ind1 = tree1->particleindex;
        ind2 = tree2->particleindex;
        if ( ind1 && ind2 ) {
/* do nothing if both are test particle arrays */
                if ( tree1->testparticles  && tree2->testparticles  ) return 0;  
/* test particles first */
                if ( tree2->testparticles ) {
                        TestParticleAcc(1, cacc, tree1, tree2, p); return 0; }
                if ( tree1->testparticles ) {
                        TestParticleAcc( -1, cacc, tree2, tree1, p); return 0; }
/* if actual particles with nonzero mass */
                Accel(tree1->bc, tree2->bc, acc);
                for ( i = 1; i <= DIM; ++i ) {
/* rewrite this to avoid bad numerical attitude */
                        tmp = acc[i]-cacc[i];
                        tree1->acc[i] += gm2*tmp; tree2->acc[i] += -gm1*tmp; }
        return 0; }
/* if not, descend but only for one branch, the other is taken
                care in turn */
        if ( !ind1 ) {
                RecPertAccSub(tree1->p2, tree2, cacc, c1, c2, p);
                RecPertAccSub(tree1->p1, tree2, cacc, c1, c2, p);
                return 0; }
        RecPertAccSub(tree1, tree2->p1, cacc, c1, c2, p);
        RecPertAccSub(tree1, tree2->p2, cacc, c1, c2, p);
        return 0; }

/* assume test particles in tree2 */
int TestParticleAccSpec(double *cm, struct tree *tree1, 
                        struct tree *tree2, double p DECLARR) {
                int i, ind, tp, num, *tpa = tree2->testparticles;
                double acc[TMPSIZE], cacc[TMPSIZE];
                double *bc = tree1->bc;
                double tmp, gm1 = tree1->tmass;
                num = tpa[0];
        for ( tp = 1; tp <= num; ++tp ) {
                ind = tpa[tp]; if ( ind < 1 ) continue;
/* rewrite this to avoid bad numerical attitude */
                        Accel(bc, &p[ind][0], acc);
                        Accel(cm, &p[ind][0], cacc);
                for ( i = 1; i <= DIM; ++i ) {
/* since test particles are in 2nd */
                        tmp = acc[i]-cacc[i];
                        p[ind][ACC+i] += -gm1*tmp; 
                                }
                        } 
        return 0; }
                        

/* c1 and c2 are for possible future use */
int RecPertAccSubTest(double *cm,  struct tree *tree1, struct tree *tree2, double p DECLARR) {
        int ind1;
        ind1 = tree1->particleindex;
        if ( ind1 ) {
/* do nothing if tree1 has only test particles */
                if ( tree1->testparticles ) return 0;
                TestParticleAccSpec(cm, tree1, tree2, p); 
                return 0; }
/* if not, descend for both branches */
        RecPertAccSubTest(cm, tree1->p1, tree2, p);
        RecPertAccSubTest(cm, tree1->p2, tree2, p);
        return 0; }


int RecPertAcc(struct tree *tree, double stepsize, double p DECLARR) {
        double cacc[TMPSIZE];
/* do nothing if all in p2 are attached to p1: pure Kepler motion */
        if ( tree->p1->particleindex &&  tree->p2->particleindex  )  return 0; 
        if ( tree->p2->testparticles == 0 ) {
/* actual nonzero masses in p2 */
        Accel(tree->p1->bc, tree->p2->bc, cacc);
        RecPertAccSub(tree->p1, tree->p2, cacc, tree->p1->bc, tree->p2->bc, p); 
        return 0;}
/* test particles in p2 */
        RecPertAccSubTest(tree->p1->bc, tree->p1, tree->p2, p); 
        return 0; }




/*
// ===============================================================================
//
//   KEPLER MOTION USING STUMPFF FUCNTIONS
//      FOLLOWING DANBY, 1992
//      (quartic Newtonian solver, Stumpff function definitions,
//       f and g functions)
//      see also: Stiefel and Sheifele, Linear and Regular
//      Celestial Mechanics
//
// ===============================================================================

// Stumpff functions (c0, c1, c2, c3)
*/

/* CHANGE these numbers to tune accuracy */
/* this one can be changed runtime */
static double StumpffQSTol = KEPLERSOLVERTOLERANCE;

/* the limit for which trig and hyp. trig functions are used */
#define STUMPFFPLAINLIM (5.0)
/* reduction limit */
#define STUMPFFRED (1.0/32.0)
#define STUMPFFSHORTTOL (1.0e-6)

/* this one is not used now */
static double StumpffSerTol = 1.0e-25;


static int StumpffErrorCount = 0;
static int StumpffErrorCountLimit = 100;

CheckStumpffErrorCount() {
        ++StumpffErrorCount;
        if ( StumpffErrorCount > StumpffErrorCountLimit ) {
                fprintf(FLOG, "\n Too many errors, aborted \n");
                fflush(FLOG);
                exit(1); }
        return 0; }

double util_fact(int k) {
        int i; double fact = 1;
        for ( i = 1; i <= k; ++i ) {
                fact *= i; }
        return fact; }

/*
In case we need the coefficients, produced by Mathematica 
static double c2_0 = 0.5;
static double c2_1 = -1.0/24;
static double c2_2 = 1.0/720;
static double c2_3 = -1.0/40320;
static double c2_4 = 1.0/3628800;
static double c2_5 = -1.0/479001600;
static double c2_6 = 1.0/87178291200;
static double c2_7 = -1.0/20922789888000;
static double c2_8 = 1.0/6402373705728000;
static double c2_9 = -1.0/2432902008176640000;                          
static double c2_10 = 1.0/1124000727777607680000;

static double c3_0 = 1.0/6;
static double c3_1 = -1.0/120;
static double c3_2 = 1.0/5040;
static double c3_3 = -1.0/362880;
static double c3_4 = 1.0/39916800;
static double c3_5 = -1.0/6227020800;
static double c3_6 = 1.0/1307674368000;
static double c3_7 = -1.0/355687428096000;
static double c3_8 = 1.0/121645100408832000;
static double c3_9 = -1.0/51090942171709440000;
static double c3_10 = 1.0/25852016738884976640000;
*/


/* This is what is actually used, generated by factoring the coefficients */
double HornerStumpff2(double x) {
        return
(1.0 - x*(1 - x*(1 - x*(1 - x*(1 - x *(1 - x
/* up to degree 10, if you need it 
*(1 - x*(1 - x*(1 - (1 - x/462)*x/380)/306)/240)        */
/182)/132)/90)/56)/30)/12)/2;
        }

double HornerStumpff3(double x) {
        return
 (1.0 - x*(1 - x*(1 - x*(1 - x*(1 - x*(1 - x
/* up to degree 10, if you need it 
*(1 - x*(1 - x*(1 - (1 - x/506)*x/420)/342)/272)        */
/210)/156)/110)/72)/42)/20)/6;
        }


/* first 4 terms, ie., up to x^3 */
double Stumpff2Short(double x) {
        return (1.0 - x*(1 - x*(1 - x/56)/30)/12)/2;
        }

double Stumpff3Short(double x) {
        return (1.0 - x*(1 - x*(1 - x/72)/42)/20)/6;
        }




/* Original for reference */
double StumpffSeriesCk(int k, double x) {
        int i;
        double term, sum = 0.0;
        double terms[100], num;
        term = 1.0;
        for ( i = 1; i <= 50; ++i ) {
                num = (k+2*i-1)*(k+2*i);
                term = term*(-x/num);
                terms[i] = term;
                if ( fabs(term) < StumpffSerTol ) break; 
                }
        for ( ; i > 0; --i )  sum += terms[i];
        sum = (1.0+sum)/util_fact(k);
        return sum; }

int StumpffFunctions(double *c, double x) {
        double x2, tmp;
        int i;
/* What this is supposed to be
c[0] = StumpffSeriesCk(0, x);
c[1] = StumpffSeriesCk(1, x);
c[2] = StumpffSeriesCk(2, x);
c[3] = StumpffSeriesCk(3, x);
return 0;
*/
/* FOR REFERENCE, THIS COULD BE USED FOR LARGE VALUES, MAYBE 
        if ( x > STUMPFFPLAINLIM ) {
                x2 = sqrt(x);
                c[0] = cos(x2);
                c[1] = sin(x2)/x2;
                c[2] = (1.0-cos(x2))/x;
                c[3] = (x2-sin(x2))/(x*x2); 
                return 0; }
        if ( x < -STUMPFFPLAINLIM ) {
                x2 = sqrt(-x);
                c[0] = cosh(x2);
                c[1] = sinh(x2)/x2;
                c[2] = (cosh(x2)-1.0)/x;
                c[3] = (sinh(x2)-x2)/(x*x2); 
                return 0; }
*/
        if ( fabs(x) < STUMPFFSHORTTOL ) {
                c[2] = Stumpff2Short(x); c[3] = Stumpff3Short(x);
                c[1] = 1.0 - x*c[3]; c[0] = 1.0 - x*c[2];
                return 0; }
        if ( fabs(x) < STUMPFFRED ) {
                c[2] = HornerStumpff2(x); c[3] = HornerStumpff3(x); 
                c[1] = 1.0 - x*c[3]; c[0] = 1.0 - x*c[2];
                return 0; }
/* reduce value to small */
        x2 = x;
        for ( i = 1; i <= 20; ++i ) {
                x2 = x2/4.0;
                if ( fabs(x2) < STUMPFFRED ) break; 
                }
        if ( i > 20 ) {
                fprintf(FLOG, "%s"," StumpfReduction No Good \n"); 
                CheckStumpffErrorCount();
                fflush(FLOG);
                return 0; } 
/* original 
        c[2] = StumpffSeriesCk(2, x2);
        c[3] = StumpffSeriesCk(3, x2);
*/
        c[2] = HornerStumpff2(x2);
        c[3] = HornerStumpff3(x2);
        c[1] = 1.0 - x2*c[3];
        c[0] = 1.0 - x2*c[2];
        for ( ; i > 0; --i ) {
                tmp = c[2];
                c[2] = (c[1]*c[1])/2.0; 
                c[3] = (tmp+c[0]*c[3])/4.0; 
                c[1] = c[0]*c[1];
                c[0] = (2.0*c[0]*c[0])-1.0;
                }
        return i; }


/* ===============================================================================
*/

/* NOTE:
It seems that there might be a systematic error here, not sure yet.
The Newton solvers tend to end up systematically on a preferred side 
of the solution.  E.g., for x^2: if starting from above the solution,
it ends at a slightly larger value than the solution.  Starting from
 below, in the first step it will get to the previous case.  Since 
the solution cannot be computed to machine precision, we could have 
a systematic error.  Since the orbit spends more time near apocenter
than near pericenter, this seems to lead an increase in energy and 
decrease in angular momentum.  Perhaps this could be randomized with 
an additional bisection at the end?
*/

double QuarticStumpffSolverIterator(double mu, double alpha, double s, 
                double r0, double u, double deltat) {
        double c[5];
        double d0, d1, d2, d3, s2, ssav = s, tmp1, tmp2, tmp3;
        double delta1, delta2, delta3;
        int i;
for ( i = 1; i <= 20; ++i ) {
        s2 = s*s;
        StumpffFunctions(c, alpha*s2);
        tmp1 = s*c[1];
        tmp2 = s2*c[2];
        tmp3 = -r0*alpha+mu;
        d0 = u*tmp2+mu*s2*s*c[3];
        d0 += r0*tmp1;
        d0 -= deltat; 
        d1 = u*tmp1+mu*tmp2; 
        d1 += r0*c[0];
        d2 = u*c[0]+tmp3*tmp1; 
        d3 = tmp3*c[0]-u*alpha*c[1]; 
        delta1 = -d0/d1;
        delta2 = -d0/(d1+delta1*d2/2);
        delta3 = -d0/(d1+delta1*d2/2+delta2*delta2*d3/6);
        ssav= s;
        s = s + delta3;
        if ( fabs(delta3/s) < StumpffQSTol ) return s; }
        fprintf(FLOG, "%s"," QUARTIC STUMPFF SOLVER NOT CONVERGENT? \n"); 
        fprintf(FLOG, "    s = %.16e \n", s); 
        fprintf(FLOG, "sprev = %.16e \n", ssav); 
        fprintf(FLOG, "delta3/s = %.16e \n", delta3/s); 
        fprintf(FLOG, "Stumpff x = %.16e \n", alpha*s2); 
        fprintf(FLOG, "deltas = %.16e ", delta1); fprintf(FLOG, "%.16e ", delta2); 
        fprintf(FLOG, "%.16e \n", delta3);
        CheckStumpffErrorCount();
        fflush(FLOG);
        return s; }


double QuarticStumpffSolver(double mu, double alpha, double s, 
                double r0, double u, double deltat) {
/* initial guess */
        double tmp = deltat/r0;
        double s0 = tmp;
        double r0dot = u/r0;
        tmp = tmp*tmp;
        s0 -= 0.5*r0dot*tmp;
/* third order in deltat */
        tmp = tmp*deltat/r0;
        s0 += (0.5*r0dot*r0dot-1/6*(mu-alpha*r0)/r0)*tmp;
        s = QuarticStumpffSolverIterator( mu, alpha, s0, r0, u, deltat);
        return s; }

/* ===============================================================================
*/

void IncfgAndDerivsStumpff(double *vals, double mu, double alpha, double s, 
                double r0, double u, double deltat) {
        double r, s2, tmp1, tmp2, tmp3;
        double c[5];
        s = QuarticStumpffSolver( mu, alpha, s, r0, u, deltat);
        s2 = s*s;
        StumpffFunctions(c, alpha*s2);
/* recompute c[1] from angular mom. */
        tmp1 = s*c[1];
        tmp2 = s2*c[2];
        tmp3 = u*tmp1+mu*tmp2;
/* r */
        r = tmp3+(r0*c[0]);
        vals[0] = r;
        tmp3 = -mu*tmp2;
/* f */
        vals[1] = tmp3/r0;
/* g, this is usually bad in terms of roundoff */
        vals[2] = r0*tmp1+u*tmp2;
/* derivs */
        vals[3] = -(mu*tmp1)/(r0*r);
        vals[4] =  tmp3/r;
        return; }

void fgAndDerivsStumpff(double *vals, double mu, double alpha, double s, 
                double r0, double u, double deltat) {
        IncfgAndDerivsStumpff(vals, mu, alpha, s, r0, u, deltat);
        vals[1] += 1.0;
        vals[4] += 1.0;
        return; }

double KeplerStumpffStep(int dim, double mu, double *xv, double deltat) {
        double fgvals[5], tmp;
        double *x = xv;
        double *v = xv+dim;
        double v02 = 0.0, u = 0.0, r0 = 0.0;
        double r0dot, alpha;
        int i;
        for ( i = 1; i <= dim; ++i ) {
                v02 += v[i]*v[i];
                u += x[i]*v[i];
                r0 += x[i]*x[i]; }
        r0 = sqrt(r0);
        r0dot = u/r0;
/* a = 1.0/(2.0/r0 - v02/mu); alpha = mu/a; */
        alpha = 2.0*mu/r0 - v02;
        fgAndDerivsStumpff(fgvals, mu, alpha, 0.0, r0, u, deltat);
/* now update */
        for ( i = 1; i <= dim; ++i ) {
                tmp = x[i];
                x[i] = fgvals[1]*x[i]+fgvals[2]*v[i];
                v[i] = fgvals[3]*tmp+fgvals[4]*v[i]; }
        return 0.0; }


/* return only the increment */
double IncrementKeplerStumpffStep(int dim, double mu, double *xv, double deltat) {
        double fg[5], tmp, tmp2;
        double *x = xv;
        double *v = xv+dim;
        double alpha, v02 = 0.0, u = 0.0, r02 = 0.0, r0;
        int i, ii;
        double tmp1, r, delta1;
        for ( i = 1; i <= dim; ++i ) {
                v02 += v[i]*v[i];
                u += x[i]*v[i];
                r02 += x[i]*x[i]; }
        r0 = sqrt(r02);
/* a = 1.0/(2.0/r0 - v02/mu); alpha = mu/a; */
        alpha = 2.0*mu/r0 - v02;
        IncfgAndDerivsStumpff(fg, mu, alpha, 0.0, r0, u, deltat);

/* using fg'-gf' = 1 ( ' = dot) or total energy do not help much */
                
/* now update */
        for ( i = 1; i <= dim; ++i ) {
                tmp = x[i];
                x[i] = fg[1]*x[i]+fg[2]*v[i]; 
                v[i] = fg[3]*tmp+fg[4]*v[i];  
                }
        return 0.0; }

/* FOR TEST 
double IncrementKeplerStumpffStep(int dim, double mu, double *xv, double deltat) {
        double tmp[7], tmp2[7];
        int i;
        for ( i = 1; i <= 2*dim; ++i ) { tmp[i] = tmp2[i] = xv[i]; }
IncrementKeplerStumpffStepSub(dim, mu, tmp, deltat);
        for ( i = 1; i <= 2*dim; ++i ) tmp2[i] += tmp[i];
IncrementKeplerStumpffStepSub(dim, mu, tmp2, -deltat);
        for ( i = 1; i <= 2*dim; ++i ) xv[i] = (tmp[i]-tmp2[i])/2.0;
        return 0.0; }
*/