Nbody6
 All Files Functions Variables
quad.f
Go to the documentation of this file.
1  SUBROUTINE quad(ISUB)
2 *
3 *
4 * Four-body regularization.
5 * -------------------------
6 *
7 * Method of Mikkola & Aarseth, Celestial Mechanics 47, 375.
8 * .........................................................
9 *
10 *
11  IMPLICIT REAL*8 (a-h,o-z)
12  parameter(ncmax=10)
13  REAL*8 m,mij,mb,mb1,savex(12),savexd(12)
14  LOGICAL switch,gtype,gtype0
15  common/creg/ m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
16  & xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
17  & ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
18  common/tpr/ switch,gtype,gtype0
19  common/config/ r2(4,4),i1,i2,i3,i4
20  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
21  common/ccoll/ qk(12),pk(12),icall,icoll,ndiss4
22  common/bssave/ ep(4),dsc,facm,tfac,itfac,jc
23  common/ksave/ k1,k2
24  common/ebsave/ ebs
25  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
26  & names(ncmax,5),isys(5)
27  SAVE
28 *
29 *
30 * Main variables for four-body regularization
31 * *******************************************
32 *
33 * -------------------------------------------------------------------
34 * CM(1-7) Coordinates, velocities & total mass of four-body system.
35 * CM(8) Total energy of N-body system (copied in routine START4).
36 * CM(9) Binary energy change (copied to BBCOL in START4).
37 * DS Regularized time-step (rescaled after switching & new TPR).
38 * ENERGY Total energy.
39 * EPS Tolerance for DIFSY4 (1.0E-08 is recommended).
40 * EPSR2 Distance criterion for stabilized time transformation.
41 * ICALL Pericentre indicator (activated if R(IMIN) < EPSR2**0.5).
42 * ICOLL Collision indicator (activated in routine DERQP4).
43 * IND Sorting index (R(IND(6)) is largest distance).
44 * I1-I4 Indices of final configuration (I1 & I2 is closest binary).
45 * KZ15 Diagnostic option (copied in START4; full output if > 1).
46 * M Particle mass (CM(7) = M(1) + M(2) + M(3) + M(4)).
47 * NAME4 Particle identity (initialized to global name).
48 * NFN Number of function calls.
49 * NREG Number of regularization switches.
50 * NSTEP4 Number of DIFSY4 calls.
51 * P Regularized momenta.
52 * Q Regularized coordinates.
53 * QPERI Minimum two-body separation.
54 * R Mutual distances (R(1), R(2), R(3) are regularized).
55 * RCOLL Minimum two-body separation.
56 * RGRAV Gravitational radius ((sum M(I)*M(J))/ABS(ENERGY)).
57 * RIJ Minimum pairwise separations.
58 * RMAXS Maximum size of unperturbed configuration.
59 * RMAX4 Maximum of unperturbed & initial size (+ 20 per cent).
60 * R2 Square mutual distances (R2(I1,I2) is closest pair).
61 * TCR Local crossing time ((sum M(I))**2.5/ABS(2*ENERGY)**1.5).
62 * TIME4 Local physical time in scaled units.
63 * TMAX Maximum integration time (based on c.m. step).
64 * TSTEP Initial time-step for regularization (physical units).
65 * X Particle coordinates (X(1), X(2), X(3) is X, Y, Z).
66 * XD Velocity components.
67 * -------------------------------------------------------------------
68 *
69 *
70 * Save termination indicator and check for restart.
71  iterm = isub
72  IF (isub.GT.0) THEN
73 * Synchronize next time interval with c.m. step unless termination.
74  tmax = time4 + steps(isub)
75  IF (steps(isub).LE.0.0d0) ds = 0.01*ds
76  go to 30
77  END IF
78 *
79 * Obtain initial conditions from N-body COMMON.
80  CALL start4(isub)
81 *
82 * Specify tolerance & DSC for DIFSY4 (relative energy error < TOL).
83  eps = 1.0d-08
84  eps = 1.0d-12
85  dsc = 1.0
86 *
87 * Initialize diagnostic & COMMON variables.
88  time4 = 0.0d0
89  rcoll = 100.0
90  DO 10 j = 1,4
91  DO 5 k = 1,4
92  rij(j,k) = 100.0
93  5 CONTINUE
94  10 CONTINUE
95  icall = 0
96  icoll = 0
97  ndiss4 = 0
98  jc = 0
99  itfac = 0
100  nstep4 = 0
101  iq = 0
102  ecoll3 = 0.0d0
103  nreg = 0
104  nfn = 0
105  itry = 0
106  isave = 0
107  ratio = 0.0
108  switch = .false.
109 *
110 * Evaluate the initial total energy (superseded; done in START4).
111 * CALL NEWSYS(X,XD,M,4,ENERGY,GAM)
112 *
113 * Initialize regularized coordinates, momenta & particle indices.
114  CALL newreg
115 *
116 * Define gravitational radius & average mass factor for routine DERQP4.
117  rgrav = (mij(1) + mij(2) + mij(3) + mij(4) + mij(5) + mij(6))/
118  & abs(energy)
119  facm = rgrav*abs(energy)/6.0
120 * Set small distance criterion for stabilized time transformation.
121  epsr2 = (0.2*rgrav)**2
122 *
123 * Ensure that unperturbed boundary exceeds system size (+ 20 per cent).
124  rmax = max(sqrt(r2(i1,i3)),sqrt(r2(i2,i4)))
125  rmax4 = max(rmaxs(isub),1.2*rmax)
126 *
127 * Specify local crossing time.
128  tcr = (m(1) + m(2) + m(3) + m(4))**2.5/abs(2.0d0*energy)**1.5
129 *
130 * Define a nominal crossing time in near-parabolic cases.
131  tstar = rmax*sqrt(rmax/cm(7))
132 *
133 * Determine the smallest two-body time-scale from parabolic orbit.
134  rm = sqrt(r2(i1,i2))
135  vp2 = 2.0*(m(i1) + m(i2))/rm
136  tp = rm/sqrt(vp2)
137  tstar = min(tp,tstar)
138 *
139 * Set initial time-step in physical units (also used by RCHAIN).
140  tstep = min(tcr,tstar)*eps**0.1
141 *
142 * Convert physical step to regularized time.
143  ds = tstep/(r(1)*r(2)*r(3))
144 *
145 * Form initial binary energies (scaled by internal energy).
146  vrel2 = 0.0d0
147  vrel21 = 0.0d0
148  DO 25 k = 1,3
149  j1 = 3*(i1-1) + k
150  j2 = 3*(i2-1) + k
151  j3 = 3*(i3-1) + k
152  j4 = 3*(i4-1) + k
153  vrel2 = vrel2 + (xd(j1) - xd(j2))**2
154  vrel21 = vrel21 + (xd(j3) - xd(j4))**2
155  25 CONTINUE
156  eb0 = 2.0d0/sqrt(r2(i1,i2)) - vrel2/(m(i1) + m(i2))
157  eb10 = 2.0d0/sqrt(r2(i3,i4)) - vrel21/(m(i3) + m(i4))
158  eb0 = -0.5d0*m(i1)*m(i2)*eb0/energy
159  eb10 = -0.5d0*m(i3)*m(i4)*eb10/energy
160 *
161 * Perform regularized integration until termination is reached.
162  30 CALL rchain(iq)
163 *
164 * See whether temporary or actual termination (otherwise R > RMAX4).
165  IF (time4.GT.tmax) THEN
166  IF (steps(isub).GT.0.0d0) go to 100
167  END IF
168 *
169 * Transform to physical variables and order particle indices.
170  switch = .false.
171  CALL endreg
172  CALL status(x,j1,j2,j3,j4)
173 *
174 * Save the current configuration.
175  DO 32 j = 1,12
176  savex(j) = x(j)
177  savexd(j) = xd(j)
178  32 CONTINUE
179  tsave = time4
180  ds0 = ds
181 *
182 * Check for collision during last step.
183  IF (icoll.LE.0) go to 40
184 *
185 * Restore the minimum configuration.
186  DO 35 k = 1,12
187  q(k) = qk(k)
188  p(k) = pk(k)
189  35 CONTINUE
190 *
191 * Transform to physical variables.
192  switch = .false.
193  CALL endreg
194 *
195 * Order particle indices of the collision configuration.
196 * CALL STATUS(X,J1,J2,J3,J4)
197 *
198 * Distinguish between tidal energy loss and collision (ICOLL holds IM).
199  im = icoll
200  icoll = 0
201  IF (qperi.LT.4.0*max(SIZE(k1),SIZE(k2))) THEN
202  IF (qperi.LT.0.75*(SIZE(k1) + SIZE(k2))) THEN
203 *
204 * Obtain global coordinates & velocities (ITERM < 0 denotes collision).
205  iterm = -1
206  isub = -isub
207  CALL start4(isub)
208 *
209 * Combine internal energy and external c.m. kinetic energy.
210  h4 = energy + 0.5d0*cm(7)*(cm(4)**2 + cm(5)**2 + cm(6)**2)
211 *
212 * Determine dominant two-body energy from non-singular expressions.
213  CALL erel4(im,ebs,semi)
214  dminc = min(rcoll,dminc)
215 *
216 * Form composite body and begin KS regularization of closest pair.
217  CALL cmbody(h4,4)
218  go to 100
219  ELSE
220 *
221 * Modify regularized variables and check stability (ITERM < 0).
222  CALL qpmod4(im,iterm)
223  switch = .true.
224  CALL newreg
225  icall = 0
226  iq = iterm
227  IF (iterm.EQ.-1) go to 100
228  go to 30
229  END IF
230  END IF
231 *
232 * Form output diagnostics (only hard binary printed if KZ15 < 2).
233  40 vrel2 = 0.0d0
234  vrel21 = 0.0d0
235  rdot = 0.0d0
236 *
237  DO 50 k = 1,3
238  j1 = 3*(i1-1) + k
239  j2 = 3*(i2-1) + k
240  j3 = 3*(i3-1) + k
241  j4 = 3*(i4-1) + k
242  vrel2 = vrel2 + (xd(j1) - xd(j2))**2
243  vrel21 = vrel21 + (xd(j3) - xd(j4))**2
244  rdot = rdot + (x(j1) - x(j2))*(xd(j1) - xd(j2))
245  50 CONTINUE
246 *
247 * Evaluate orbital elements.
248  rb = sqrt(r2(i1,i2))
249  rb1 = sqrt(r2(i3,i4))
250  r13 = sqrt(r2(i1,i3))
251  r24 = sqrt(r2(i2,i4))
252  mb = m(i1) + m(i2)
253  semi = 2.0d0/rb - vrel2/mb
254  semi = 1.0/semi
255  ecc = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
256  eb = -m(i1)*m(i2)/(2.0d0*semi*energy)
257  et = -m(i1)*m(i2)/(2.0d0*semi*cm(8))
258  mb1 = m(i3) + m(i4)
259  semi1 = 2.0d0/rb1 - vrel21/mb1
260  semi1 = 1.0/semi1
261  eb1 = -m(i3)*m(i4)/(2.0d0*semi1*energy)
262 * Scale binding energy of relative motion by initial value.
263  e1 = (1.0 - eb - eb1)/eb0
264 *
265 * Determine nearest and most distant binary perturber.
266  IF (r13.LT.r24) THEN
267  im = i3
268  rm = r13
269  imax = i4
270  rmax = r24
271  ELSE
272  im = i4
273  rm = r24
274  imax = i3
275  rmax = r13
276  END IF
277 *
278 * Estimate perturbations on smallest binary.
279  gb = 2.0*mb1/mb*(rb/rm)**3
280  IF (eb1.LT.0.0) gb = gb*m(im)/mb1
281  g4 = 2.0*m(imax)/(mb + m(im))*(rm/rmax)**3
282  db = (eb - eb0)/eb0
283  IF (itry.LT.0) go to 70
284 *
285 * Find the optimum configuration for two binaries or compact triple.
286  itry = itry + 1
287  IF (rm/semi.LT.ratio) go to 60
288  IF (rb.LT.semi.OR.rb1.LT.semi1) go to 60
289 * Avoid new regularization of any binary near small pericentre.
290  ratio = rm/semi
291  isave = itry
292  60 IF (ratio.LT.5.0.AND.itry.LE.5.
293  & or.ratio.LT.3.0.AND.itry.LE.10) THEN
294  go to 30
295  END IF
296 *
297  IF (isave.EQ.itry) go to 70
298  itry = -1
299 * Restore the maximum triple configuration and obtain indices.
300  DO 65 j = 1,12
301  x(j) = savex(j)
302  xd(j) = savexd(j)
303  65 CONTINUE
304  time4 = tsave
305  CALL newreg
306  ds = 0.5*ds0
307  go to 40
308 *
309 * Postpone temination by one small step unless restart case.
310  70 IF (steps(isub).GT.0.0d0) THEN
311  steps(isub) = 0.0d0
312  go to 100
313  END IF
314 *
315 * Print final binary if relative energy increase exceeds 0.1.
316  IF (db.GT.0.1) THEN
317  WRITE (6,80) mb, semi, ecc, eb, gb, g4, rb1, eb1, e1, et
318  80 FORMAT (/,' QUAD BINARY',' MB =',f7.4,' A =',1p,e8.1,
319  & ' E =',0p,f5.2,' EB =',f7.4,' GB =',1p,e8.1,' G4 =',e8.1,
320  & ' RB1 =',e8.1,' EB1 =',0p,f5.2,' E1 =',f5.2,' ET =',f6.3)
321  END IF
322 *
323 * Check optional print diagnostics of four-body regularization.
324  IF (kz15.GT.1) THEN
325  tc = time4/tcr
326  ec = energy/cm(8)
327  WRITE (6,90) i1, i2, i3, i4, rb, r13, r24, rgrav, tc, nstep4,
328  & nreg, db, ec
329  90 FORMAT (/,' END QUAD ',4i3,' RB =',1p,e8.1,' R13 =',e8.1,
330  & ' R24 =',e8.1,' RG =',e8.1,' TC =',0p,f5.1,' #',
331  & i5,i4,' DB =',f5.2,' EC =',f6.3)
332  END IF
333 *
334 * Form binary energy change (set in BBCOLL; routine START4).
335  cm(9) = (eb - eb0 - eb10)*energy
336  IF (eb1.GT.0.0) cm(9) = cm(9) + eb1*energy
337 *
338 * Transform to global variables and begin new KS (I1 & I2), I3 & I4.
339  CALL start4(isub)
340 *
341 * Activate termination index for routine INTGRT.
342  iterm = -1
343 *
344 * Update current time unless termination and set subsystem index.
345  100 IF (iterm.GE.0) ts(isub) = t0s(isub) + time4
346  isub = iterm
347 *
348  RETURN
349 *
350  END