-
Notifications
You must be signed in to change notification settings - Fork 0
/
hermite.cpp
605 lines (548 loc) · 21.4 KB
/
hermite.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
// Time-stamp: <2002-01-18 21:51:36 piet>
//================================================================
// |
// /__----__ ........ |
// . \ ....: :. |
// : _\|/_ ..: |
// : /|\ : _\|/_ |
// ___ ___ _____ ___ /|\ |
// / | \ /\ \ / | | \ / | /\ | \ |
// | __ |___/ / \ \ / | | \/ | / \ |___/ |
// | | | \ /____\ \ / | | / | /____\ | \ \/ |
// \___| | \ / \ V | | / |____ / \ |___/ | |
// / |
// : _/| :.. |/ |
// :.. ____/ :.... .. |
/* o // :. _\|/_ / :........: |
* O `//\ /|\ |
* | /\ |
*=============================================================================
*
* nbody_sh1.C: an N-body integrator with a shared but variable time step
* (the same for all particles but changing in time), using
* the Hermite integration scheme.
*
* ref.: Hut, P., Makino, J. & McMillan, S., 1995,
* Astrophysical Journal Letters 443, L93-L96.
*
* note: in this first version, all functions are included in one file,
* without any use of a special library or header files.
*_____________________________________________________________________________
*
* usage: nbody_sh1 [-h (for help)] [-d step_size_control_parameter]
* [-e diagnostics_interval] [-o output_interval]
* [-t total_duration] [-i (start output at t = 0)]
* [-x (extra debugging diagnostics)]
*
* "step_size_control_parameter" is a coefficient determining the
* the size of the shared but variable time step for all particles
*
* "diagnostics_interval" is the time between output of diagnostics,
* in the form of kinetic, potential, and total energy; with the
* -x option, a dump of the internal particle data is made as well
*
* "output_interval" is the time between successive snapshot outputs
*
* "total_duration" is the integration time, until the program stops
*
* Input and output are written from the standard i/o streams. Since
* all options have sensible defaults, the simplest way to run the code
* is by only specifying the i/o files for the N-body snapshots:
*
* nbody_sh1 < data.in > data.out
*
* The diagnostics information will then appear on the screen.
* To capture the diagnostics information in a file, capture the
* standard error stream as follows:
*
* (nbody_sh1 < data.in > data.out) >& data.err
*
* Note: if any of the times specified in the -e, -o, or -t options are not an
* an integer multiple of "step", output will occur slightly later than
* predicted, after a full time step has been taken. And even if they
* are integer multiples, round-off error may induce one extra step.
*_____________________________________________________________________________
*
* External data format:
*
* The program expects input of a single snapshot of an N-body system,
* in the following format: the number of particles in the snapshot n;
* the time t; mass mi, position ri and velocity vi for each particle i,
* with position and velocity given through their three Cartesian
* coordinates, divided over separate lines as follows:
*
* n
* t
* m1 r1_x r1_y r1_z v1_x v1_y v1_z
* m2 r2_x r2_y r2_z v2_x v2_y v2_z
* ...
* mn rn_x rn_y rn_z vn_x vn_y vn_z
*
* Output of each snapshot is written according to the same format.
*
* Internal data format:
*
* The data for an N-body system is stored internally as a 1-dimensional
* array for the masses, and 2-dimensional arrays for the positions,
* velocities, accelerations and jerks of all particles.
*_____________________________________________________________________________
*
* version 1: Jan 2002 Piet Hut, Jun Makino
*_____________________________________________________________________________
*/
#include <iostream>
#include <cmath> // to include sqrt(), etc.
#include <cstdlib> // for atoi() and atof()
#include <unistd.h> // for getopt()
using namespace std;
typedef double real; // "real" as a general name for the
// standard floating-point data type
const int NDIM = 3; // number of spatial dimensions
void correct_step(real pos[][NDIM], real vel[][NDIM],
const real acc[][NDIM], const real jerk[][NDIM],
const real old_pos[][NDIM], const real old_vel[][NDIM],
const real old_acc[][NDIM], const real old_jerk[][NDIM],
int n, real dt);
void evolve(const real mass[], real pos[][NDIM], real vel[][NDIM],
int n, real & t, real dt_param, real dt_dia, real dt_out,
real dt_tot, bool init_out, bool x_flag);
void evolve_step(const real mass[], real pos[][NDIM], real vel[][NDIM],
real acc[][NDIM], real jerk[][NDIM], int n, real & t,
real dt, real & epot, real & coll_time);
void get_acc_jerk_pot_coll(const real mass[], const real pos[][NDIM],
const real vel[][NDIM], real acc[][NDIM],
real jerk[][NDIM], int n, real & epot,
real & coll_time);
void get_snapshot(real mass[], real pos[][NDIM], real vel[][NDIM], int n);
void predict_step(real pos[][NDIM], real vel[][NDIM],
const real acc[][NDIM], const real jerk[][NDIM],
int n, real dt);
void put_snapshot(const real mass[], const real pos[][NDIM],
const real vel[][NDIM], int n, real t);
bool read_options(int argc, char *argv[], real & dt_param, real & dt_dia,
real & dt_out, real & dt_tot, bool & i_flag, bool & x_flag);
void write_diagnostics(const real mass[], const real pos[][NDIM],
const real vel[][NDIM], const real acc[][NDIM],
const real jerk[][NDIM], int n, real t, real epot,
int nsteps, real & einit, bool init_flag,
bool x_flag);
/*-----------------------------------------------------------------------------
* main -- reads option values, reads a snapshot, and launches the
* integrator
*-----------------------------------------------------------------------------
*/
int main(int argc, char *argv[])
{
real dt_param = 0.03; // control parameter to determine time step size
real dt_dia = 1; // time interval between diagnostics output
real dt_out = 1; // time interval between output of snapshots
real dt_tot = 10; // duration of the integration
bool init_out = false; // if true: snapshot output with start at t = 0
// with an echo of the input snapshot
bool x_flag = false; // if true: extra debugging diagnostics output
if (! read_options(argc, argv, dt_param, dt_dia, dt_out, dt_tot, init_out,
x_flag))
return 1; // halt criterion detected by read_options()
int n; // N, number of particles in the N-body system
cin >> n;
real t; // time
cin >> t;
real * mass = new real[n]; // masses for all particles
real (* pos)[NDIM] = new real[n][NDIM]; // positions for all particles
real (* vel)[NDIM] = new real[n][NDIM]; // velocities for all particles
get_snapshot(mass, pos, vel, n);
evolve(mass, pos, vel, n, t, dt_param, dt_dia, dt_out, dt_tot, init_out,
x_flag);
delete[] mass;
delete[] pos;
delete[] vel;
}
/*-----------------------------------------------------------------------------
* read_options -- reads the command line options, and implements them.
*
* note: when the help option -h is invoked, the return value is set to false,
* to prevent further execution of the main program; similarly, if an
* unknown option is used, the return value is set to false.
*-----------------------------------------------------------------------------
*/
bool read_options(int argc, char *argv[], real & dt_param, real & dt_dia,
real & dt_out, real & dt_tot, bool & i_flag, bool & x_flag)
{
int c;
while ((c = getopt(argc, argv, "hd:e:o:t:ix")) != -1)
switch(c){
case 'h': cerr << "usage: " << argv[0]
<< " [-h (for help)]"
<< " [-d step_size_control_parameter]\n"
<< " [-e diagnostics_interval]"
<< " [-o output_interval]\n"
<< " [-t total_duration]"
<< " [-i (start output at t = 0)]\n"
<< " [-x (extra debugging diagnostics)]"
<< endl;
return false; // execution should stop after help
case 'd': dt_param = atof(optarg);
break;
case 'e': dt_dia = atof(optarg);
break;
case 'i': i_flag = true;
break;
case 'o': dt_out = atof(optarg);
break;
case 't': dt_tot = atof(optarg);
break;
case 'x': x_flag = true;
break;
case '?': cerr << "usage: " << argv[0]
<< " [-h (for help)]"
<< " [-d step_size_control_parameter]\n"
<< " [-e diagnostics_interval]"
<< " [-o output_interval]\n"
<< " [-t total_duration]"
<< " [-i (start output at t = 0)]\n"
<< " [-x (extra debugging diagnostics)]"
<< endl;
return false; // execution should stop after error
}
return true; // ready to continue program execution
}
/*-----------------------------------------------------------------------------
* get_snapshot -- reads a single snapshot from the input stream cin.
*
* note: in this implementation, only the particle data are read in, and it
* is left to the main program to first read particle number and time
*-----------------------------------------------------------------------------
*/
void get_snapshot(real mass[], real pos[][NDIM], real vel[][NDIM], int n)
{
for (int i = 0; i < n ; i++){
cin >> mass[i]; // mass of particle i
for (int k = 0; k < NDIM; k++)
cin >> pos[i][k]; // position of particle i
for (int k = 0; k < NDIM; k++)
cin >> vel[i][k]; // velocity of particle i
}
}
/*-----------------------------------------------------------------------------
* put_snapshot -- writes a single snapshot on the output stream cout.
*
* note: unlike get_snapshot(), put_snapshot handles particle number and time
*-----------------------------------------------------------------------------
*/
void put_snapshot(const real mass[], const real pos[][NDIM],
const real vel[][NDIM], int n, real t)
{
cout.precision(16); // full double precision
cout << n << endl; // N, total particle number
cout << t << endl; // current time
for (int i = 0; i < n ; i++){
cout << mass[i]; // mass of particle i
for (int k = 0; k < NDIM; k++)
cout << ' ' << pos[i][k]; // position of particle i
for (int k = 0; k < NDIM; k++)
cout << ' ' << vel[i][k]; // velocity of particle i
cout << endl;
}
}
/*-----------------------------------------------------------------------------
* write_diagnostics -- writes diagnostics on the error stream cerr:
* current time; number of integration steps so far;
* kinetic, potential, and total energy; absolute and
* relative energy errors since the start of the run.
* If x_flag (x for eXtra data) is true, all internal
* data are dumped for each particle (mass, position,
* velocity, acceleration, and jerk).
*
* note: the kinetic energy is calculated here, while the potential energy is
* calculated in the function get_acc_jerk_pot_coll().
*-----------------------------------------------------------------------------
*/
void write_diagnostics(const real mass[], const real pos[][NDIM],
const real vel[][NDIM], const real acc[][NDIM],
const real jerk[][NDIM], int n, real t, real epot,
int nsteps, real & einit, bool init_flag,
bool x_flag)
{
real ekin = 0; // kinetic energy of the n-body system
for (int i = 0; i < n ; i++)
for (int k = 0; k < NDIM ; k++)
ekin += 0.5 * mass[i] * vel[i][k] * vel[i][k];
real etot = ekin + epot; // total energy of the n-body system
if (init_flag) // at first pass, pass the initial
einit = etot; // energy back to the calling function
cerr << "at time t = " << t << " , after " << nsteps
<< " steps :\n E_kin = " << ekin
<< " , E_pot = " << epot
<< " , E_tot = " << etot << endl;
cerr << " "
<< "absolute energy error: E_tot - E_init = "
<< etot - einit << endl;
cerr << " "
<< "relative energy error: (E_tot - E_init) / E_init = "
<< (etot - einit) / einit << endl;
if (x_flag){
cerr << " for debugging purposes, here is the internal data "
<< "representation:\n";
for (int i = 0; i < n ; i++){
cerr << " internal data for particle " << i+1 << " : " << endl;
cerr << " ";
cerr << mass[i];
for (int k = 0; k < NDIM; k++)
cerr << ' ' << pos[i][k];
for (int k = 0; k < NDIM; k++)
cerr << ' ' << vel[i][k];
for (int k = 0; k < NDIM; k++)
cerr << ' ' << acc[i][k];
for (int k = 0; k < NDIM; k++)
cerr << ' ' << jerk[i][k];
cerr << endl;
}
}
}
/*-----------------------------------------------------------------------------
* evolve -- integrates an N-body system, for a total duration dt_tot.
* Snapshots are sent to the standard output stream once every
* time interval dt_out. Diagnostics are sent to the standard
* error stream once every time interval dt_dia.
*
* note: the integration time step, shared by all particles at any given time,
* is variable. Before each integration step we use coll_time (short
* for collision time, an estimate of the time scale for any significant
* change in configuration to happen), multiplying it by dt_param (the
* accuracy parameter governing the size of dt in units of coll_time),
* to obtain the new time step size.
*
* Before moving any particles, we start with an initial diagnostics output
* and snapshot output if desired. In order to write the diagnostics, we
* first have to calculate the potential energy, with get_acc_jerk_pot_coll().
* That function also calculates accelerations, jerks, and an estimate for the
* collision time scale, all of which are needed before we can enter the main
* integration loop below.
* In the main loop, we take as many integration time steps as needed to
* reach the next output time, do the output required, and continue taking
* integration steps and invoking output this way until the final time is
* reached, which triggers a `break' to jump out of the infinite loop set up
* with `while(true)'.
*-----------------------------------------------------------------------------
*/
void evolve(const real mass[], real pos[][NDIM], real vel[][NDIM],
int n, real & t, real dt_param, real dt_dia, real dt_out,
real dt_tot, bool init_out, bool x_flag)
{
cerr << "Starting a Hermite integration for a " << n
<< "-body system,\n from time t = " << t
<< " with time step control parameter dt_param = " << dt_param
<< " until time " << t + dt_tot
<< " ,\n with diagnostics output interval dt_dia = "
<< dt_dia << ",\n and snapshot output interval dt_out = "
<< dt_out << "." << endl;
real (* acc)[NDIM] = new real[n][NDIM]; // accelerations and jerks
real (* jerk)[NDIM] = new real[n][NDIM]; // for all particles
real epot; // potential energy of the n-body system
real coll_time; // collision (close encounter) time scale
get_acc_jerk_pot_coll(mass, pos, vel, acc, jerk, n, epot, coll_time);
int nsteps = 0; // number of integration time steps completed
real einit; // initial total energy of the system
write_diagnostics(mass, pos, vel, acc, jerk, n, t, epot, nsteps, einit,
true, x_flag);
if (init_out) // flag for initial output
put_snapshot(mass, pos, vel, n, t);
real t_dia = t + dt_dia; // next time for diagnostics output
real t_out = t + dt_out; // next time for snapshot output
real t_end = t + dt_tot; // final time, to finish the integration
while (true){
while (t < t_dia && t < t_out && t < t_end){
real dt = dt_param * coll_time;
evolve_step(mass, pos, vel, acc, jerk, n, t, dt, epot, coll_time);
nsteps++;
}
if (t >= t_dia){
write_diagnostics(mass, pos, vel, acc, jerk, n, t, epot, nsteps,
einit, false, x_flag);
t_dia += dt_dia;
}
if (t >= t_out){
put_snapshot(mass, pos, vel, n, t);
t_out += dt_out;
}
if (t >= t_end)
break;
}
delete[] acc;
delete[] jerk;
}
/*-----------------------------------------------------------------------------
* evolve_step -- takes one integration step for an N-body system, using the
* Hermite algorithm.
*-----------------------------------------------------------------------------
*/
void evolve_step(const real mass[], real pos[][NDIM], real vel[][NDIM],
real acc[][NDIM], real jerk[][NDIM], int n, real & t,
real dt, real & epot, real & coll_time)
{
real (* old_pos)[NDIM] = new real[n][NDIM];
real (* old_vel)[NDIM] = new real[n][NDIM];
real (* old_acc)[NDIM] = new real[n][NDIM];
real (* old_jerk)[NDIM] = new real[n][NDIM];
for (int i = 0; i < n ; i++)
for (int k = 0; k < NDIM ; k++){
old_pos[i][k] = pos[i][k];
old_vel[i][k] = vel[i][k];
old_acc[i][k] = acc[i][k];
old_jerk[i][k] = jerk[i][k];
}
predict_step(pos, vel, acc, jerk, n, dt);
get_acc_jerk_pot_coll(mass, pos, vel, acc, jerk, n, epot, coll_time);
correct_step(pos, vel, acc, jerk, old_pos, old_vel, old_acc, old_jerk,
n, dt);
t += dt;
delete[] old_pos;
delete[] old_vel;
delete[] old_acc;
delete[] old_jerk;
}
/*-----------------------------------------------------------------------------
* predict_step -- takes the first approximation of one Hermite integration
* step, advancing the positions and velocities through a
* Taylor series development up to the order of the jerks.
*-----------------------------------------------------------------------------
*/
void predict_step(real pos[][NDIM], real vel[][NDIM],
const real acc[][NDIM], const real jerk[][NDIM],
int n, real dt)
{
for (int i = 0; i < n ; i++)
for (int k = 0; k < NDIM ; k++){
pos[i][k] += vel[i][k]*dt + acc[i][k]*dt*dt/2
+ jerk[i][k]*dt*dt*dt/6;
vel[i][k] += acc[i][k]*dt + jerk[i][k]*dt*dt/2;
}
}
/*-----------------------------------------------------------------------------
* correct_step -- takes one iteration to improve the new values of position
* and velocities, effectively by using a higher-order
* Taylor series constructed from the terms up to jerk at
* the beginning and the end of the time step.
*-----------------------------------------------------------------------------
*/
void correct_step(real pos[][NDIM], real vel[][NDIM],
const real acc[][NDIM], const real jerk[][NDIM],
const real old_pos[][NDIM], const real old_vel[][NDIM],
const real old_acc[][NDIM], const real old_jerk[][NDIM],
int n, real dt)
{
for (int i = 0; i < n ; i++)
for (int k = 0; k < NDIM ; k++){
vel[i][k] = old_vel[i][k] + (old_acc[i][k] + acc[i][k])*dt/2
+ (old_jerk[i][k] - jerk[i][k])*dt*dt/12;
pos[i][k] = old_pos[i][k] + (old_vel[i][k] + vel[i][k])*dt/2
+ (old_acc[i][k] - acc[i][k])*dt*dt/12;
}
}
/*-----------------------------------------------------------------------------
* get_acc_jerk_pot_coll -- calculates accelerations and jerks, and as side
* effects also calculates potential energy and
* the time scale coll_time for significant changes
* in local configurations to occur.
* __ __
* | --> --> |
* M M | r . v |
* --> j --> --> j | --> ji ji --> |
* a == -------- r ; j == -------- | v - 3 --------- r |
* ji |--> |3 ji ji |--> |3 | ji |--> |2 ji |
* | r | | r | | | r | |
* | ji | | ji | |__ | ji | __|
*
* note: it would be cleaner to calculate potential energy and collision time
* in a separate function. However, the current function is by far the
* most time consuming part of the whole program, with a double loop
* over all particles that is executed every time step. Splitting off
* some of the work to another function would significantly increase
* the total computer time (by an amount close to a factor two).
*
* We determine the values of all four quantities of interest by walking
* through the system in a double {i,j} loop. The first three, acceleration,
* jerk, and potential energy, are calculated by adding successive terms;
* the last, the estimate for the collision time, is found by determining the
* minimum value over all particle pairs and over the two choices of collision
* time, position/velocity and sqrt(position/acceleration), where position and
* velocity indicate their relative values between the two particles, while
* acceleration indicates their pairwise acceleration. At the start, the
* first three quantities are set to zero, to prepare for accumulation, while
* the last one is set to a very large number, to prepare for minimization.
* The integration loops only over half of the pairs, with j > i, since
* the contributions to the acceleration and jerk of particle j on particle i
* is the same as those of particle i on particle j, apart from a minus sign
* and a different mass factor.
*-----------------------------------------------------------------------------
*/
void get_acc_jerk_pot_coll(const real mass[], const real pos[][NDIM],
const real vel[][NDIM], real acc[][NDIM],
real jerk[][NDIM], int n, real & epot,
real & coll_time)
{
for (int i = 0; i < n ; i++)
for (int k = 0; k < NDIM ; k++)
acc[i][k] = jerk[i][k] = 0;
epot = 0;
const real VERY_LARGE_NUMBER = 1e300;
real coll_time_q = VERY_LARGE_NUMBER; // collision time to 4th power
real coll_est_q; // collision time scale estimate
// to 4th power (quartic)
for (int i = 0; i < n ; i++){
for (int j = i+1; j < n ; j++){ // rji[] is the vector from
real rji[NDIM]; // particle i to particle j
real vji[NDIM]; // vji[] = d rji[] / d t
for (int k = 0; k < NDIM ; k++){
rji[k] = pos[j][k] - pos[i][k];
vji[k] = vel[j][k] - vel[i][k];
}
real r2 = 0; // | rji |^2
real v2 = 0; // | vji |^2
real rv_r2 = 0; // ( rij . vij ) / | rji |^2
for (int k = 0; k < NDIM ; k++){
r2 += rji[k] * rji[k];
v2 += vji[k] * vji[k];
rv_r2 += rji[k] * vji[k];
}
rv_r2 /= r2;
real r = sqrt(r2); // | rji |
real r3 = r * r2; // | rji |^3
// add the {i,j} contribution to the total potential energy for the system:
epot -= mass[i] * mass[j] / r;
// add the {j (i)} contribution to the {i (j)} values of acceleration and jerk:
real da[NDIM]; // main terms in pairwise
real dj[NDIM]; // acceleration and jerk
for (int k = 0; k < NDIM ; k++){
da[k] = rji[k] / r3; // see equations
dj[k] = (vji[k] - 3 * rv_r2 * rji[k]) / r3; // in the header
}
for (int k = 0; k < NDIM ; k++){
acc[i][k] += mass[j] * da[k]; // using symmetry
acc[j][k] -= mass[i] * da[k]; // find pairwise
jerk[i][k] += mass[j] * dj[k]; // acceleration
jerk[j][k] -= mass[i] * dj[k]; // and jerk
}
// first collision time estimate, based on unaccelerated linear motion:
coll_est_q = (r2*r2) / (v2*v2);
if (coll_time_q > coll_est_q)
coll_time_q = coll_est_q;
// second collision time estimate, based on free fall:
real da2 = 0; // da2 becomes the
for (int k = 0; k < NDIM ; k++) // square of the
da2 += da[k] * da[k]; // pair-wise accel-
double mij = mass[i] + mass[j]; // eration between
da2 *= mij * mij; // particles i and j
coll_est_q = r2/da2;
if (coll_time_q > coll_est_q)
coll_time_q = coll_est_q;
}
} // from q for quartic back
coll_time = sqrt(sqrt(coll_time_q)); // to linear collision time
}
/*-----------------------------------------------------------------------------
* \\ o
* end of file: nbody_sh1.C /\\' O
* /\ |
*=============================================================================
*/