Index: proiecte/NBody/NBI/NBIparalel.c
===================================================================
--- proiecte/NBody/NBI/NBIparalel.c	(revision 124)
+++ proiecte/NBody/NBI/NBIparalel.c	(revision 124)
@@ -0,0 +1,2113 @@
+
+
+#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; }
+*/
+
+	
+
+
Index: proiecte/NBody/NBI/NBIsampleinput
===================================================================
--- proiecte/NBody/NBI/NBIsampleinput	(revision 124)
+++ proiecte/NBody/NBI/NBIsampleinput	(revision 124)
@@ -0,0 +1,25 @@
+
+
+NBIsampleoutput HWH 2 100.0 1.0e+3 100.0 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 
+
Index: proiecte/NBody/NBI/NBIsampleinputCSN
===================================================================
--- proiecte/NBody/NBI/NBIsampleinputCSN	(revision 124)
+++ proiecte/NBody/NBI/NBIsampleinputCSN	(revision 124)
@@ -0,0 +1,25 @@
+
+
+NBIsampleoutputCSN CSN 14 5.0 1.0e+3 100.0 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 
+
Index: proiecte/NBody/NBI/NBIsampleoutput
===================================================================
--- proiecte/NBody/NBI/NBIsampleoutput	(revision 124)
+++ proiecte/NBody/NBI/NBIsampleoutput	(revision 124)
@@ -0,0 +1,131 @@
+0.000000000000   -0.0000000000000000e+00  0.0000000000000000e+00
+4.50417099317595e-03  7.62961724685589e-04  2.64217371485701e-04  
+-2.68697979929186e-07  5.22529622296852e-06  2.24893094591555e-06  
+-5.37970523578698e+00  -8.30484073974418e-01  -2.24831631285813e-01  
+1.09209442155944e-03  -6.51806804633966e-03  -2.82076550685720e-03  
+7.89439586897901e+00  4.59647081929467e+00  1.55869642252380e+00  
+-3.21747131481692e-03  4.33585784900737e-03  1.92866675819079e-03  
+-1.82653939237009e+01  -1.16195110729092e+00  -2.50107720935802e-01  
+2.21271749628263e-04  -3.76242860345373e-03  -1.65099556049815e-03  
+-1.60550335112139e+01  -2.39421866167271e+01  -9.40016532150946e+00  
+2.64285432917179e-03  -1.49826690091408e-03  -6.79022140848015e-04  
+end of step, time = 100.000000
+100.000000000000   3.1833330095145376e-06  0.0000000000000000e+00
+4.43470163799496e-03  1.27881617027082e-03  4.87169548090825e-04  
+-1.11650683486053e-06  5.07395274516697e-06  2.20235889538463e-06  
+-5.22177706769564e+00  -1.47273457596714e+00  -5.03999255723156e-01  
+2.06158745571407e-03  -6.30740421968653e-03  -2.75409259857205e-03  
+7.55812238591514e+00  5.02121264801879e+00  1.74853560928254e+00  
+-3.50631128748689e-03  4.15607858332428e-03  1.86684693053287e-03  
+-1.82388583447985e+01  -1.53788291438590e+00  -4.15133513249719e-01  
+3.09408166782867e-04  -3.75590576521328e-03  -1.64938766719185e-03  
+-1.57898995885472e+01  -2.40907383549782e+01  -9.46756680406041e+00  
+2.65977732482993e-03  -1.47274129942280e-03  -6.68995472394411e-04  
+end of step, time = 200.000000
+200.000000000000   5.6476979739305321e-06  0.0000000000000000e+00
+4.28193883600039e-03  1.77424193093240e-03  7.03157606181118e-04  
+-1.93224968021170e-06  4.81733658902443e-06  2.10985689790137e-06  
+-4.96861452881590e+00  -2.08812008838431e+00  -7.73972057404004e-01  
+2.99428784504893e-03  -5.98143590405836e-03  -2.63709446480907e-03  
+7.19349357674360e+00  5.42711437992709e+00  1.93181408062007e+00  
+-3.78435266814374e-03  3.95913801957762e-03  1.79747667516786e-03  
+-1.82035194244026e+01  -1.91307201387244e+00  -5.79958766677753e-01  
+3.97330571484281e-04  -3.74757557370095e-03  -1.64698516089011e-03  
+-1.55230874357681e+01  -2.42367296075453e+01  -9.53396202565169e+00  
+2.67641863579295e-03  -1.44705762648138e-03  -6.58897094085098e-04  
+end of step, time = 300.000000
+300.000000000000   7.3865420854433705e-06  0.0000000000000000e+00
+4.04977926155266e-03  2.23891489506285e-03  9.07659296215067e-04  
+-2.70217587368023e-06  4.45970140834114e-06  1.97292283353421e-06  
+-4.62463967265561e+00  -2.66533079887110e+00  -1.02979398959264e+00  
+3.87539384423796e-03  -5.54478047979900e-03  -2.47138966745866e-03  
+6.80165661700951e+00  5.81248621792361e+00  2.10778477369029e+00  
+-4.05025181993996e-03  3.74557987249766e-03  1.72072275087504e-03  
+-1.81594009612966e+01  -2.28733802736740e+00  -7.44504142277774e-01  
+4.84991063077416e-04  -3.73744536956626e-03  -1.64379057782127e-03  
+-1.52546253095784e+01  -2.43801446993267e+01  -9.59934386714829e+00  
+2.69277650197591e-03  -1.42121852563903e-03  -6.48728044142499e-04  
+end of step, time = 400.000000
+400.000000000000   8.4014629122640643e-06  2.0004763764791396e-16
+3.74347935687927e-03  2.66299423450318e-03  1.09632737559165e-03  
+-3.41284209252147e-06  4.00646796653433e-06  1.79356273648099e-06  
+-4.19573841296085e+00  -3.19358082974420e+00  -1.26669792167572e+00  
+4.69043198560084e-03  -5.00330847735433e-03  -2.25914196690840e-03  
+6.38389206250193e+00  6.17569787638520e+00  2.27571946621414e+00  
+-4.30268655785009e-03  3.51605083979174e-03  1.63679536241057e-03  
+-1.81065315285780e+01  -2.66050137086335e+00  -9.08690581368448e-01  
+5.72342063111074e-04  -3.72552369742407e-03  -1.63980698608947e-03  
+-1.49845416413168e+01  -2.45209682203765e+01  -9.66370531359938e+00  
+2.70884919076286e-03  -1.39522665642566e-03  -6.38489368013643e-04  
+end of step, time = 500.000000
+500.000000000000   8.6976174729724906e-06  2.0004763764791396e-16
+3.36962366834482e-03  3.03724228814317e-03  1.26504180408840e-03  
+-4.05113651632841e-06  3.46429285858795e-06  1.57432078829002e-06  
+-3.68922655739975e+00  -3.66273688685910e+00  -1.48016161591597e+00  
+5.42528568466351e-03  -4.36421571732426e-03  -2.00309194416796e-03  
+5.94161116049950e+00  6.51518891608248e+00  2.43491315956388e+00  
+-4.54036710216602e-03  3.27130114836189e-03  1.54594885686208e-03  
+-1.80449444392974e+01  -3.03238337442468e+00  -1.07243935844040e+00  
+6.59336358016544e-04  -3.71182029159280e-03  -1.63503798003310e-03  
+-1.47128650342026e+01  -2.46591850274917e+01  -9.72703945495643e+00  
+2.72463499717430e-03  -1.36908469412638e-03  -6.28182118282058e-04  
+end of step, time = 600.000000
+600.000000000000   8.2807943701113841e-06  2.0004763764791396e-16
+2.93608889801093e-03  3.35315278045891e-03  1.40996561078366e-03  
+-4.60434235294456e-06  2.84116522490936e-06  1.31832226209558e-06  
+-3.11381101924540e+00  -4.06345563725238e+00  -1.66596756357677e+00  
+6.06626517692627e-03  -3.63612441830140e-03  -1.70660204828907e-03  
+5.47635204349755e+00  6.82947907616276e+00  2.58468850866821e+00  
+-4.76204745011304e-03  3.01218393305239e-03  1.44848195590137e-03  
+-1.79746777078709e+01  -3.40280640038206e+00  -1.23567213346957e+00  
+7.45927140961730e-04  -3.69634606032894e-03  -1.62948767391305e-03  
+-1.44396242605648e+01  -2.47947802457757e+01  -9.78933948678235e+00  
+2.74013224403298e-03  -1.34279532955478e-03  -6.17807354579075e-04  
+end of step, time = 700.000000
+700.000000000000   7.1564067506393124e-06  0.0000000000000000e+00
+2.45199913734639e-03  3.60308975067932e-03  1.52760556063849e-03  
+-5.06025261559541e-06  2.14652251879342e-06  1.02932587695307e-06  
+-2.47954144914123e+00  -4.38733226741861e+00  -1.82026784744986e+00  
+6.60023120894513e-03  -2.82920388460447e-03  -1.37371134751867e-03  
+4.98977477156768e+00  7.11717848187924e+00  2.72440025009881e+00  
+-4.96653700255652e-03  2.73965335984200e-03  1.34473748055526e-03  
+-1.78957740073504e+01  -3.77159395987703e+00  -1.39831100356664e+00  
+8.32068052806902e-04  -3.67911306860361e-03  -1.62316069494958e-03  
+-1.41648482590527e+01  -2.49277392701806e+01  -9.85059871095189e+00  
+2.75533928212815e-03  -1.31636126882479e-03  -6.07366143494546e-04  
+end of step, time = 800.000000
+800.000000000000   5.3300551415920242e-06  0.0000000000000000e+00
+1.92766658849532e-03  3.78043852249901e-03  1.61487832898167e-03  
+-5.40734736839962e-06  1.39137323430739e-06  7.11780652300366e-07  
+-1.79774631779193e+00  -4.62706157647161e+00  -1.93965475682201e+00  
+7.01478417686661e-03  -1.95529808833808e-03  -1.00919483119124e-03  
+4.48365520688614e+00  7.37699759900516e+00  2.85343957575419e+00  
+-5.15271226164252e-03  2.45476142740796e-03  1.23510153329619e-03  
+-1.78082806226468e+01  -4.13857082761334e+00  -1.56027855389842e+00  
+9.17713222091436e-04  -3.66013451946995e-03  -1.61606217572835e-03  
+-1.38885661318321e+01  -2.50580477670251e+01  -9.91081053634296e+00  
+2.77025449037808e-03  -1.28978523312122e-03  -5.96859558486974e-04  
+end of step, time = 900.000000
+900.000000000000   2.8090517521581794e-06  0.0000000000000000e+00
+1.37451102565442e-03  3.87976864822592e-03  1.66918231943803e-03  
+-5.63504308914401e-06  5.88410657355560e-07  3.70880666697070e-07  
+-1.08094638350058e+00  -4.77661149471615e+00  -2.02123728018911e+00  
+7.29852875779565e-03  -1.02804344201647e-03  -6.18620346692741e-04  
+3.95987772329120e+00  7.60775680054539e+00  2.97123839652598e+00  
+-5.31952839775412e-03  2.15865340727430e-03  1.12000211235745e-03  
+-1.77122493998059e+01  -4.50356315467677e+00  -1.72149790781927e+00  
+1.00281730398437e-03  -3.63942473407376e-03  -1.60819774599750e-03  
+-1.36108071417633e+01  -2.51856916754904e+01  -9.96996847951898e+00  
+2.78487627599101e-03  -1.26306995846903e-03  -5.86288679793060e-04  
+end of step, time = 1000.000000
+1000.000000000000   -3.9565189531119833e-07  0.0000000000000000e+00
+8.04950420409227e-04  3.89700695119693e-03  1.68847450520941e-03  
+-5.73402078130268e-06  -2.47902691920044e-07  1.26091677916143e-08  
+-3.42737616773387e-01  -4.83140702164559e+00  -2.06272280903873e+00  
+7.44142091582259e-03  -6.29554642815864e-05  -2.08394209498509e-04  
+3.42042677457069e+00  7.80839540884993e+00  3.07727343840199e+00  
+-5.46603047146411e-03  1.85256191272964e-03  9.99907145008186e-04  
+-1.76077366914428e+01  -4.86639857929047e+00  -1.88189277615268e+00  
+1.08733551813521e-03  -3.61699913036238e-03  -1.59957352387901e-03  
+-1.33316007095650e+01  -2.53106572090924e+01  -1.00280661654022e+01  
+2.79920307462434e-03  -1.23621819550115e-03  -5.75654594336657e-04  
Index: proiecte/NBody/NBI/NBIsampleoutput.log
===================================================================
--- proiecte/NBody/NBI/NBIsampleoutput.log	(revision 124)
+++ proiecte/NBody/NBI/NBIsampleoutput.log	(revision 124)
@@ -0,0 +1,39 @@
+ID = NBIsampleoutput
+ORDER =  2
+STEPSIZE = 100.000000000000
+FINAL TIME = 1000.000000000000
+PRINTING TIME = 100.000000000000
+NUMBER OF PARTS =  5
+GMs
+2.95912208285591e-04  
+2.82534590952423e-07  
+8.45971518568066e-08  
+1.29202491678197e-08  
+1.52435890078428e-08  
+Inits 
+0.00000000000000e+00  0.00000000000000e+00  0.00000000000000e+00  
+0.00000000000000e+00  0.00000000000000e+00  0.00000000000000e+00  
+-5.38420940678015e+00  -8.31247035699104e-01  -2.25095848657299e-01  
+1.09236311953937e-03  -6.52329334256263e-03  -2.82301443780312e-03  
+7.88989169798584e+00  4.59570785756998e+00  1.55843220515232e+00  
+-3.21720261683699e-03  4.33063255278441e-03  1.92641782724487e-03  
+-1.82698980946941e+01  -1.16271406901561e+00  -2.50371938307288e-01  
+2.21540447608192e-04  -3.76765389967670e-03  -1.65324449144407e-03  
+-1.60595376822070e+01  -2.39429495784518e+01  -9.40042953888094e+00  
+2.64312302715172e-03  -1.50349219713705e-03  -6.81271071793931e-04  
+Actual Barycentric Inits 
+4.50417099317595e-03  7.62961724685589e-04  2.64217371485701e-04  
+-2.68697979929186e-07  5.22529622296852e-06  2.24893094591555e-06  
+-5.37970523578698e+00  -8.30484073974418e-01  -2.24831631285813e-01  
+1.09209442155944e-03  -6.51806804633966e-03  -2.82076550685720e-03  
+7.89439586897901e+00  4.59647081929467e+00  1.55869642252380e+00  
+-3.21747131481692e-03  4.33585784900737e-03  1.92866675819079e-03  
+-1.82653939237009e+01  -1.16195110729092e+00  -2.50107720935802e-01  
+2.21271749628263e-04  -3.76242860345373e-03  -1.65099556049815e-03  
+-1.60550335112139e+01  -2.39421866167271e+01  -9.40016532150946e+00  
+2.64285432917179e-03  -1.49826690091408e-03  -6.79022140848015e-04  
+Actual integrator is HWH 
+CODE: NBI 
+VERSION: 1 
+nin = 1 nps = 10
+ Finished 
Index: proiecte/NBody/NBI/NBIsampleoutput1
===================================================================
--- proiecte/NBody/NBI/NBIsampleoutput1	(revision 124)
+++ proiecte/NBody/NBI/NBIsampleoutput1	(revision 124)
@@ -0,0 +1,131 @@
+0.000000000000   -0.0000000000000000e+00  0.0000000000000000e+00
+4.50417099317595e-03  7.62961724685589e-04  2.64217371485701e-04  
+-2.68697979929186e-07  5.22529622296852e-06  2.24893094591555e-06  
+-5.37970523578698e+00  -8.30484073974418e-01  -2.24831631285813e-01  
+1.09209442155944e-03  -6.51806804633966e-03  -2.82076550685720e-03  
+7.89439586897901e+00  4.59647081929467e+00  1.55869642252380e+00  
+-3.21747131481692e-03  4.33585784900737e-03  1.92866675819079e-03  
+-1.82653939237009e+01  -1.16195110729092e+00  -2.50107720935802e-01  
+2.21271749628263e-04  -3.76242860345373e-03  -1.65099556049815e-03  
+-1.60550335112139e+01  -2.39421866167271e+01  -9.40016532150946e+00  
+2.64285432917179e-03  -1.49826690091408e-03  -6.79022140848015e-04  
+end of step, time = 100.000000
+100.000000000000   3.1833330095145376e-06  0.0000000000000000e+00
+4.43470163799496e-03  1.27881617027082e-03  4.87169548090825e-04  
+-1.11650683486053e-06  5.07395274516697e-06  2.20235889538463e-06  
+-5.22177706769564e+00  -1.47273457596714e+00  -5.03999255723156e-01  
+2.06158745571407e-03  -6.30740421968653e-03  -2.75409259857205e-03  
+7.55812238591514e+00  5.02121264801879e+00  1.74853560928254e+00  
+-3.50631128748689e-03  4.15607858332428e-03  1.86684693053287e-03  
+-1.82388583447985e+01  -1.53788291438590e+00  -4.15133513249719e-01  
+3.09408166782867e-04  -3.75590576521328e-03  -1.64938766719185e-03  
+-1.57898995885472e+01  -2.40907383549782e+01  -9.46756680406041e+00  
+2.65977732482993e-03  -1.47274129942280e-03  -6.68995472394411e-04  
+end of step, time = 200.000000
+200.000000000000   5.6476979739305321e-06  0.0000000000000000e+00
+4.28193883600039e-03  1.77424193093240e-03  7.03157606181118e-04  
+-1.93224968021170e-06  4.81733658902443e-06  2.10985689790137e-06  
+-4.96861452881590e+00  -2.08812008838431e+00  -7.73972057404004e-01  
+2.99428784504893e-03  -5.98143590405836e-03  -2.63709446480907e-03  
+7.19349357674360e+00  5.42711437992709e+00  1.93181408062007e+00  
+-3.78435266814374e-03  3.95913801957762e-03  1.79747667516786e-03  
+-1.82035194244026e+01  -1.91307201387244e+00  -5.79958766677753e-01  
+3.97330571484281e-04  -3.74757557370095e-03  -1.64698516089011e-03  
+-1.55230874357681e+01  -2.42367296075453e+01  -9.53396202565169e+00  
+2.67641863579295e-03  -1.44705762648138e-03  -6.58897094085098e-04  
+end of step, time = 300.000000
+300.000000000000   7.3865420854433705e-06  0.0000000000000000e+00
+4.04977926155266e-03  2.23891489506285e-03  9.07659296215067e-04  
+-2.70217587368023e-06  4.45970140834114e-06  1.97292283353421e-06  
+-4.62463967265561e+00  -2.66533079887110e+00  -1.02979398959264e+00  
+3.87539384423796e-03  -5.54478047979900e-03  -2.47138966745866e-03  
+6.80165661700951e+00  5.81248621792361e+00  2.10778477369029e+00  
+-4.05025181993996e-03  3.74557987249766e-03  1.72072275087504e-03  
+-1.81594009612966e+01  -2.28733802736740e+00  -7.44504142277774e-01  
+4.84991063077416e-04  -3.73744536956626e-03  -1.64379057782127e-03  
+-1.52546253095784e+01  -2.43801446993267e+01  -9.59934386714829e+00  
+2.69277650197591e-03  -1.42121852563903e-03  -6.48728044142499e-04  
+end of step, time = 400.000000
+400.000000000000   8.4014629120944070e-06  2.0004763764791396e-16
+3.74347935687927e-03  2.66299423450318e-03  1.09632737559165e-03  
+-3.41284209252147e-06  4.00646796653433e-06  1.79356273648099e-06  
+-4.19573841296085e+00  -3.19358082974420e+00  -1.26669792167572e+00  
+4.69043198560084e-03  -5.00330847735433e-03  -2.25914196690840e-03  
+6.38389206250193e+00  6.17569787638520e+00  2.27571946621414e+00  
+-4.30268655785009e-03  3.51605083979174e-03  1.63679536241057e-03  
+-1.81065315285780e+01  -2.66050137086335e+00  -9.08690581368448e-01  
+5.72342063111074e-04  -3.72552369742407e-03  -1.63980698608947e-03  
+-1.49845416413168e+01  -2.45209682203765e+01  -9.66370531359938e+00  
+2.70884919076286e-03  -1.39522665642566e-03  -6.38489368013643e-04  
+end of step, time = 500.000000
+500.000000000000   8.6976174731421496e-06  2.0004763764791396e-16
+3.36962366834482e-03  3.03724228814317e-03  1.26504180408840e-03  
+-4.05113651632841e-06  3.46429285858795e-06  1.57432078829002e-06  
+-3.68922655739975e+00  -3.66273688685910e+00  -1.48016161591597e+00  
+5.42528568466351e-03  -4.36421571732426e-03  -2.00309194416796e-03  
+5.94161116049950e+00  6.51518891608248e+00  2.43491315956388e+00  
+-4.54036710216602e-03  3.27130114836189e-03  1.54594885686208e-03  
+-1.80449444392974e+01  -3.03238337442468e+00  -1.07243935844040e+00  
+6.59336358016544e-04  -3.71182029159280e-03  -1.63503798003310e-03  
+-1.47128650342026e+01  -2.46591850274917e+01  -9.72703945495643e+00  
+2.72463499717430e-03  -1.36908469412638e-03  -6.28182118282058e-04  
+end of step, time = 600.000000
+600.000000000000   8.2807943701113841e-06  2.0004763764791396e-16
+2.93608889801093e-03  3.35315278045891e-03  1.40996561078366e-03  
+-4.60434235294456e-06  2.84116522490936e-06  1.31832226209558e-06  
+-3.11381101924540e+00  -4.06345563725238e+00  -1.66596756357677e+00  
+6.06626517692627e-03  -3.63612441830140e-03  -1.70660204828907e-03  
+5.47635204349755e+00  6.82947907616276e+00  2.58468850866821e+00  
+-4.76204745011304e-03  3.01218393305239e-03  1.44848195590137e-03  
+-1.79746777078709e+01  -3.40280640038206e+00  -1.23567213346957e+00  
+7.45927140961730e-04  -3.69634606032894e-03  -1.62948767391305e-03  
+-1.44396242605648e+01  -2.47947802457757e+01  -9.78933948678235e+00  
+2.74013224403298e-03  -1.34279532955478e-03  -6.17807354579075e-04  
+end of step, time = 700.000000
+700.000000000000   7.1564067508089714e-06  0.0000000000000000e+00
+2.45199913734639e-03  3.60308975067932e-03  1.52760556063849e-03  
+-5.06025261559541e-06  2.14652251879342e-06  1.02932587695307e-06  
+-2.47954144914123e+00  -4.38733226741861e+00  -1.82026784744986e+00  
+6.60023120894513e-03  -2.82920388460447e-03  -1.37371134751867e-03  
+4.98977477156768e+00  7.11717848187924e+00  2.72440025009881e+00  
+-4.96653700255652e-03  2.73965335984200e-03  1.34473748055526e-03  
+-1.78957740073504e+01  -3.77159395987703e+00  -1.39831100356664e+00  
+8.32068052806902e-04  -3.67911306860361e-03  -1.62316069494958e-03  
+-1.41648482590527e+01  -2.49277392701806e+01  -9.85059871095189e+00  
+2.75533928212815e-03  -1.31636126882479e-03  -6.07366143494546e-04  
+end of step, time = 800.000000
+800.000000000000   5.3300551417616823e-06  0.0000000000000000e+00
+1.92766658849532e-03  3.78043852249901e-03  1.61487832898167e-03  
+-5.40734736839962e-06  1.39137323430739e-06  7.11780652300366e-07  
+-1.79774631779193e+00  -4.62706157647161e+00  -1.93965475682201e+00  
+7.01478417686661e-03  -1.95529808833808e-03  -1.00919483119124e-03  
+4.48365520688614e+00  7.37699759900516e+00  2.85343957575419e+00  
+-5.15271226164252e-03  2.45476142740796e-03  1.23510153329619e-03  
+-1.78082806226468e+01  -4.13857082761334e+00  -1.56027855389842e+00  
+9.17713222091436e-04  -3.66013451946995e-03  -1.61606217572835e-03  
+-1.38885661318321e+01  -2.50580477670251e+01  -9.91081053634296e+00  
+2.77025449037808e-03  -1.28978523312122e-03  -5.96859558486974e-04  
+end of step, time = 900.000000
+900.000000000000   2.8090517521581794e-06  0.0000000000000000e+00
+1.37451102565442e-03  3.87976864822592e-03  1.66918231943803e-03  
+-5.63504308914401e-06  5.88410657355560e-07  3.70880666697070e-07  
+-1.08094638350058e+00  -4.77661149471615e+00  -2.02123728018911e+00  
+7.29852875779565e-03  -1.02804344201647e-03  -6.18620346692741e-04  
+3.95987772329120e+00  7.60775680054539e+00  2.97123839652598e+00  
+-5.31952839775412e-03  2.15865340727430e-03  1.12000211235745e-03  
+-1.77122493998059e+01  -4.50356315467677e+00  -1.72149790781927e+00  
+1.00281730398437e-03  -3.63942473407376e-03  -1.60819774599750e-03  
+-1.36108071417633e+01  -2.51856916754904e+01  -9.96996847951898e+00  
+2.78487627599101e-03  -1.26306995846903e-03  -5.86288679793060e-04  
+end of step, time = 1000.000000
+1000.000000000000   -3.9565189531119833e-07  0.0000000000000000e+00
+8.04950420409227e-04  3.89700695119693e-03  1.68847450520941e-03  
+-5.73402078130268e-06  -2.47902691920044e-07  1.26091677916143e-08  
+-3.42737616773387e-01  -4.83140702164559e+00  -2.06272280903873e+00  
+7.44142091582259e-03  -6.29554642815864e-05  -2.08394209498509e-04  
+3.42042677457069e+00  7.80839540884993e+00  3.07727343840199e+00  
+-5.46603047146411e-03  1.85256191272964e-03  9.99907145008186e-04  
+-1.76077366914428e+01  -4.86639857929047e+00  -1.88189277615268e+00  
+1.08733551813521e-03  -3.61699913036238e-03  -1.59957352387901e-03  
+-1.33316007095650e+01  -2.53106572090924e+01  -1.00280661654022e+01  
+2.79920307462434e-03  -1.23621819550115e-03  -5.75654594336657e-04  
Index: proiecte/NBody/NBI/NBIserial.c
===================================================================
--- proiecte/NBody/NBI/NBIserial.c	(revision 124)
+++ proiecte/NBody/NBI/NBIserial.c	(revision 124)
@@ -0,0 +1,1922 @@
+
+
+#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>
+
+/* 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; }
+*/
+
+	
+
+
Index: proiecte/NBody/NBI/makefile
===================================================================
--- proiecte/NBody/NBI/makefile	(revision 124)
+++ proiecte/NBody/NBI/makefile	(revision 124)
@@ -0,0 +1,6 @@
+build:
+	gcc NBIserial.c -o NBIserial -lm
+	gcc -fopenmp NBIparalel.c -o NBIparalel -lm
+clean:
+	rm -f NBIserial NBIparalel
+	
Index: proiecte/NBody/NBI/scr.sh
===================================================================
--- proiecte/NBody/NBI/scr.sh	(revision 124)
+++ proiecte/NBody/NBI/scr.sh	(revision 124)
@@ -0,0 +1,6 @@
+set THREADS=$1
+set SCHED=$2
+export OMP_NUM_THREADS=$THREADS
+export OMP_SCHEDULE=$SCHED
+echo $SCHED $THREADS
+./NBIparalel NBIsampleinput
Index: proiecte/NBody/NBI/scr_serial.sh
===================================================================
--- proiecte/NBody/NBI/scr_serial.sh	(revision 124)
+++ proiecte/NBody/NBI/scr_serial.sh	(revision 124)
@@ -0,0 +1,1 @@
+~valentina.sandulescu/NBI/NBIserial ~valentina.sandulescu/NBI/NBIsampleinput
