Nbody6
 All Files Functions Variables
adjust.f
Go to the documentation of this file.
1  SUBROUTINE adjust
2 *
3 *
4 * Parameter adjustment and energy check.
5 * --------------------------------------
6 *
7  include 'common6.h'
8  common/echain/ ech
9  SAVE dtoff
10  DATA dtoff /100.0d0/
11 *
12 *
13 * Predict X & XDOT for all particles (except unperturbed pairs).
14  CALL xvpred(ifirst,ntot)
15 *
16 * Obtain the total energy at current time (resolve all KS pairs).
17  CALL energy
18 *
19 * Initialize c.m. terms.
20  DO 10 k = 1,3
21  cmr(k) = 0.0d0
22  cmrdot(k) = 0.0d0
23  10 CONTINUE
24 *
25 * Obtain c.m. & angular momentum integrals and Z-moment of inertia.
26  az = 0.0d0
27  zm = 0.0d0
28  zmass = 0.0d0
29  isun = 1
30  DO 20 i = 1,n
31  IF (name(i).EQ.1) isun = i
32  zmass = zmass + body(i)
33  DO 15 k = 1,3
34  cmr(k) = cmr(k) + body(i)*x(k,i)
35  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
36  15 CONTINUE
37 * RI2 = (X(1,I) - RDENS(1))**2 + (X(2,I) - RDENS(2))**2 +
38 * & (X(3,I) - RDENS(3))**2
39  az = az + body(i)*(x(1,i)*xdot(2,i) - x(2,i)*xdot(1,i))
40  zm = zm + body(i)*(x(1,i)**2 + x(2,i)**2)
41  20 CONTINUE
42 *
43 * Form c.m. coordinates & velocities (vectors & scalars).
44  DO 25 k = 1,3
45  cmr(k) = cmr(k)/zmass
46  cmrdot(k) = cmrdot(k)/zmass
47  25 CONTINUE
48 *
49  cmr(4) = sqrt(cmr(1)**2 + cmr(2)**2 + cmr(3)**2)
50  cmrdot(4) = sqrt(cmrdot(1)**2 + cmrdot(2)**2 + cmrdot(3)**2)
51 *
52 * Subtract the kinetic energy of c.m. due to possible cloud effects.
53  IF (kz(13).GT.0) zkin = zkin - 0.5*zmass*cmrdot(4)**2
54 *
55 * Include non-zero virial energy for Plummer potential and/or 3D case.
56  vir = pot - vir
57 * Note angular momentum term is contained in virial energy (#14=1/2).
58  q = zkin/vir
59  e(3) = zkin - pot
60 * Modify single particle energy by tidal energy (except pure 3D).
61  IF (kz(14).NE.3) THEN
62  e(3) = e(3) + etide
63  END IF
64 *
65 * Modify angular momentum integral using Chandrasekhar eq. (5.530).
66  IF (kz(14).EQ.1.OR.kz(14).EQ.2) THEN
67  az = az + 0.5*tidal(4)*zm
68  END IF
69 *
70 * Define crossing time using single particle energy (cf. option 14).
71  tcr = zmass**2.5/(2.0*abs(e(3)))**1.5
72 * Note: see below for better definition in Plummer or 3D orbit case.
73  IF (q.GT.1.0.AND.kz(14).LT.3) THEN
74  tcr = tcr*sqrt(2.0*q)
75  END IF
76 * Form provisional total energy.
77  etot = zkin - pot + etide
78 *
79 * Include KS pairs, triple, quad, mergers, collisions & chain.
80  etot = etot + ebin + esub + emerge + ecoll + emdot + ecdot
81 * Add any contributions from chain regularization.
82  IF (nch.GT.0) THEN
83  etot = etot + ech
84  END IF
85 *
86 * Update energies and form the relative error (divide by ZKIN or E(3)).
87  IF (time.LE.0.0d0) THEN
88  de = 0.0d0
89  be(1) = etot
90  be(3) = etot
91  delta1 = 0.0d0
92  ELSE
93  be(2) = be(3)
94  be(3) = etot
95  de = be(3) - be(2)
96  delta1 = de
97  detot = detot + de
98  de = de/max(zkin,abs(e(3)))
99 * Save sum of relative energy error for main output and accumulate DE.
100  error = error + de
101  errtot = errtot + de
102  END IF
103 *
104 * Set provisional half-mass radius.
105  rscale = 0.5*zmass**2/pot
106 *
107 * Determine average neighbour number and smallest neighbour sphere.
108  nnb = 0
109  rs0 = rscale
110  DO 30 i = ifirst,ntot
111  nnb = nnb + list(1,i)
112  IF (list(1,i).GT.0) rs0 = min(rs0,rs(i))
113  30 CONTINUE
114  nnb = nnb/(n - npairs)
115 *
116 * Use current value if minimum neighbour sphere not implemented.
117  IF (rsmin.EQ.0.0d0) rsmin = rs0
118 *
119 * Find density centre & core radius (Casertano & Hut, Ap.J. 298, 80).
120  IF (n-npairs.GE.20.AND.kz(29).EQ.0.AND.kz(5).NE.3) THEN
121  CALL core
122  ELSE
123  nc = n
124  zmc = zmass
125  rhod = 1.0
126  rhom = 1.0
127  rc = rscale
128  rc2 = rc**2
129  rc2in = 1.0/rc2
130  END IF
131 *
132 * Take the Sun as reference for plotting planetesimal disk members.
133  IF (kz(5).EQ.3) THEN
134  DO 35 k = 1,3
135  rdens(k) = x(k,isun)
136  35 CONTINUE
137 *
138 * Determine the eccentricity dispersion and total energy of disk.
139  disp2 = 0.0
140  edisk = 0.0
141  DO 40 i = 1,n
142  IF (name(i).EQ.1.OR.name(i).EQ.nzero) go to 40
143  ri2 = (x(1,i) - x(1,isun))**2 + (x(2,i) - x(2,isun))**2
144  vi2 = (xdot(1,i) - xdot(1,isun))**2 +
145  & (xdot(2,i) - xdot(2,isun))**2
146  rrdot = (x(1,i) - x(1,isun))*(xdot(1,i) - xdot(1,isun)) +
147  & (x(2,i) - x(2,isun))*(xdot(2,i) - xdot(2,isun))
148  ri = sqrt(ri2)
149  semi = 2.0/ri - vi2/(body(isun) + body(i))
150  semi = 1.0/semi
151  ecc2 = (1.0 - ri/semi)**2 +
152  & rrdot**2/(semi*(body(i) + body(isun)))
153  disp2 = disp2 + ecc2
154  edisk = edisk - 0.5*body(i)/semi
155  40 CONTINUE
156  disp = sqrt(disp2/float(n-2))
157  WRITE (35,42) ttot, disp, edisk
158  42 FORMAT (' ',f8.1,1p,e10.2,e12.4)
159  END IF
160 *
161 * Check optional sorting of Lagrangian radii & half-mass radius.
162  IF (kz(7).GT.0) THEN
163  CALL lagr(rdens)
164  END IF
165 *
166 * Scale average & maximum core density by the mean value.
167  rhod = 4.0*twopi*rhod*rscale**3/(3.0*zmass)
168  rhom = 4.0*twopi*rhom*rscale**3/(3.0*zmass)
169 *
170 * Adopt density contrasts of unity for hot system.
171  IF (kz(29).GT.0.AND.zkin.GT.pot) THEN
172  rhod = 1.0
173  rhom = 1.0
174  END IF
175 *
176 * Check optional determination of regularization parameters.
177  IF (kz(16).GT.0) THEN
178  rmin0 = rmin
179 *
180 * Form close encounter distance from scale factor & density contrast.
181  rmin = 4.0*rscale/(float(n)*rhod**0.3333)
182 * Include alternative expression based on core radius (experimental).
183  IF (kz(16).GT.1.AND.nc.LT.0.01*n) THEN
184  rmin = 0.05*rc/float(nc)**0.3333
185  END IF
186 * Use harmonic mean to reduce fluctuations (avoid initial value).
187  IF (time.GT.0.0d0) rmin = sqrt(rmin0*rmin)
188 * Impose maximum value for sufficient perturbers.
189  rmin = min(rmin,rsmin*gmin**0.3333)
190 * Define scaled DTMIN by RMIN & <M> and include ETAI for consistency.
191  dtmin = 0.04*sqrt(etai/0.02d0)*sqrt(rmin**3/bodym)
192 * Specify binding energy per unit mass of hard binary (impose Q = 0.5).
193  eclose = 4.0*max(zkin,abs(zkin - pot))/zmass
194 * Adopt central velocity as upper limit (avoids large kick velocities).
195  IF (2.0*zkin/zmass.GT.vc**2.AND.vc.GT.0.0) eclose = 2.0*vc**2
196  IF (q.GT.0.5) THEN
197  eclose = 0.5*eclose/q
198  END IF
199  eclose = min(eclose,1.0d0)
200  ksmag = 0
201  END IF
202 *
203 * Check optional modification of DTMIN, ECLOSE & TCR for hot system.
204  IF (kz(29).GT.0.AND.q.GT.1.0) THEN
205  dtmin = 0.04*sqrt(etai/0.02d0)*sqrt(rmin**3/bodym)
206  sigma2 = 2.0*zkin/zmass
207  vp2 = 4.0*bodym/rmin
208  dtmin = dtmin*sqrt((vp2 + sigma2/q)/(vp2 + 2.0d0*sigma2))
209  eclose = sigma2
210  tcr = 2.0*rscale/sqrt(sigma2)
211  END IF
212 *
213 * Set useful scalars for the integrator.
214  smin = 2.0*dtmin
215  rmin2 = rmin**2
216  rmin22 = 4.0*rmin2
217  ebh = -0.5*bodym*eclose
218  IF (time.LE.0.0d0) THEN
219  stepj = 0.01*(60000.0/float(n))**0.3333
220  IF (dmod(smax,dtk(10)).NE.0.0d0) THEN
221  WRITE (6,43) smax, smax/dtk(10)
222  43 FORMAT (' FATAL ERROR! SMAX SMAX/DTK(10) ',1p,2e10.2)
223  stop
224  END IF
225  END IF
226 * Adopt 2*RSMIN for neighbour sphere volume factor in routine REGINT.
227  rsfac = max(25.0/tcr,3.0d0*vc/(2.0d0*rsmin))
228 *
229 * Update density contrast factor for neighbour sphere modification.
230  IF (time.LE.0.0d0.OR.kz(40).EQ.0) THEN
231  alpha = float(nnbmax)*sqrt(0.08d0*rscale**3/float(n-npairs))
232  END IF
233 * Include optional stabilization to increase neighbour number.
234  IF (kz(40).EQ.1.AND.float(nnb).LT.0.5*nnbmax) THEN
235  fac = 1.0 + (0.5*nnbmax - nnb)/(0.5*float(nnbmax))
236  alpha = fac*alpha
237  END IF
238 *
239 * Define tidal radius for isolated system (2*RTIDE used in ESCAPE).
240  IF (kz(14).EQ.0) rtide = 10.0*rscale
241 * Redefine the crossing time for 3D cluster orbit or Plummer model.
242  IF ((kz(14).EQ.3.OR.kz(14).EQ.4).AND.zkin.GT.0.0) THEN
243  tcr = 2.0*rscale/sqrt(2.0*zkin/zmass)
244  END IF
245 *
246 * Print energy diagnostics & KS parameters.
247  icr = ttot/tcr
248  WRITE (6,45) ttot, q, de, be(3), rmin, dtmin, icr, delta1, e(3),
249  & detot
250  45 FORMAT (/,' ADJUST: TIME =',f8.2,' Q =',f5.2,' DE =',1p,e10.2,
251  & ' E =',0p,f10.6,' RMIN =',1p,e8.1,' DTMIN =',e8.1,
252  & ' TC =',0p,i5,' DELTA =',1p,e9.1,' E(3) =',0p,f10.6,
253  & ' DETOT =',f10.6)
254  CALL flush(6)
255 *
256 * Perform automatic error control (RETURN on restart with KZ(2) > 1).
257  CALL check(de)
258  IF (abs(de).GT.5.0*qe) go to 70
259 *
260 * Check for escaper removal.
261  IF (kz(23).GT.0) THEN
262  CALL escape
263  END IF
264 *
265 * Check correction for c.m. displacements.
266  IF (kz(31).GT.0) THEN
267  CALL cmcorr
268  END IF
269 *
270 * Include diagnostics for massive binary (bound or unbound initially).
271  IF (kz(5).EQ.4) THEN
272  ip = 0
273  DO 50 ipair = 1,npairs
274  IF (name(2*ipair-1).LE.2.OR.name(2*ipair).LE.2) THEN
275  ip = ipair
276  END IF
277  50 CONTINUE
278  IF (ip.GT.0) THEN
279  i1 = 2*ip - 1
280  i2 = i1 + 1
281  semi = -0.5*body(n+ip)/h(ip)
282  ecc2 = (1.0 - r(ip)/semi)**2 +
283  & tdot2(ip)**2/(semi*body(n+ip))
284  eb = body(i1)*body(i2)*h(ip)/body(n+ip)
285  WRITE (35,52) ttot, semi, eb, e(3), sqrt(ecc2),
286  & name(i1), name(i2)
287  52 FORMAT (' ',f8.1,1p,3e12.4,0p,f7.3,2i5)
288 * 52 FORMAT (' T A E EB ECL NAME ',F8.1,1P,3E12.4,0P,F7.3,2I5)
289  CALL flush(35)
290  END IF
291  END IF
292 *
293 * Allow a gradual decrease of NNBMAX due to escaper removal.
294  IF (kz(40).EQ.3) THEN
295  nnbmax = nbzero*sqrt(float(n)/float(nzero))
296  znbmax = 0.9*nnbmax
297  znbmin = max(0.2*nnbmax,1.0)
298  END IF
299 *
300 * Include optional fine-tuning of neighbour number (#40 >= 2).
301  IF (kz(40).GE.2) THEN
302  IF (nnb.LT.nnbmax/5.0) THEN
303  alpha = 1.1*alpha
304  ELSE
305  alpha = 0.9*alpha
306  END IF
307  END IF
308 *
309 * See whether standard output is due.
310  iout = 0
311  IF (time.GE.tnext) THEN
312  CALL output
313  iout = 1
314  END IF
315 *
316 * Include optional diagnostics for the hardest binary below ECLOSE.
317  IF (kz(33).GE.2.AND.iout.GT.0) THEN
318  hp = 0.0
319  ip = 0
320  DO 60 ipair = 1,npairs
321 * Skip outer (ghost) binary of quadruple system.
322  IF (h(ipair).LT.hp.AND.body(n+ipair).GT.0.0d0) THEN
323  hp = h(ipair)
324  ip = ipair
325  END IF
326  60 CONTINUE
327  IF (ip.GT.0.AND.hp.LT.-eclose) THEN
328  i1 = 2*ip - 1
329  i2 = i1 + 1
330  semi = -0.5*body(n+ip)/h(ip)
331  pb = days*semi*sqrt(semi/body(n+ip))
332  ecc2 = (1.0 - r(ip)/semi)**2 +
333  & tdot2(ip)**2/(semi*body(n+ip))
334  eb = body(i1)*body(i2)*h(ip)/body(n+ip)
335  WRITE (39,62) ttot, name(i1), name(i2), kstar(n+ip),
336  & list(1,i1), sqrt(ecc2), semi, pb, eb, e(3)
337  62 FORMAT (' BINARY: T NAME K* NP E A P EB E3 ',
338  & f8.1,2i6,2i4,f7.3,1p,2e10.2,0p,2f9.4)
339  CALL flush(39)
340  END IF
341  END IF
342 *
343 * Include optional KS reg of binaries outside standard criterion.
344  IF (kz(8).GT.0.AND.dmod(time,dtk(10)).EQ.0.0d0) THEN
345 * Note DMOD condition needed for CALL KSREG and CALL STEPS.
346  dtcl = 30.0*dtmin
347  rcl = 10.0*rmin
348  CALL sweep(dtcl,rcl)
349  END IF
350 *
351 * Update time for next adjustment.
352  tadj = tadj + dtadj
353 * Re-initialize marginal stability counter to avoid including old case.
354  nmtry = 0
355 *
356 * Check optional truncation of time.
357  IF (kz(35).GT.0.AND.time.GE.dtoff) THEN
358  CALL offset(dtoff)
359  END IF
360 *
361 * Obtain elapsed CPU time and update total since last output/restart.
362  CALL cputim(tcomp)
363  cputot = cputot + tcomp - cpu0
364  cpu0 = tcomp
365 *
366 * Save COMMON after energy check (skip TRIPLE, QUAD, CHAIN).
367  tdump = time
368  IF (kz(2).GE.1.AND.nsub.EQ.0) CALL mydump(1,2)
369 *
370 * Check termination criteria (TIME > TCRIT & N <= NCRIT).
371  IF (ttot.GE.tcrit.OR.n.LE.ncrit.OR.ttot+dtadj.GT.tcrit) THEN
372 * Terminate after optional COMMON save.
373  WRITE (6,65) ttot, cputot/60.0, errtot, detot
374  65 FORMAT (//,9x,'END RUN',3x,'TIME =',f7.1,' CPUTOT =',f7.1,
375  & ' ERRTOT =',f10.6,' DETOT =',f10.6)
376  IF (kz(1).GT.0.AND.nsub.EQ.0) CALL mydump(1,1)
377  stop
378  END IF
379 *
380  70 RETURN
381 *
382  END