- Timestamp:
- Jan 17, 2010, 8:46:37 PM (14 years ago)
- Location:
- proiecte/NBody/NBI
- Files:
-
- 2 added
- 3 deleted
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
proiecte/NBody/NBI/NBIparalel.c
r124 r162 3 3 #define CODENAME "NBI" 4 4 #define VERSION "1" 5 6 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++7 8 NBI9 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 collaboration21 with M. Ghil, W. I. Newman, W. M. Kaula, K. Grazier, D. Goldstein22 and M. Lessnick. The partial financial support of NSF Grant23 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.edu27 WWW: http://www.atmos.ucla.edu/~varadi28 29 30 For the impatient:31 Save this file as NBI.c32 Compile the code as33 cc -o NBI NBI.c -lm34 35 Download the sample input file NBIsampleinput and save it as36 NBIsampleinput. It is an input file for the outer Solar System, with37 the planets from Jupiter through Neptune based on the Jet Propulsion38 Laboratory's DE403 (Digital Ephemeris 403). This was obtained39 from the ftp site navigator.jpl.nasa.gov, updates can be accessed40 there.41 42 Run the code:43 NBI NBIsampleinput44 Take a look at the output in NBIsampleoutput and the log file45 NBIsampleoutput.log Change the MAXNBODIES macro, i.e., the46 maximum number of particles, to suit your needs. It is set to47 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 bodies52 are assumed to be point masses. Particles with zero mass are called53 test particles; they are assumed not to influence each other and54 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 = 460 This can deal with more general arrangements than the standard61 Wisdom-Holman mapping. Hierarchical N-body systems are described62 by Roy (1988). These are represented in the code by binary trees.63 The integrator takes advantage of certain properties of64 symplectic forms and Jacobi coordinates, these will be described65 in a paper (in preparation). We also use singularly weighted66 symplectic forms, these are discussed by Varadi et al. (1995).67 68 b) A modified Cowell-Stormer integrator, with modifications by69 W. I. Newman and his students.70 Code name: CSN, max global order = 1571 Note: (global order) = (local order) - 2 (Goldstein, 1996)72 73 c) A Gragg-Bulirsch-Stoer integrator, as it is described by74 Hairer et al., (1993).75 Code name: GBS, max order = 20 (max stage = 10)76 IMPORTANT: The user has to specify the number of stages77 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 splitting81 Code name: SKP, max global order = 482 This is included for the sake of completeness.83 84 85 2) References:86 87 Goldstein, D.: 1996, The Near-Optimality of Stormer88 Methods for Long Time Integrations of y''=g(y),89 Ph.D. Dissertation, Univ. of California, Los Angeles,90 Dept. of Mathematics91 92 Hairer, E., Norsett, S. P. and Wanner, G.: 1993,93 Solving Ordinary Differential Equations I. Nonstiff Problems.94 Second Revised Edition95 96 Lessnick, M.: 1996, Stability Analysis of Symplectic97 Integration Schemes, Ph.D. Dissertation, Univ. of California,98 Los Angeles, Dept. of Mathematics99 100 Roy, A. E.: 1988, Orbital Motion, Institute of Physics101 Publishing, Bristol102 103 Varadi, F. De la Barre, Kaula, W. M. and Ghil, M.: 1995,104 Singularly Weighted Symplectic Forms and Applications to105 Asteroid Motion, Celestial Mechanics and Dynamical Astronomy,106 vol. 62, pp. 23-41107 108 Yoshida, H.: 1993, Recent Progress in The Theory and109 Application of Symplectic Integrators,110 Celestial Mechanics and Dynamical Astronomy, vol. 56, pp. 27-43111 112 113 114 3) Notes115 116 All long-term integrations suffer from roundoff errors.117 Symplectic schemes guarantee symplecticity for short118 integrations but they are not strictly symplectic119 for very long integrations. The implementation of120 the Wisdom-Holman mapping in this code was aimed at121 reducing the effects of roundoff errors and thus make122 the scheme as symplectic as possible.123 124 The Cowell-Stormer integrator, being a multistep125 integrator, has to be intialized. In theory, the starting126 points should be computed with machine precision. This127 can be done using quadrupole precision arithmetic but128 the Gragg-Bulirsch-Stoer scheme in double precision129 appers to be quite adequate. It should be noted that130 the stability of the Cowell-Stormer scheme ensures131 that initial errors do not grow due to computational132 limitations (Goldstein, 1996).133 134 In the case of the Gragg-Bulirsch-Stoer scheme one135 specifies the number of stages k in the extrapolation136 scheme. The order of the method is then 2*k.137 This factor of 2 comes from the fact that the138 underlying method (Gragg's) has only even powers139 of the step size in the expansion of the local error140 (cf. Hairer et al., 1993, especially equation 9.22).141 142 All integrators use only G*mass, there is no need for143 the masses themselves nor the value of G. Any set of units144 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 intend151 to provide regularized code in the future. For the Wisdom-Holman152 mapping this involves changing the tree of Jacobian reference153 frames and dealing with the transition region accurately.154 155 In the case of the Wisdom-Holman mapping, the computation of156 acceleration differences between Jacobi and inertial coordinates157 could be improved in order to further reduce roundoff error.158 These differences can be represented by multi-pole expansions159 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 cases164 one has to adjust the step size to the shortest orbital165 period T in the system. The Wisdom-Holman mapping at second166 order is usually used with step size of T/10-T/100.167 We recommend T/1000 for Cowell-Stormer integrator168 since this renders the calculation EXACT to double169 precision, even for high (e=0.5) eccentricities.170 The cost of using this short time step is substantially171 compensated by reduced computational complexity (as172 compared to the Wisdom-Holman mapping). For much larger173 step sizes the integrator is not stable.174 175 The Gragg-Bulirsch-Stoer method appears to be stable176 for large orders and step sizes. For order 20 one can177 even use T/5 as the step size and still get accurate results.178 Since this is an accurate but compuationally very expensive179 integrator, it is mainly used by us for short integrations.180 We have not investigated the propagation of roundoff181 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 to186 compile the code with whatever C compiler is available. The187 code was compiled and tested on Sun SPARCstations (Solaris 2.5.1188 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 -lm192 193 should work with Sun's C compiler. The code uses only the standard194 C library, more exactly: I/O functions, exit, malloc and free for195 the tree, sqrt, cbrt, sin, cos, sinh, cosh.196 197 Caution: optimization may introduce problems with some compilers198 and it is prudent to check the correctness of the optimized code199 by compiling without optimization, e.g.,200 cc -g -o NBI_no_op NBI.c -lm201 By comparing the output from NBI and NBI_no_op one can detect if202 the optimization leads to any problems. (It should not but in203 practice the situation is not so clear.)204 205 One also can play with various compiler switches which control206 roundoff but we have not explored these.207 208 5) How to run209 210 NBI inputfile211 212 (Everything has to be specified in inputfile.)213 214 6) On hierarchical N-body systems215 216 We can give only a brief introduction into these; a detailed technical217 description is in preparation. Our approach is somewhat different218 from that of Roy (1988).219 220 One starts with the concept of joining subsystems and applies this221 concept recursively to form new systems. When two subsystems are222 joined, the two subsystems are assumed to behave as point masses223 and thus define the appropriate two-body problem for the Wisdom-Holman224 mapping. The subsystems are, in general, not point masses and one225 takes this into account in the perturbation step of the Wisdom-Holman226 mapping.227 228 Instead of a simple chain of Jacobi coordinates, the code uses a229 tree of Jacobi coordinates. This makes it possible to use the code230 e.g., in situations when the massive particles have sattelites, either231 massive ones or not. It also reduces, when carefully specified, the232 integration error since more of the perturbations are absorbed into233 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 which235 works like a simple finite automaton but we did not use lex, yacc or236 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, and248 Sun, Jupiter, Saturn, Uranus and Neptune (1 - 5)249 can described as250 (((((1 [6 - 2006]) 2) 3) 4) 5)251 252 Another example for the tree specification when test particles between253 Jupiter and Saturn are added with indeces from 6 to 1000 and another set254 of test particles is added between Saturn and Uranus with indices from255 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 indeces260 1,2,3,4 and 5, respectively, for which the tree is261 ((1 (2 3)) (4 5))262 263 264 7) Format of the input file265 Id method order step_size total_integration_time printing_time266 number_of_particles267 268 tree_description269 270 GMs of particles, possibly zeros271 272 Initial conditions: x,y,z, vx, vy, vz273 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 used281 to initialize CSN. One can put a dummy tree description of the form282 (1 2)283 in the input file when HWH is not used. This will satisfy the parser284 which processes the tree description but will not play any role in the285 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+4293 5294 ((((1 2) 3) 4) 5)295 296 0.295912208285591095e-03297 0.282534590952422643e-06 0.845971518568065874e-07298 0.129202491678196939e-07 0.152435890078427628e-07299 300 0.0 0.0 0.0 0.0 0.0 0.0301 302 -0.538420940678015203e+01 -0.831247035699103631e+00 -0.225095848657298592e+00303 0.109236311953937086e-02 -0.652329334256263223e-02 -0.282301443780311710e-02304 305 0.788989169798583756e+01 0.459570785756998301e+01 0.155843220515231762e+01306 -0.321720261683698921e-02 0.433063255278440598e-02 0.192641782724487106e-02307 308 -0.182698980946940850e+02 -0.116271406901560681e+01 -0.250371938307287545e+00309 0.221540447608191936e-03 -0.376765389967669917e-02 -0.165324449144407023e-02310 311 -0.160595376822070470e+02 -0.239429495784517528e+02 -0.940042953888094424e+01312 0.264312302715172384e-02 -0.150349219713705076e-02 -0.681271071793930938e-03313 314 315 There are 2 output files:316 1) SomeId contains the actual integration data317 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 user319 has to redirect it in order to capture them.320 321 In the case of the Cowell-Stormer integrator, an additional file, SomeId.lerr322 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 printing324 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 printing328 time. The following is printed:329 time relative_change_in_total_energy relative_change_in_ang.mom.z.comp.330 after which331 x y z332 dx dy dz333 for each particle.334 335 336 9) The uset can change the code, of course. The section MAIN LOOP is where337 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 change344 it everywhere to DIM. It is still better to keep345 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 */353 5 354 6 /* dimension of physical space */ … … 383 35 384 36 385 /* VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV */386 387 /* The actual code.388 Change things below this line at your own peril.389 390 */391 392 37 /* MACROS TO DISTINGUISH BETWEEN TEST PARTICLES AND 393 38 PARTICLES WITH FINITE MASS 394 You can add more types here, these values are stored395 in tpstore for each particle396 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 ARRAYS400 WHICH HOLD TEST PARTICLE INDICES.401 39 */ 402 40 -
proiecte/NBody/NBI/NBIserial.c
r124 r162 3 3 #define CODENAME "NBI" 4 4 #define VERSION "1" 5 6 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++7 8 NBI9 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 collaboration21 with M. Ghil, W. I. Newman, W. M. Kaula, K. Grazier, D. Goldstein22 and M. Lessnick. The partial financial support of NSF Grant23 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.edu27 WWW: http://www.atmos.ucla.edu/~varadi28 29 30 For the impatient:31 Save this file as NBI.c32 Compile the code as33 cc -o NBI NBI.c -lm34 35 Download the sample input file NBIsampleinput and save it as36 NBIsampleinput. It is an input file for the outer Solar System, with37 the planets from Jupiter through Neptune based on the Jet Propulsion38 Laboratory's DE403 (Digital Ephemeris 403). This was obtained39 from the ftp site navigator.jpl.nasa.gov, updates can be accessed40 there.41 42 Run the code:43 NBI NBIsampleinput44 Take a look at the output in NBIsampleoutput and the log file45 NBIsampleoutput.log Change the MAXNBODIES macro, i.e., the46 maximum number of particles, to suit your needs. It is set to47 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 bodies52 are assumed to be point masses. Particles with zero mass are called53 test particles; they are assumed not to influence each other and54 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 = 460 This can deal with more general arrangements than the standard61 Wisdom-Holman mapping. Hierarchical N-body systems are described62 by Roy (1988). These are represented in the code by binary trees.63 The integrator takes advantage of certain properties of64 symplectic forms and Jacobi coordinates, these will be described65 in a paper (in preparation). We also use singularly weighted66 symplectic forms, these are discussed by Varadi et al. (1995).67 68 b) A modified Cowell-Stormer integrator, with modifications by69 W. I. Newman and his students.70 Code name: CSN, max global order = 1571 Note: (global order) = (local order) - 2 (Goldstein, 1996)72 73 c) A Gragg-Bulirsch-Stoer integrator, as it is described by74 Hairer et al., (1993).75 Code name: GBS, max order = 20 (max stage = 10)76 IMPORTANT: The user has to specify the number of stages77 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 splitting81 Code name: SKP, max global order = 482 This is included for the sake of completeness.83 84 85 2) References:86 87 Goldstein, D.: 1996, The Near-Optimality of Stormer88 Methods for Long Time Integrations of y''=g(y),89 Ph.D. Dissertation, Univ. of California, Los Angeles,90 Dept. of Mathematics91 92 Hairer, E., Norsett, S. P. and Wanner, G.: 1993,93 Solving Ordinary Differential Equations I. Nonstiff Problems.94 Second Revised Edition95 96 Lessnick, M.: 1996, Stability Analysis of Symplectic97 Integration Schemes, Ph.D. Dissertation, Univ. of California,98 Los Angeles, Dept. of Mathematics99 100 Roy, A. E.: 1988, Orbital Motion, Institute of Physics101 Publishing, Bristol102 103 Varadi, F. De la Barre, Kaula, W. M. and Ghil, M.: 1995,104 Singularly Weighted Symplectic Forms and Applications to105 Asteroid Motion, Celestial Mechanics and Dynamical Astronomy,106 vol. 62, pp. 23-41107 108 Yoshida, H.: 1993, Recent Progress in The Theory and109 Application of Symplectic Integrators,110 Celestial Mechanics and Dynamical Astronomy, vol. 56, pp. 27-43111 112 113 114 3) Notes115 116 All long-term integrations suffer from roundoff errors.117 Symplectic schemes guarantee symplecticity for short118 integrations but they are not strictly symplectic119 for very long integrations. The implementation of120 the Wisdom-Holman mapping in this code was aimed at121 reducing the effects of roundoff errors and thus make122 the scheme as symplectic as possible.123 124 The Cowell-Stormer integrator, being a multistep125 integrator, has to be intialized. In theory, the starting126 points should be computed with machine precision. This127 can be done using quadrupole precision arithmetic but128 the Gragg-Bulirsch-Stoer scheme in double precision129 appers to be quite adequate. It should be noted that130 the stability of the Cowell-Stormer scheme ensures131 that initial errors do not grow due to computational132 limitations (Goldstein, 1996).133 134 In the case of the Gragg-Bulirsch-Stoer scheme one135 specifies the number of stages k in the extrapolation136 scheme. The order of the method is then 2*k.137 This factor of 2 comes from the fact that the138 underlying method (Gragg's) has only even powers139 of the step size in the expansion of the local error140 (cf. Hairer et al., 1993, especially equation 9.22).141 142 All integrators use only G*mass, there is no need for143 the masses themselves nor the value of G. Any set of units144 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 intend151 to provide regularized code in the future. For the Wisdom-Holman152 mapping this involves changing the tree of Jacobian reference153 frames and dealing with the transition region accurately.154 155 In the case of the Wisdom-Holman mapping, the computation of156 acceleration differences between Jacobi and inertial coordinates157 could be improved in order to further reduce roundoff error.158 These differences can be represented by multi-pole expansions159 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 cases164 one has to adjust the step size to the shortest orbital165 period T in the system. The Wisdom-Holman mapping at second166 order is usually used with step size of T/10-T/100.167 We recommend T/1000 for Cowell-Stormer integrator168 since this renders the calculation EXACT to double169 precision, even for high (e=0.5) eccentricities.170 The cost of using this short time step is substantially171 compensated by reduced computational complexity (as172 compared to the Wisdom-Holman mapping). For much larger173 step sizes the integrator is not stable.174 175 The Gragg-Bulirsch-Stoer method appears to be stable176 for large orders and step sizes. For order 20 one can177 even use T/5 as the step size and still get accurate results.178 Since this is an accurate but compuationally very expensive179 integrator, it is mainly used by us for short integrations.180 We have not investigated the propagation of roundoff181 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 to186 compile the code with whatever C compiler is available. The187 code was compiled and tested on Sun SPARCstations (Solaris 2.5.1188 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 -lm192 193 should work with Sun's C compiler. The code uses only the standard194 C library, more exactly: I/O functions, exit, malloc and free for195 the tree, sqrt, cbrt, sin, cos, sinh, cosh.196 197 Caution: optimization may introduce problems with some compilers198 and it is prudent to check the correctness of the optimized code199 by compiling without optimization, e.g.,200 cc -g -o NBI_no_op NBI.c -lm201 By comparing the output from NBI and NBI_no_op one can detect if202 the optimization leads to any problems. (It should not but in203 practice the situation is not so clear.)204 205 One also can play with various compiler switches which control206 roundoff but we have not explored these.207 208 5) How to run209 210 NBI inputfile211 212 (Everything has to be specified in inputfile.)213 214 6) On hierarchical N-body systems215 216 We can give only a brief introduction into these; a detailed technical217 description is in preparation. Our approach is somewhat different218 from that of Roy (1988).219 220 One starts with the concept of joining subsystems and applies this221 concept recursively to form new systems. When two subsystems are222 joined, the two subsystems are assumed to behave as point masses223 and thus define the appropriate two-body problem for the Wisdom-Holman224 mapping. The subsystems are, in general, not point masses and one225 takes this into account in the perturbation step of the Wisdom-Holman226 mapping.227 228 Instead of a simple chain of Jacobi coordinates, the code uses a229 tree of Jacobi coordinates. This makes it possible to use the code230 e.g., in situations when the massive particles have sattelites, either231 massive ones or not. It also reduces, when carefully specified, the232 integration error since more of the perturbations are absorbed into233 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 which235 works like a simple finite automaton but we did not use lex, yacc or236 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, and248 Sun, Jupiter, Saturn, Uranus and Neptune (1 - 5)249 can described as250 (((((1 [6 - 2006]) 2) 3) 4) 5)251 252 Another example for the tree specification when test particles between253 Jupiter and Saturn are added with indeces from 6 to 1000 and another set254 of test particles is added between Saturn and Uranus with indices from255 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 indeces260 1,2,3,4 and 5, respectively, for which the tree is261 ((1 (2 3)) (4 5))262 263 264 7) Format of the input file265 Id method order step_size total_integration_time printing_time266 number_of_particles267 268 tree_description269 270 GMs of particles, possibly zeros271 272 Initial conditions: x,y,z, vx, vy, vz273 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 used281 to initialize CSN. One can put a dummy tree description of the form282 (1 2)283 in the input file when HWH is not used. This will satisfy the parser284 which processes the tree description but will not play any role in the285 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+4293 5294 ((((1 2) 3) 4) 5)295 296 0.295912208285591095e-03297 0.282534590952422643e-06 0.845971518568065874e-07298 0.129202491678196939e-07 0.152435890078427628e-07299 300 0.0 0.0 0.0 0.0 0.0 0.0301 302 -0.538420940678015203e+01 -0.831247035699103631e+00 -0.225095848657298592e+00303 0.109236311953937086e-02 -0.652329334256263223e-02 -0.282301443780311710e-02304 305 0.788989169798583756e+01 0.459570785756998301e+01 0.155843220515231762e+01306 -0.321720261683698921e-02 0.433063255278440598e-02 0.192641782724487106e-02307 308 -0.182698980946940850e+02 -0.116271406901560681e+01 -0.250371938307287545e+00309 0.221540447608191936e-03 -0.376765389967669917e-02 -0.165324449144407023e-02310 311 -0.160595376822070470e+02 -0.239429495784517528e+02 -0.940042953888094424e+01312 0.264312302715172384e-02 -0.150349219713705076e-02 -0.681271071793930938e-03313 314 315 There are 2 output files:316 1) SomeId contains the actual integration data317 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 user319 has to redirect it in order to capture them.320 321 In the case of the Cowell-Stormer integrator, an additional file, SomeId.lerr322 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 printing324 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 printing328 time. The following is printed:329 time relative_change_in_total_energy relative_change_in_ang.mom.z.comp.330 after which331 x y z332 dx dy dz333 for each particle.334 335 336 9) The uset can change the code, of course. The section MAIN LOOP is where337 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 change344 it everywhere to DIM. It is still better to keep345 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 */353 5 354 6 /* dimension of physical space */ … … 381 33 /* cange this line accordingly */ 382 34 EXSEQBUL 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 AND393 PARTICLES WITH FINITE MASS394 You can add more types here, these values are stored395 in tpstore for each particle396 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 ARRAYS400 WHICH HOLD TEST PARTICLE INDICES.401 */402 35 403 36 #define NOTTESTPARTICLE 1
Note: See TracChangeset
for help on using the changeset viewer.