Changeset 162


Ignore:
Timestamp:
Jan 17, 2010, 8:46:37 PM (14 years ago)
Author:
(none)
Message:
 
Location:
proiecte/NBody/NBI
Files:
2 added
3 deleted
2 edited

Legend:

Unmodified
Added
Removed

  • proiecte/NBody/NBI/NBIparalel.c

    r124 r162  
    33#define CODENAME "NBI"
    44#define VERSION    "1"
    5 
    6 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    7 
    8                         NBI
    9   A set of numerical integrators for the gravitational N-body problem.
    10 
    11         (C) Copyright 1996 Ferenc Varadi.  All rights reserved.
    12   This copyright notice is intended to prevent the commercial use of this code.
    13   The code is free and can be freely distributed.
    14   Do not remove this copyright notice unless you make changes in the code.
    15   If you make changes in the code, add your name to the copyright notice.
    16 
    17   The usual disclaimers apply, i.e., use the code at your own risk.
    18   The author does not guarantee anything.
    19 
    20   This code was developed by F. Varadi in collaboration
    21   with M. Ghil, W. I. Newman, W. M. Kaula, K. Grazier, D. Goldstein
    22   and M. Lessnick. The partial financial support of NSF Grant
    23   ATM95-23787 (M. Ghil and F. Varadi) is gratefully acknowledged.
    24 
    25   If you have any questions, suggestions etc., contact F. Varadi,
    26   e-mail: varadi@ucla.edu
    27   WWW: http://www.atmos.ucla.edu/~varadi
    28 
    29  
    30   For the impatient:
    31         Save this file as NBI.c
    32         Compile the code as
    33         cc -o NBI NBI.c -lm
    34 
    35         Download the sample input file NBIsampleinput and save it as
    36         NBIsampleinput.  It is an input file for the outer Solar System, with
    37         the planets from Jupiter through Neptune based on the Jet Propulsion
    38         Laboratory's DE403 (Digital Ephemeris 403).  This was obtained
    39         from the ftp site navigator.jpl.nasa.gov, updates can be accessed
    40         there.
    41 
    42         Run the code:
    43         NBI NBIsampleinput
    44         Take a look at the output in NBIsampleoutput and the log file
    45         NBIsampleoutput.log Change the MAXNBODIES macro, i.e., the
    46         maximum number of particles, to suit your needs.  It is set to
    47         2100 in this version.
    48 
    49 1) Description:
    50         Four numerical integrators for the gravitational (nonrelativistic)
    51         N-body problem are bundled into a single source code. The N bodies
    52         are assumed to be point masses. Particles with zero mass are called
    53         test particles; they are assumed not to influence each other and
    54         the particles with nonzero mass.
    55 
    56         The integrators are:
    57 
    58         a) A Wisdom-Holman-type mapping for hierarchical N-body systems.
    59                 Code name: HWH, maximum global order = 4
    60            This can deal with more general arrangements than the standard
    61            Wisdom-Holman mapping.  Hierarchical N-body systems are described
    62            by Roy (1988). These are represented in the code by binary trees.
    63            The integrator takes advantage of certain properties of
    64            symplectic forms and Jacobi coordinates, these will be described
    65            in a paper (in preparation). We also use singularly weighted
    66            symplectic forms, these are discussed by Varadi et al. (1995).
    67 
    68         b) A modified Cowell-Stormer integrator, with modifications by
    69            W. I. Newman and his students.
    70                 Code name: CSN, max global order = 15
    71            Note: (global order) = (local order) - 2 (Goldstein, 1996)
    72 
    73         c) A Gragg-Bulirsch-Stoer integrator, as it is described by
    74            Hairer et al., (1993).
    75                 Code name: GBS, max order = 20 (max stage = 10)
    76         IMPORTANT: The user has to specify the number of stages
    77         in the input file and not the order. Order = 2*(number of stages),
    78         see the Notes below.
    79 
    80         d) A symplectic mapping with Kinetic-Potential energy splitting
    81                 Code name: SKP, max global order = 4
    82            This is included for the sake of completeness.
    83 
    84        
    85 2) References:
    86        
    87         Goldstein, D.: 1996, The Near-Optimality of Stormer
    88         Methods for Long Time Integrations of y''=g(y),
    89         Ph.D. Dissertation, Univ. of California, Los Angeles,
    90         Dept. of Mathematics
    91 
    92         Hairer, E., Norsett, S. P. and Wanner, G.: 1993,
    93         Solving Ordinary Differential Equations I.  Nonstiff Problems.
    94         Second Revised Edition
    95        
    96         Lessnick, M.: 1996, Stability Analysis of Symplectic
    97         Integration Schemes, Ph.D. Dissertation, Univ. of California,
    98         Los Angeles, Dept. of Mathematics
    99 
    100         Roy, A. E.: 1988, Orbital Motion, Institute of Physics
    101         Publishing, Bristol
    102        
    103         Varadi, F. De la Barre, Kaula, W. M. and Ghil, M.: 1995,
    104         Singularly Weighted Symplectic Forms and Applications to
    105         Asteroid Motion, Celestial Mechanics and Dynamical Astronomy,
    106         vol. 62, pp. 23-41
    107 
    108         Yoshida, H.: 1993, Recent Progress in The Theory and
    109         Application of Symplectic Integrators,
    110         Celestial Mechanics and Dynamical Astronomy, vol. 56, pp. 27-43
    111        
    112        
    113 
    114 3) Notes
    115        
    116         All long-term integrations suffer from roundoff errors.
    117         Symplectic schemes guarantee symplecticity for short
    118         integrations but they are not strictly symplectic
    119         for very long integrations.  The implementation of
    120         the Wisdom-Holman mapping in this code was aimed at
    121         reducing the effects of roundoff errors and thus make
    122         the scheme as symplectic as possible.
    123        
    124         The Cowell-Stormer integrator, being a multistep
    125         integrator, has to be intialized.  In theory, the starting
    126         points should be computed with machine precision. This
    127         can be done using quadrupole precision arithmetic but
    128         the Gragg-Bulirsch-Stoer scheme in double precision
    129         appers to be quite adequate. It should be noted that
    130         the stability of the Cowell-Stormer scheme ensures
    131         that initial errors do not grow due to computational
    132         limitations (Goldstein, 1996).
    133 
    134         In the case of the Gragg-Bulirsch-Stoer scheme one
    135         specifies the number of stages k in the extrapolation
    136         scheme. The order of the method is then 2*k.
    137         This factor of 2 comes from the fact that the
    138         underlying method (Gragg's) has only even powers
    139         of the step size in the expansion of the local error
    140         (cf. Hairer et al., 1993, especially equation 9.22).
    141 
    142         All integrators use only G*mass, there is no need for
    143         the masses themselves nor the value of G.  Any set of units
    144         can be used:
    145         distance_unit for position,
    146         time_unit for time and step size,
    147         distance_unit/time_unit for velocity.
    148         E.g., one can use AU, day and AU/day as units.
    149 
    150         None of the codes is regularized in any sense. We intend
    151         to provide regularized code in the future. For the Wisdom-Holman
    152         mapping this involves changing the tree of Jacobian reference
    153         frames and dealing with the transition region accurately.
    154 
    155         In the case of the Wisdom-Holman mapping, the computation of
    156         acceleration differences between Jacobi and inertial coordinates
    157         could be improved in order to further reduce roundoff error.
    158         These differences can be represented by multi-pole expansions
    159         of the potential but we have not explored the details yet.
    160 
    161         The present version of the integrator does not monitor error.
    162 
    163         The choice of step size is up to the user. In most cases
    164         one has to adjust the step size to the shortest orbital
    165         period T in the system. The Wisdom-Holman mapping at second
    166         order is usually used with step size of T/10-T/100.
    167         We recommend T/1000 for Cowell-Stormer integrator
    168         since this renders the calculation EXACT to double
    169         precision, even for high (e=0.5) eccentricities.
    170         The cost of using this short time step is substantially
    171         compensated by reduced computational complexity (as
    172         compared to the Wisdom-Holman mapping).  For much larger
    173         step sizes the integrator is not stable.
    174 
    175         The Gragg-Bulirsch-Stoer method appears to be stable
    176         for large orders and step sizes. For order 20 one can
    177         even use T/5 as the step size and still get accurate results.
    178         Since this is an accurate but compuationally very expensive
    179         integrator, it is mainly used by us for short integrations.
    180         We have not investigated the propagation of roundoff
    181         error in the GBS method.
    182        
    183 
    184 4) How to compile:
    185         Since there is only a single source file, it is very easy to
    186         compile the code with whatever C compiler is available. The
    187         code was compiled and tested on Sun SPARCstations (Solaris 2.5.1
    188         and earlier) and DEC Alpha workstations (Digital UNIX V3.2).
    189         For instance, on the SPARCstations,
    190        
    191         cc -fast -xO5 -xdepend -o NBI NBI.c -lm
    192 
    193         should work with Sun's C compiler.  The code uses only the standard
    194         C library, more exactly: I/O functions, exit, malloc and free for
    195         the tree, sqrt, cbrt, sin, cos, sinh, cosh.
    196 
    197         Caution: optimization may introduce problems with some compilers
    198         and it is prudent to check the correctness of the optimized code
    199         by compiling without optimization, e.g.,
    200         cc -g -o NBI_no_op NBI.c -lm
    201         By comparing the output from NBI and NBI_no_op one can detect if
    202         the optimization leads to any problems.  (It should not but in
    203         practice the situation is not so clear.)
    204 
    205         One also can play with various compiler switches which control
    206         roundoff but we have not explored these.
    207 
    208 5) How to run
    209 
    210         NBI inputfile
    211 
    212         (Everything has to be specified in inputfile.)
    213 
    214 6) On hierarchical N-body systems
    215 
    216         We can give only a brief introduction into these; a detailed technical
    217         description is in preparation.  Our approach is somewhat different
    218         from that of Roy (1988).
    219 
    220         One starts with the concept of joining subsystems and applies this
    221         concept recursively to form new systems. When two subsystems are
    222         joined, the two subsystems are assumed to behave as point masses
    223         and thus define the appropriate two-body problem for the Wisdom-Holman
    224         mapping. The subsystems are, in general, not point masses and one
    225         takes this into account in the perturbation step of the Wisdom-Holman
    226         mapping.
    227 
    228         Instead of a simple chain of Jacobi coordinates, the code uses a
    229         tree of Jacobi coordinates. This makes it possible to use the code
    230         e.g., in situations when the massive particles have sattelites, either
    231         massive ones or not. It also reduces, when carefully specified, the
    232         integration error since more of the perturbations are absorbed into
    233         two-body interactions. The tree has to be specified in the input file. 
    234         One can define this tree formally as it is processed by the code which
    235         works like a simple finite automaton but we did not use lex, yacc or
    236         any other compiler-compiler tools.
    237 
    238  Tree description rules:
    239  Rule 1) particles are indexed by consecutive natural numbers, from 1 to N,
    240  Rule 2) particles are subsystems,
    241  Rule 3) (x y) : join the subsystems x and y to create a new subsystem.,
    242  Rule 4) [i1 - i2] denotes a set of test particles, from index i1 to index i2,
    243         which are treated the same way,
    244         Note: the space before and after - is necessary!
    245  Rule 5) ([i1 - i2] i3) is forbidden, use instead (i3 [i1 - i2])
    246 
    247  Example: 2000 main belt asteroids, from  6 to 2006, and
    248         Sun, Jupiter, Saturn, Uranus and Neptune (1 - 5)
    249         can described as
    250  (((((1 [6 - 2006]) 2) 3) 4) 5)
    251 
    252  Another example for the tree specification when test particles between
    253  Jupiter and Saturn are added with indeces from 6 to 1000 and another set
    254  of test particles is added between Saturn and Uranus with indices from
    255  1001 to 2000. The indices for the massive bodies are:
    256  Sun 1, Jupiter, 2, Saturn 3, Uranus 4, Neptune 5.
    257 ((((((1 2) [6 - 1000]) 3) [1001 - 2000]) 4) 5)
    258 
    259  Yet another example: the Sun, Earth, Moon, Jupiter and Io with indeces
    260  1,2,3,4 and 5, respectively, for which the tree is
    261  ((1 (2 3)) (4 5))
    262 
    263 
    264 7) Format of the input file
    265 Id method order step_size  total_integration_time printing_time
    266         number_of_particles
    267 
    268  tree_description
    269 
    270 GMs of particles, possibly zeros
    271 
    272 Initial conditions: x,y,z, vx, vy, vz
    273 
    274 Explanation:
    275 
    276         Id is an arbitrary string which specifies the name of the experiment.
    277 
    278         Method can be SKP, HWH, CSN and GBS.
    279 
    280         Tree_description is used only when the method is HWH or the HWH is used
    281         to initialize CSN. One can put a dummy tree description of the form
    282         (1 2)
    283         in the input file when HWH is not used.  This will satisfy the parser
    284         which processes the tree description but will not play any role in the
    285         integration.
    286 
    287 
    288 Example:
    289  For Sun, Jupiter, Saturn, Uranus and Neptune,
    290  with indeces 1,2,3,4 and 5, in heliocentric coordinates.
    291 
    292 SomeId HWH 4 100.0 365.0e+4 1.0e+4
    293 5
    294 ((((1 2) 3) 4) 5)
    295 
    296 0.295912208285591095e-03   
    297 0.282534590952422643e-06  0.845971518568065874e-07 
    298 0.129202491678196939e-07   0.152435890078427628e-07 
    299 
    300 0.0 0.0 0.0 0.0 0.0 0.0
    301 
    302 -0.538420940678015203e+01   -0.831247035699103631e+00 -0.225095848657298592e+00 
    303 0.109236311953937086e-02   -0.652329334256263223e-02 -0.282301443780311710e-02 
    304 
    305 0.788989169798583756e+01   0.459570785756998301e+01  0.155843220515231762e+01
    306 -0.321720261683698921e-02   0.433063255278440598e-02  0.192641782724487106e-02
    307 
    308 -0.182698980946940850e+02   -0.116271406901560681e+01 -0.250371938307287545e+00 
    309 0.221540447608191936e-03   -0.376765389967669917e-02 -0.165324449144407023e-02
    310 
    311 -0.160595376822070470e+02   -0.239429495784517528e+02 -0.940042953888094424e+01 
    312 0.264312302715172384e-02   -0.150349219713705076e-02 -0.681271071793930938e-03
    313 
    314 
    315 There are 2 output files:
    316         1) SomeId contains the actual integration data
    317         2) SomeId.log contains a copy of the input data and error messages.
    318         Some error messages might be sent to the standard output and the user
    319         has to redirect it in order to capture them.
    320 
    321 In the case of the Cowell-Stormer integrator, an additional file, SomeId.lerr
    322 is created. This contains the LOCAL truncation error of the integration,
    323 computed for each step but only printed for the last step in the printing
    324 cycle. This also should help to determine the appropriate step size.
    325 
    326 The output is in the original coordinate system, i.e., inertial and not Jacobian.
    327 Output is produced when time is an integer multiple of the specified printing
    328 time. The following is printed:
    329 time    relative_change_in_total_energy relative_change_in_ang.mom.z.comp.
    330 after which
    331 x y z
    332 dx dy dz
    333 for each particle.
    334 
    335 
    336 9) The uset can change the code, of course. The section MAIN LOOP is where
    337    this can be done easily. You can insert code to do more things, e.g.,
    338    detecting close approaches.
    339 
    340 10) MACROS TO GIVE SOME CONTROL OVER THINGS.
    341 Changes in the code that you can make HERE:
    342         CHANGE the DIM macro to 2 for planar motions.
    343         (Note: a variables dim is also used, easy to change
    344                 it everywhere to DIM. It is still better to keep
    345                 the Kepler solver in the more general 3-dimensional form.)
    346         CHANGE the MAXNBODIES macro to whatever you want,
    347                 including test particles.
    348 
    349 **************** END OF DESCRIPTION ********************** */
    350 
    351 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    352         USER-CHANGEABLE MACROS    */
    3535
    3546/* dimension of physical space */
    … …  
    38335
    38436
    385 /* VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV */
    386 
    387 /*      The actual code.
    388         Change things below this line at your own peril.
    389        
    390  */
    391 
    39237/* MACROS TO DISTINGUISH BETWEEN TEST PARTICLES AND
    39338   PARTICLES WITH FINITE MASS
    394    You can add more types here, these values are stored
    395         in tpstore for each particle
    396    e.g., you could define stopped test particles,
    397    but you have to change the code accordingly.
    398    It is perhaps easier to add more integer arrays.
    399    IMPORTANT: THE TREE MAINTAINS ITS OWN INTEGER ARRAYS
    400    WHICH HOLD TEST PARTICLE INDICES.
    40139*/
    40240

  • proiecte/NBody/NBI/NBIserial.c

    r124 r162  
    33#define CODENAME "NBI"
    44#define VERSION    "1"
    5 
    6 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    7 
    8                         NBI
    9   A set of numerical integrators for the gravitational N-body problem.
    10 
    11         (C) Copyright 1996 Ferenc Varadi.  All rights reserved.
    12   This copyright notice is intended to prevent the commercial use of this code.
    13   The code is free and can be freely distributed.
    14   Do not remove this copyright notice unless you make changes in the code.
    15   If you make changes in the code, add your name to the copyright notice.
    16 
    17   The usual disclaimers apply, i.e., use the code at your own risk.
    18   The author does not guarantee anything.
    19 
    20   This code was developed by F. Varadi in collaboration
    21   with M. Ghil, W. I. Newman, W. M. Kaula, K. Grazier, D. Goldstein
    22   and M. Lessnick. The partial financial support of NSF Grant
    23   ATM95-23787 (M. Ghil and F. Varadi) is gratefully acknowledged.
    24 
    25   If you have any questions, suggestions etc., contact F. Varadi,
    26   e-mail: varadi@ucla.edu
    27   WWW: http://www.atmos.ucla.edu/~varadi
    28 
    29  
    30   For the impatient:
    31         Save this file as NBI.c
    32         Compile the code as
    33         cc -o NBI NBI.c -lm
    34 
    35         Download the sample input file NBIsampleinput and save it as
    36         NBIsampleinput.  It is an input file for the outer Solar System, with
    37         the planets from Jupiter through Neptune based on the Jet Propulsion
    38         Laboratory's DE403 (Digital Ephemeris 403).  This was obtained
    39         from the ftp site navigator.jpl.nasa.gov, updates can be accessed
    40         there.
    41 
    42         Run the code:
    43         NBI NBIsampleinput
    44         Take a look at the output in NBIsampleoutput and the log file
    45         NBIsampleoutput.log Change the MAXNBODIES macro, i.e., the
    46         maximum number of particles, to suit your needs.  It is set to
    47         2100 in this version.
    48 
    49 1) Description:
    50         Four numerical integrators for the gravitational (nonrelativistic)
    51         N-body problem are bundled into a single source code. The N bodies
    52         are assumed to be point masses. Particles with zero mass are called
    53         test particles; they are assumed not to influence each other and
    54         the particles with nonzero mass.
    55 
    56         The integrators are:
    57 
    58         a) A Wisdom-Holman-type mapping for hierarchical N-body systems.
    59                 Code name: HWH, maximum global order = 4
    60            This can deal with more general arrangements than the standard
    61            Wisdom-Holman mapping.  Hierarchical N-body systems are described
    62            by Roy (1988). These are represented in the code by binary trees.
    63            The integrator takes advantage of certain properties of
    64            symplectic forms and Jacobi coordinates, these will be described
    65            in a paper (in preparation). We also use singularly weighted
    66            symplectic forms, these are discussed by Varadi et al. (1995).
    67 
    68         b) A modified Cowell-Stormer integrator, with modifications by
    69            W. I. Newman and his students.
    70                 Code name: CSN, max global order = 15
    71            Note: (global order) = (local order) - 2 (Goldstein, 1996)
    72 
    73         c) A Gragg-Bulirsch-Stoer integrator, as it is described by
    74            Hairer et al., (1993).
    75                 Code name: GBS, max order = 20 (max stage = 10)
    76         IMPORTANT: The user has to specify the number of stages
    77         in the input file and not the order. Order = 2*(number of stages),
    78         see the Notes below.
    79 
    80         d) A symplectic mapping with Kinetic-Potential energy splitting
    81                 Code name: SKP, max global order = 4
    82            This is included for the sake of completeness.
    83 
    84        
    85 2) References:
    86        
    87         Goldstein, D.: 1996, The Near-Optimality of Stormer
    88         Methods for Long Time Integrations of y''=g(y),
    89         Ph.D. Dissertation, Univ. of California, Los Angeles,
    90         Dept. of Mathematics
    91 
    92         Hairer, E., Norsett, S. P. and Wanner, G.: 1993,
    93         Solving Ordinary Differential Equations I.  Nonstiff Problems.
    94         Second Revised Edition
    95        
    96         Lessnick, M.: 1996, Stability Analysis of Symplectic
    97         Integration Schemes, Ph.D. Dissertation, Univ. of California,
    98         Los Angeles, Dept. of Mathematics
    99 
    100         Roy, A. E.: 1988, Orbital Motion, Institute of Physics
    101         Publishing, Bristol
    102        
    103         Varadi, F. De la Barre, Kaula, W. M. and Ghil, M.: 1995,
    104         Singularly Weighted Symplectic Forms and Applications to
    105         Asteroid Motion, Celestial Mechanics and Dynamical Astronomy,
    106         vol. 62, pp. 23-41
    107 
    108         Yoshida, H.: 1993, Recent Progress in The Theory and
    109         Application of Symplectic Integrators,
    110         Celestial Mechanics and Dynamical Astronomy, vol. 56, pp. 27-43
    111        
    112        
    113 
    114 3) Notes
    115        
    116         All long-term integrations suffer from roundoff errors.
    117         Symplectic schemes guarantee symplecticity for short
    118         integrations but they are not strictly symplectic
    119         for very long integrations.  The implementation of
    120         the Wisdom-Holman mapping in this code was aimed at
    121         reducing the effects of roundoff errors and thus make
    122         the scheme as symplectic as possible.
    123        
    124         The Cowell-Stormer integrator, being a multistep
    125         integrator, has to be intialized.  In theory, the starting
    126         points should be computed with machine precision. This
    127         can be done using quadrupole precision arithmetic but
    128         the Gragg-Bulirsch-Stoer scheme in double precision
    129         appers to be quite adequate. It should be noted that
    130         the stability of the Cowell-Stormer scheme ensures
    131         that initial errors do not grow due to computational
    132         limitations (Goldstein, 1996).
    133 
    134         In the case of the Gragg-Bulirsch-Stoer scheme one
    135         specifies the number of stages k in the extrapolation
    136         scheme. The order of the method is then 2*k.
    137         This factor of 2 comes from the fact that the
    138         underlying method (Gragg's) has only even powers
    139         of the step size in the expansion of the local error
    140         (cf. Hairer et al., 1993, especially equation 9.22).
    141 
    142         All integrators use only G*mass, there is no need for
    143         the masses themselves nor the value of G.  Any set of units
    144         can be used:
    145         distance_unit for position,
    146         time_unit for time and step size,
    147         distance_unit/time_unit for velocity.
    148         E.g., one can use AU, day and AU/day as units.
    149 
    150         None of the codes is regularized in any sense. We intend
    151         to provide regularized code in the future. For the Wisdom-Holman
    152         mapping this involves changing the tree of Jacobian reference
    153         frames and dealing with the transition region accurately.
    154 
    155         In the case of the Wisdom-Holman mapping, the computation of
    156         acceleration differences between Jacobi and inertial coordinates
    157         could be improved in order to further reduce roundoff error.
    158         These differences can be represented by multi-pole expansions
    159         of the potential but we have not explored the details yet.
    160 
    161         The present version of the integrator does not monitor error.
    162 
    163         The choice of step size is up to the user. In most cases
    164         one has to adjust the step size to the shortest orbital
    165         period T in the system. The Wisdom-Holman mapping at second
    166         order is usually used with step size of T/10-T/100.
    167         We recommend T/1000 for Cowell-Stormer integrator
    168         since this renders the calculation EXACT to double
    169         precision, even for high (e=0.5) eccentricities.
    170         The cost of using this short time step is substantially
    171         compensated by reduced computational complexity (as
    172         compared to the Wisdom-Holman mapping).  For much larger
    173         step sizes the integrator is not stable.
    174 
    175         The Gragg-Bulirsch-Stoer method appears to be stable
    176         for large orders and step sizes. For order 20 one can
    177         even use T/5 as the step size and still get accurate results.
    178         Since this is an accurate but compuationally very expensive
    179         integrator, it is mainly used by us for short integrations.
    180         We have not investigated the propagation of roundoff
    181         error in the GBS method.
    182        
    183 
    184 4) How to compile:
    185         Since there is only a single source file, it is very easy to
    186         compile the code with whatever C compiler is available. The
    187         code was compiled and tested on Sun SPARCstations (Solaris 2.5.1
    188         and earlier) and DEC Alpha workstations (Digital UNIX V3.2).
    189         For instance, on the SPARCstations,
    190        
    191         cc -fast -xO5 -xdepend -o NBI NBI.c -lm
    192 
    193         should work with Sun's C compiler.  The code uses only the standard
    194         C library, more exactly: I/O functions, exit, malloc and free for
    195         the tree, sqrt, cbrt, sin, cos, sinh, cosh.
    196 
    197         Caution: optimization may introduce problems with some compilers
    198         and it is prudent to check the correctness of the optimized code
    199         by compiling without optimization, e.g.,
    200         cc -g -o NBI_no_op NBI.c -lm
    201         By comparing the output from NBI and NBI_no_op one can detect if
    202         the optimization leads to any problems.  (It should not but in
    203         practice the situation is not so clear.)
    204 
    205         One also can play with various compiler switches which control
    206         roundoff but we have not explored these.
    207 
    208 5) How to run
    209 
    210         NBI inputfile
    211 
    212         (Everything has to be specified in inputfile.)
    213 
    214 6) On hierarchical N-body systems
    215 
    216         We can give only a brief introduction into these; a detailed technical
    217         description is in preparation.  Our approach is somewhat different
    218         from that of Roy (1988).
    219 
    220         One starts with the concept of joining subsystems and applies this
    221         concept recursively to form new systems. When two subsystems are
    222         joined, the two subsystems are assumed to behave as point masses
    223         and thus define the appropriate two-body problem for the Wisdom-Holman
    224         mapping. The subsystems are, in general, not point masses and one
    225         takes this into account in the perturbation step of the Wisdom-Holman
    226         mapping.
    227 
    228         Instead of a simple chain of Jacobi coordinates, the code uses a
    229         tree of Jacobi coordinates. This makes it possible to use the code
    230         e.g., in situations when the massive particles have sattelites, either
    231         massive ones or not. It also reduces, when carefully specified, the
    232         integration error since more of the perturbations are absorbed into
    233         two-body interactions. The tree has to be specified in the input file. 
    234         One can define this tree formally as it is processed by the code which
    235         works like a simple finite automaton but we did not use lex, yacc or
    236         any other compiler-compiler tools.
    237 
    238  Tree description rules:
    239  Rule 1) particles are indexed by consecutive natural numbers, from 1 to N,
    240  Rule 2) particles are subsystems,
    241  Rule 3) (x y) : join the subsystems x and y to create a new subsystem.,
    242  Rule 4) [i1 - i2] denotes a set of test particles, from index i1 to index i2,
    243         which are treated the same way,
    244         Note: the space before and after - is necessary!
    245  Rule 5) ([i1 - i2] i3) is forbidden, use instead (i3 [i1 - i2])
    246 
    247  Example: 2000 main belt asteroids, from  6 to 2006, and
    248         Sun, Jupiter, Saturn, Uranus and Neptune (1 - 5)
    249         can described as
    250  (((((1 [6 - 2006]) 2) 3) 4) 5)
    251 
    252  Another example for the tree specification when test particles between
    253  Jupiter and Saturn are added with indeces from 6 to 1000 and another set
    254  of test particles is added between Saturn and Uranus with indices from
    255  1001 to 2000. The indices for the massive bodies are:
    256  Sun 1, Jupiter, 2, Saturn 3, Uranus 4, Neptune 5.
    257 ((((((1 2) [6 - 1000]) 3) [1001 - 2000]) 4) 5)
    258 
    259  Yet another example: the Sun, Earth, Moon, Jupiter and Io with indeces
    260  1,2,3,4 and 5, respectively, for which the tree is
    261  ((1 (2 3)) (4 5))
    262 
    263 
    264 7) Format of the input file
    265 Id method order step_size  total_integration_time printing_time
    266         number_of_particles
    267 
    268  tree_description
    269 
    270 GMs of particles, possibly zeros
    271 
    272 Initial conditions: x,y,z, vx, vy, vz
    273 
    274 Explanation:
    275 
    276         Id is an arbitrary string which specifies the name of the experiment.
    277 
    278         Method can be SKP, HWH, CSN and GBS.
    279 
    280         Tree_description is used only when the method is HWH or the HWH is used
    281         to initialize CSN. One can put a dummy tree description of the form
    282         (1 2)
    283         in the input file when HWH is not used.  This will satisfy the parser
    284         which processes the tree description but will not play any role in the
    285         integration.
    286 
    287 
    288 Example:
    289  For Sun, Jupiter, Saturn, Uranus and Neptune,
    290  with indeces 1,2,3,4 and 5, in heliocentric coordinates.
    291 
    292 SomeId HWH 4 100.0 365.0e+4 1.0e+4
    293 5
    294 ((((1 2) 3) 4) 5)
    295 
    296 0.295912208285591095e-03   
    297 0.282534590952422643e-06  0.845971518568065874e-07 
    298 0.129202491678196939e-07   0.152435890078427628e-07 
    299 
    300 0.0 0.0 0.0 0.0 0.0 0.0
    301 
    302 -0.538420940678015203e+01   -0.831247035699103631e+00 -0.225095848657298592e+00 
    303 0.109236311953937086e-02   -0.652329334256263223e-02 -0.282301443780311710e-02 
    304 
    305 0.788989169798583756e+01   0.459570785756998301e+01  0.155843220515231762e+01
    306 -0.321720261683698921e-02   0.433063255278440598e-02  0.192641782724487106e-02
    307 
    308 -0.182698980946940850e+02   -0.116271406901560681e+01 -0.250371938307287545e+00 
    309 0.221540447608191936e-03   -0.376765389967669917e-02 -0.165324449144407023e-02
    310 
    311 -0.160595376822070470e+02   -0.239429495784517528e+02 -0.940042953888094424e+01 
    312 0.264312302715172384e-02   -0.150349219713705076e-02 -0.681271071793930938e-03
    313 
    314 
    315 There are 2 output files:
    316         1) SomeId contains the actual integration data
    317         2) SomeId.log contains a copy of the input data and error messages.
    318         Some error messages might be sent to the standard output and the user
    319         has to redirect it in order to capture them.
    320 
    321 In the case of the Cowell-Stormer integrator, an additional file, SomeId.lerr
    322 is created. This contains the LOCAL truncation error of the integration,
    323 computed for each step but only printed for the last step in the printing
    324 cycle. This also should help to determine the appropriate step size.
    325 
    326 The output is in the original coordinate system, i.e., inertial and not Jacobian.
    327 Output is produced when time is an integer multiple of the specified printing
    328 time. The following is printed:
    329 time    relative_change_in_total_energy relative_change_in_ang.mom.z.comp.
    330 after which
    331 x y z
    332 dx dy dz
    333 for each particle.
    334 
    335 
    336 9) The uset can change the code, of course. The section MAIN LOOP is where
    337    this can be done easily. You can insert code to do more things, e.g.,
    338    detecting close approaches.
    339 
    340 10) MACROS TO GIVE SOME CONTROL OVER THINGS.
    341 Changes in the code that you can make HERE:
    342         CHANGE the DIM macro to 2 for planar motions.
    343         (Note: a variables dim is also used, easy to change
    344                 it everywhere to DIM. It is still better to keep
    345                 the Kepler solver in the more general 3-dimensional form.)
    346         CHANGE the MAXNBODIES macro to whatever you want,
    347                 including test particles.
    348 
    349 **************** END OF DESCRIPTION ********************** */
    350 
    351 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    352         USER-CHANGEABLE MACROS    */
    3535
    3546/* dimension of physical space */
    … …  
    38133/* cange this line accordingly */
    38234EXSEQBUL
    383 
    384 
    385 /* VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV */
    386 
    387 /*      The actual code.
    388         Change things below this line at your own peril.
    389        
    390  */
    391 
    392 /* MACROS TO DISTINGUISH BETWEEN TEST PARTICLES AND
    393    PARTICLES WITH FINITE MASS
    394    You can add more types here, these values are stored
    395         in tpstore for each particle
    396    e.g., you could define stopped test particles,
    397    but you have to change the code accordingly.
    398    It is perhaps easier to add more integer arrays.
    399    IMPORTANT: THE TREE MAINTAINS ITS OWN INTEGER ARRAYS
    400    WHICH HOLD TEST PARTICLE INDICES.
    401 */
    40235
    40336#define NOTTESTPARTICLE         1
Note: See TracChangeset for help on using the changeset viewer.