source: proiecte/NBody/NBI/NBIparalel.c @ 124

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

                        NBI 
  A set of numerical integrators for the gravitational N-body problem.

        (C) Copyright 1996 Ferenc Varadi.  All rights reserved.
  This copyright notice is intended to prevent the commercial use of this code.
  The code is free and can be freely distributed.
  Do not remove this copyright notice unless you make changes in the code.
  If you make changes in the code, add your name to the copyright notice.

  The usual disclaimers apply, i.e., use the code at your own risk.
  The author does not guarantee anything.

  This code was developed by F. Varadi in collaboration
  with M. Ghil, W. I. Newman, W. M. Kaula, K. Grazier, D. Goldstein 
  and M. Lessnick. The partial financial support of NSF Grant 
  ATM95-23787 (M. Ghil and F. Varadi) is gratefully acknowledged.

  If you have any questions, suggestions etc., contact F. Varadi, 
  e-mail: varadi@ucla.edu
  WWW: http://www.atmos.ucla.edu/~varadi

 
  For the impatient: 
        Save this file as NBI.c
        Compile the code as
        cc -o NBI NBI.c -lm

        Download the sample input file NBIsampleinput and save it as 
        NBIsampleinput.  It is an input file for the outer Solar System, with
        the planets from Jupiter through Neptune based on the Jet Propulsion 
        Laboratory's DE403 (Digital Ephemeris 403).  This was obtained 
        from the ftp site navigator.jpl.nasa.gov, updates can be accessed 
        there.

        Run the code:
        NBI NBIsampleinput
        Take a look at the output in NBIsampleoutput and the log file 
        NBIsampleoutput.log Change the MAXNBODIES macro, i.e., the 
        maximum number of particles, to suit your needs.  It is set to 
        2100 in this version.

1) Description:
        Four numerical integrators for the gravitational (nonrelativistic) 
        N-body problem are bundled into a single source code. The N bodies 
        are assumed to be point masses. Particles with zero mass are called
        test particles; they are assumed not to influence each other and 
        the particles with nonzero mass.

        The integrators are:

        a) A Wisdom-Holman-type mapping for hierarchical N-body systems.
                Code name: HWH, maximum global order = 4
           This can deal with more general arrangements than the standard
           Wisdom-Holman mapping.  Hierarchical N-body systems are described
           by Roy (1988). These are represented in the code by binary trees.
           The integrator takes advantage of certain properties of 
           symplectic forms and Jacobi coordinates, these will be described
           in a paper (in preparation). We also use singularly weighted
           symplectic forms, these are discussed by Varadi et al. (1995).

        b) A modified Cowell-Stormer integrator, with modifications by 
           W. I. Newman and his students. 
                Code name: CSN, max global order = 15
           Note: (global order) = (local order) - 2 (Goldstein, 1996)

        c) A Gragg-Bulirsch-Stoer integrator, as it is described by
           Hairer et al., (1993).
                Code name: GBS, max order = 20 (max stage = 10)
        IMPORTANT: The user has to specify the number of stages
        in the input file and not the order. Order = 2*(number of stages),
        see the Notes below.

        d) A symplectic mapping with Kinetic-Potential energy splitting
                Code name: SKP, max global order = 4 
           This is included for the sake of completeness.

        
2) References:
        
        Goldstein, D.: 1996, The Near-Optimality of Stormer
        Methods for Long Time Integrations of y''=g(y), 
        Ph.D. Dissertation, Univ. of California, Los Angeles,
        Dept. of Mathematics

        Hairer, E., Norsett, S. P. and Wanner, G.: 1993, 
        Solving Ordinary Differential Equations I.  Nonstiff Problems.
        Second Revised Edition
        
        Lessnick, M.: 1996, Stability Analysis of Symplectic
        Integration Schemes, Ph.D. Dissertation, Univ. of California, 
        Los Angeles, Dept. of Mathematics

        Roy, A. E.: 1988, Orbital Motion, Institute of Physics
        Publishing, Bristol
        
        Varadi, F. De la Barre, Kaula, W. M. and Ghil, M.: 1995,
        Singularly Weighted Symplectic Forms and Applications to 
        Asteroid Motion, Celestial Mechanics and Dynamical Astronomy,
        vol. 62, pp. 23-41

        Yoshida, H.: 1993, Recent Progress in The Theory and 
        Application of Symplectic Integrators,
        Celestial Mechanics and Dynamical Astronomy, vol. 56, pp. 27-43
        
        

3) Notes
        
        All long-term integrations suffer from roundoff errors.
        Symplectic schemes guarantee symplecticity for short 
        integrations but they are not strictly symplectic
        for very long integrations.  The implementation of 
        the Wisdom-Holman mapping in this code was aimed at
        reducing the effects of roundoff errors and thus make
        the scheme as symplectic as possible.
        
        The Cowell-Stormer integrator, being a multistep
        integrator, has to be intialized.  In theory, the starting 
        points should be computed with machine precision. This
        can be done using quadrupole precision arithmetic but
        the Gragg-Bulirsch-Stoer scheme in double precision 
        appers to be quite adequate. It should be noted that
        the stability of the Cowell-Stormer scheme ensures
        that initial errors do not grow due to computational
        limitations (Goldstein, 1996).

        In the case of the Gragg-Bulirsch-Stoer scheme one
        specifies the number of stages k in the extrapolation
        scheme. The order of the method is then 2*k.
        This factor of 2 comes from the fact that the 
        underlying method (Gragg's) has only even powers
        of the step size in the expansion of the local error
        (cf. Hairer et al., 1993, especially equation 9.22).

        All integrators use only G*mass, there is no need for
        the masses themselves nor the value of G.  Any set of units 
        can be used:
        distance_unit for position,
        time_unit for time and step size,
        distance_unit/time_unit for velocity.
        E.g., one can use AU, day and AU/day as units.

        None of the codes is regularized in any sense. We intend
        to provide regularized code in the future. For the Wisdom-Holman 
        mapping this involves changing the tree of Jacobian reference 
        frames and dealing with the transition region accurately. 

        In the case of the Wisdom-Holman mapping, the computation of 
        acceleration differences between Jacobi and inertial coordinates 
        could be improved in order to further reduce roundoff error.
        These differences can be represented by multi-pole expansions
        of the potential but we have not explored the details yet.

        The present version of the integrator does not monitor error.

        The choice of step size is up to the user. In most cases
        one has to adjust the step size to the shortest orbital
        period T in the system. The Wisdom-Holman mapping at second 
        order is usually used with step size of T/10-T/100.
        We recommend T/1000 for Cowell-Stormer integrator
        since this renders the calculation EXACT to double
        precision, even for high (e=0.5) eccentricities.
        The cost of using this short time step is substantially
        compensated by reduced computational complexity (as
        compared to the Wisdom-Holman mapping).  For much larger 
        step sizes the integrator is not stable.

        The Gragg-Bulirsch-Stoer method appears to be stable
        for large orders and step sizes. For order 20 one can
        even use T/5 as the step size and still get accurate results.
        Since this is an accurate but compuationally very expensive
        integrator, it is mainly used by us for short integrations.
        We have not investigated the propagation of roundoff 
        error in the GBS method.
        

4) How to compile:
        Since there is only a single source file, it is very easy to 
        compile the code with whatever C compiler is available. The 
        code was compiled and tested on Sun SPARCstations (Solaris 2.5.1 
        and earlier) and DEC Alpha workstations (Digital UNIX V3.2).
        For instance, on the SPARCstations,
        
        cc -fast -xO5 -xdepend -o NBI NBI.c -lm

        should work with Sun's C compiler.  The code uses only the standard 
        C library, more exactly: I/O functions, exit, malloc and free for 
        the tree, sqrt, cbrt, sin, cos, sinh, cosh.

        Caution: optimization may introduce problems with some compilers 
        and it is prudent to check the correctness of the optimized code 
        by compiling without optimization, e.g., 
        cc -g -o NBI_no_op NBI.c -lm
        By comparing the output from NBI and NBI_no_op one can detect if 
        the optimization leads to any problems.  (It should not but in 
        practice the situation is not so clear.)

        One also can play with various compiler switches which control 
        roundoff but we have not explored these.

5) How to run

        NBI inputfile

        (Everything has to be specified in inputfile.)

6) On hierarchical N-body systems

        We can give only a brief introduction into these; a detailed technical 
        description is in preparation.  Our approach is somewhat different 
        from that of Roy (1988). 

        One starts with the concept of joining subsystems and applies this 
        concept recursively to form new systems. When two subsystems are 
        joined, the two subsystems are assumed to behave as point masses
        and thus define the appropriate two-body problem for the Wisdom-Holman 
        mapping. The subsystems are, in general, not point masses and one 
        takes this into account in the perturbation step of the Wisdom-Holman 
        mapping. 

        Instead of a simple chain of Jacobi coordinates, the code uses a 
        tree of Jacobi coordinates. This makes it possible to use the code
        e.g., in situations when the massive particles have sattelites, either 
        massive ones or not. It also reduces, when carefully specified, the 
        integration error since more of the perturbations are absorbed into 
        two-body interactions. The tree has to be specified in the input file.  
        One can define this tree formally as it is processed by the code which 
        works like a simple finite automaton but we did not use lex, yacc or 
        any other compiler-compiler tools.

 Tree description rules:
 Rule 1) particles are indexed by consecutive natural numbers, from 1 to N,
 Rule 2) particles are subsystems,
 Rule 3) (x y) : join the subsystems x and y to create a new subsystem.,
 Rule 4) [i1 - i2] denotes a set of test particles, from index i1 to index i2,
        which are treated the same way,
        Note: the space before and after - is necessary!
 Rule 5) ([i1 - i2] i3) is forbidden, use instead (i3 [i1 - i2]) 

 Example: 2000 main belt asteroids, from  6 to 2006, and
        Sun, Jupiter, Saturn, Uranus and Neptune (1 - 5) 
        can described as
 (((((1 [6 - 2006]) 2) 3) 4) 5)

 Another example for the tree specification when test particles between 
 Jupiter and Saturn are added with indeces from 6 to 1000 and another set 
 of test particles is added between Saturn and Uranus with indices from 
 1001 to 2000. The indices for the massive bodies are:
 Sun 1, Jupiter, 2, Saturn 3, Uranus 4, Neptune 5.
((((((1 2) [6 - 1000]) 3) [1001 - 2000]) 4) 5)

 Yet another example: the Sun, Earth, Moon, Jupiter and Io with indeces 
 1,2,3,4 and 5, respectively, for which the tree is
 ((1 (2 3)) (4 5))


7) Format of the input file
Id method order step_size  total_integration_time printing_time 
        number_of_particles

 tree_description

GMs of particles, possibly zeros

Initial conditions: x,y,z, vx, vy, vz

Explanation:

        Id is an arbitrary string which specifies the name of the experiment.

        Method can be SKP, HWH, CSN and GBS.

        Tree_description is used only when the method is HWH or the HWH is used
        to initialize CSN. One can put a dummy tree description of the form
        (1 2)
        in the input file when HWH is not used.  This will satisfy the parser 
        which processes the tree description but will not play any role in the 
        integration.


Example: 
 For Sun, Jupiter, Saturn, Uranus and Neptune, 
 with indeces 1,2,3,4 and 5, in heliocentric coordinates.

SomeId HWH 4 100.0 365.0e+4 1.0e+4 
5
((((1 2) 3) 4) 5)

0.295912208285591095e-03   
0.282534590952422643e-06  0.845971518568065874e-07  
0.129202491678196939e-07   0.152435890078427628e-07  

0.0 0.0 0.0 0.0 0.0 0.0

-0.538420940678015203e+01   -0.831247035699103631e+00 -0.225095848657298592e+00  
0.109236311953937086e-02   -0.652329334256263223e-02 -0.282301443780311710e-02  

0.788989169798583756e+01   0.459570785756998301e+01  0.155843220515231762e+01 
-0.321720261683698921e-02   0.433063255278440598e-02  0.192641782724487106e-02 

-0.182698980946940850e+02   -0.116271406901560681e+01 -0.250371938307287545e+00  
0.221540447608191936e-03   -0.376765389967669917e-02 -0.165324449144407023e-02 

-0.160595376822070470e+02   -0.239429495784517528e+02 -0.940042953888094424e+01  
0.264312302715172384e-02   -0.150349219713705076e-02 -0.681271071793930938e-03 


There are 2 output files:
        1) SomeId contains the actual integration data
        2) SomeId.log contains a copy of the input data and error messages.
        Some error messages might be sent to the standard output and the user 
        has to redirect it in order to capture them.

In the case of the Cowell-Stormer integrator, an additional file, SomeId.lerr 
is created. This contains the LOCAL truncation error of the integration, 
computed for each step but only printed for the last step in the printing 
cycle. This also should help to determine the appropriate step size.

The output is in the original coordinate system, i.e., inertial and not Jacobian.
Output is produced when time is an integer multiple of the specified printing 
time. The following is printed:
time    relative_change_in_total_energy relative_change_in_ang.mom.z.comp.
after which
x y z 
dx dy dz 
for each particle.


9) The uset can change the code, of course. The section MAIN LOOP is where 
   this can be done easily. You can insert code to do more things, e.g., 
   detecting close approaches.

10) MACROS TO GIVE SOME CONTROL OVER THINGS.
Changes in the code that you can make HERE:
        CHANGE the DIM macro to 2 for planar motions.
        (Note: a variables dim is also used, easy to change
                it everywhere to DIM. It is still better to keep
                the Kepler solver in the more general 3-dimensional form.)
        CHANGE the MAXNBODIES macro to whatever you want,
                including test particles.

**************** END OF DESCRIPTION ********************** */

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        USER-CHANGEABLE MACROS    */

/* 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


/* VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV */

/*      The actual code.
        Change things below this line at your own peril.
        
 */

/* MACROS TO DISTINGUISH BETWEEN TEST PARTICLES AND
   PARTICLES WITH FINITE MASS
   You can add more types here, these values are stored
        in tpstore for each particle
   e.g., you could define stopped test particles,
   but you have to change the code accordingly.
   It is perhaps easier to add more integer arrays.
   IMPORTANT: THE TREE MAINTAINS ITS OWN INTEGER ARRAYS
   WHICH HOLD TEST PARTICLE INDICES. 
*/

#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>
#include "omp.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; }

// X
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; 
}

// X
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; 
}

//X 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;
                
                #pragma omp parallel for reduction(+:tmp) schedule(runtime)
                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;
                        #pragma omp parallel for reduction(+:tmp) schedule(runtime)
                        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; 
}

// X 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;
                // x*vy - y*vx  
                az = az + (p[i][1]*p[i][2+dim] - p[i][2]*p[i][1+dim]) * p[i][GMASSCO]; 
        }
        return az; 
}

// X transforms all particles to coordinates relative to barycenter 
int tobarycenter(int nump, int dim, double p DECLARR) {
        double bc[10];
        int i, j;
        
        #pragma omp parallel private(j)
        #pragma omp for schedule(runtime)
        for (j = 0; j <= 2*dim; ++j ) 
                bc[j] = 0.0;
        
        for ( i = nump; i >= 1; --i ) {
                bc[0] += p[i][GMASSCO];
                #pragma omp parallel private(j)
                #pragma omp for schedule(runtime)
                for (j = 1; j <= 2*dim; ++j ) 
                        bc[j] += p[i][j]*p[i][GMASSCO]; 
        }
        
        #pragma omp parallel private(j)
        #pragma omp for schedule(runtime)
        for (j = 1; j <= 2*dim; ++j ) 
                bc[j] = bc[j]/bc[0];
                
    
        for ( i = nump; i >= 1; --i ) {
                #pragma omp parallel private(j)
                #pragma omp for schedule(runtime)
                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
        // id - name of output 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
        */
        
        // the length of the inner loop, i.e.,  between printing 
        nin = (int) (fabs(ptime)/fabs(stepsize));  
        
        if ( nin == 0 ) {
                fprintf(flog, "%s\n", " Stepsize is larger than printing time, aborted");
                exit(1); 
        }
        
        // the length of the outer loop, i.e., printings from beginning to end 
        nps = (int) (fabs(ftime)/fabs(ptime));  

        fprintf(flog, "nin = %d nps = %d", nin, nps);
        
        /* 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;
                fprintf(f, "end of step, time = %f\n", time);
                
                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 = cbrt(2.0);
        Coeff4x0 = -cbrt2/(2.0-cbrt2);
        Coeff4x1 = 1.0/(2.0-cbrt2);
        SympCalled = 1; }
        

static int HWHCalled = 0;

// X
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
*/

// X 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;
        
        #pragma omp parallel for private(j) schedule(runtime)
        for ( j = 1; j <= DIM; ++j ) 
                tree->acc[j] = 0.0;
        
        tree->tmass = 0.0; tree->rm1 = 0.0; tree->rm2 = 0.0;
        #pragma omp parallel for private(j) schedule(runtime)
        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;
                                #pragma omp parallel for private(i) schedule(runtime)
                                for ( i = 1; i <= DIM; ++i ) 
                                        p[ind][ACC+i] = 0.0;
                                } 
                        }  
            else { 
                        tree->tmass = p[ind][GMASSCO];
                        #pragma omp parallel for private(j) schedule(runtime)
                        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);
        
        #pragma omp parallel for private(j) schedule(runtime)
        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);

// X 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;
        #pragma omp parallel for private(i) schedule(runtime)
        for ( i = 1; i <= 2*DIM; ++i ) 
                cmove[i] = 0.0;

        // move by 1/2*(Tmass)*v^2 */
        #pragma omp parallel for private(i) schedule(runtime)
        for ( i = 1; i <= DIM; ++i ) 
                cmove[i] = (tree->bc[DIM+i])*stepsize; 

        MoveCMandKepler(tree, cmove, stepsize, p);
        return 0; 
}


// X 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;
                        #pragma omp parallel for private(i) schedule(runtime)
                        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; 
}

// X 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;
        #pragma omp parallel for private(i) schedule(runtime)
        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); 
        #pragma omp parallel for private(i) schedule(runtime)
        // 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);
        #pragma omp parallel for private(i) schedule(runtime)
        for ( i = 1; i <= DIM2; ++i ) { 
                tmp[i] = cmove[i] + tmpkep[i]*(tree->rm1); 
        }
        MoveCMandKepler(tree->p2, tmp, stepsize, p);
        return 0; 
}
        
/*  ++++++++++++++++++++++++++++++++++++++++++++++++++++++
        PERTURBATION ACCELERATIONS
*/

/* X 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
*/
// X 
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 ) {
                #pragma omp parallel private(j)
                #pragma omp for schedule(runtime)
                for ( j = 1; j <= DIM2; ++j ) { 
                        bc[j] = tree->p1->bc[j]; 
                }
                return 0; 
        }
        
        #pragma omp parallel private(j)
        #pragma omp for schedule(runtime)
        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; }
                        

/* X 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; }

// X
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?
*/
// X
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; 
}

// X
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; 
}

/* ===============================================================================
*/
// X
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; }


// X 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; }
*/