Index: proiecte/NBody/NBI/NBIparalel.c
===================================================================
--- proiecte/NBody/NBI/NBIparalel.c	(revision 124)
+++ proiecte/NBody/NBI/NBIparalel.c	(revision 162)
@@ -3,352 +3,4 @@
 #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 */
@@ -383,20 +35,6 @@
 
 
-/* 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. 
 */
 
