Nbody6
recoil.f
Go to the documentation of this file.
1  SUBROUTINE recoil(IEND)
2 *
3 *
4 * Binary analysis and final output.
5 * ---------------------------------
6 *
7  include 'commonc.h'
8  include 'common2.h'
9  REAL*8 r2(nmx,nmx),rc(3),vc(3),xx(3,3),vv(3,3)
10  INTEGER ij(nmx)
11  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
12  common/cpert/ rgrav,gpert,ipert,npert
13  common/calls/ tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
14  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
15  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
16  common/clump/ bodys(nmx,5),t0s(5),ts(5),steps(5),rmaxs(5),
17  & names(nmx,5),isys(5)
18  SAVE
19 *
20 *
21 * Sort particle separations (I1 & I2 form closest pair).
22  CALL r2sort(ij,r2)
23  i1 = ij(1)
24  i2 = ij(2)
25  i3 = ij(3)
26  i4 = ij(4)
27  IF (n.EQ.2) THEN
28  i3 = i2
29  i4 = i1
30  ELSE IF (n.EQ.3) THEN
31  i4 = i1
32  ELSE IF (n.GE.4) THEN
33 * Determine indices of second closest pair (avoid pair I1-I2).
34  rx1 = 1.0
35  rx0 = r2(i1,i2)
36  DO 2 j1 = 1,n
37  IF (j1.EQ.i1.OR.j1.EQ.i2) go to 2
38  DO 1 j2 = j1+1,n
39  IF (j2.EQ.i1.OR.j2.EQ.i2) go to 1
40  IF (r2(j1,j2).LT.rx1.AND.r2(j1,j2).GT.rx0) THEN
41  rx1 = r2(j1,j2)
42  i3 = j1
43  i4 = j2
44  END IF
45  1 CONTINUE
46  2 CONTINUE
47  END IF
48  k1 = i3
49  k2 = i4
50 *
51 * Ensure that original KS pair is considered (H > 0 is possible).
52  IF (n.LE.3.AND.iend.EQ.0) THEN
53  DO 5 k = 1,n
54  IF (iname(k).EQ.1) k1 = k
55  IF (iname(k).EQ.2) k2 = k
56  5 CONTINUE
57  END IF
58 *
59 * Form output diagnostics.
60  vrel2 = 0.0d0
61  vrel21 = 0.0d0
62  vrel23 = 0.0d0
63  rdot = 0.0d0
64  rc2 = 0.0d0
65  vc2 = 0.0d0
66  rcdot = 0.0d0
67  mb = m(i1) + m(i2)
68 *
69  DO 10 k = 1,3
70  j1 = 3*(i1-1) + k
71  j2 = 3*(i2-1) + k
72  j3 = 3*(i3-1) + k
73  j4 = 3*(i4-1) + k
74  l1 = 3*(k1-1) + k
75  l2 = 3*(k2-1) + k
76  vrel2 = vrel2 + (v(j1) - v(j2))**2
77  vrel21 = vrel21 + (v(j3) - v(j4))**2
78  vrel23 = vrel23 + (v(l1) - v(l2))**2
79  rdot = rdot + (x(j1) - x(j2))*(v(j1) - v(j2))
80  rc(k) = (m(i1)*x(j1) + m(i2)*x(j2))/mb
81  vc(k) = (m(i1)*v(j1) + m(i2)*v(j2))/mb
82  rc2 = rc2 + (rc(k) - x(j3))**2
83  vc2 = vc2 + (vc(k) - v(j3))**2
84  rcdot = rcdot + (rc(k) - x(j3))*(vc(k) - v(j3))
85  xx(k,1) = x(j1)
86  xx(k,2) = x(j2)
87  xx(k,3) = x(j3)
88  vv(k,1) = v(j1)
89  vv(k,2) = v(j2)
90  vv(k,3) = v(j3)
91  10 CONTINUE
92 *
93 * Evaluate orbital elements.
94  rb = sqrt(r2(i1,i2))
95  rb1 = sqrt(r2(i3,i4))
96  r13 = sqrt(r2(i1,i3))
97  r24 = sqrt(r2(i2,i4))
98  semi = 2.0d0/rb - vrel2/mb
99  semi = 1.0/semi
100  ecc = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
101  eb = -m(i1)*m(i2)/(2.0d0*semi)
102  et = -m(i1)*m(i2)/(2.0d0*semi*cm(8))
103 *
104 * Determine parameters for possible hierarchy of body #I3.
105  rcp = sqrt(rc2)
106  a1 = 2.0/rcp - vc2/(mb + m(i3))
107  a1 = 1.0/a1
108  ecc1 = sqrt((1.0 - rcp/a1)**2 + rcdot**2/(a1*(mb + m(i3))))
109  pmin = a1*(1.0 - ecc1)
110 *
111 * Obtain binding energy of K1 & K2 (in case I1 & I2 not bound).
112  semi2 = 2.0d0/sqrt(r2(k1,k2)) - vrel23/(m(k1) + m(k2))
113  semi2 = 1.0/semi2
114  eb2 = -m(k1)*m(k2)/(2.0d0*semi2)
115 *
116 * Include possible second binary (suppress positive energy).
117  IF (n.GT.3) THEN
118  mb1 = m(i3) + m(i4)
119  semi1 = 2.0d0/rb1 - vrel21/mb1
120  semi1 = 1.0/semi1
121  eb1 = -m(i3)*m(i4)/(2.0d0*semi1)
122  eb1 = min(eb1,0.0d0)
123  ELSE
124  r24 = 1.0e+10
125  IF (timec.LE.0.0d0) eb1 = 0.0
126  END IF
127 *
128 * Ensure that perturbed boundary exceeds system size.
129  rmax = max(sqrt(r2(i1,i3)),sqrt(r2(i2,i4)))
130  rmaxc = max(rmaxc,1.2*rmax)
131 *
132 * Save initial energies (note: K1 & K2 may be most energetic).
133  IF (iend.EQ.0) THEN
134  e0 = energy
135  de = 0.0
136 * IF (TIMEC.GT.0.0D0) GO TO 50
137  ebch0 = cm(9)
138  eb0 = min(eb,eb2)
139  eb10 = eb1
140  a0 = semi
141  ecc0 = ecc
142  IF (eb2.LT.eb) a0 = semi2
143 * Note closest pair (#I1/I2 may be intruder or wide binary at peri).
144  IF (n.EQ.3) THEN
145  name1 = namec(1)
146  name2 = namec(2)
147  ELSE
148  name1 = namec(i1)
149  name2 = namec(i2)
150  END IF
151  go to 50
152  ELSE IF (iend.EQ.1) THEN
153  de = de + abs(energy - e0)
154  e0 = energy
155 * Initialize second binary on subsequent absorption.
156  IF (eb10.EQ.0.0d0) eb10 = eb1
157  go to 50
158  END IF
159 *
160 * Determine nearest and most distant binary perturber (N = 4).
161  IF (r13.LT.r24) THEN
162  im = i3
163  rm = r13
164  imax = i4
165  rmax = r24
166  ELSE
167  im = i4
168  rm = r24
169  imax = i3
170  rmax = r13
171  END IF
172 *
173 * Estimate perturbations on smallest binary.
174  IF (n.GT.3) THEN
175  IF (n.GT.4) rm = max(r13,r24)
176  gb = 2.0*mb1/mb*(rb/rm)**3
177  IF (eb1.LT.0.0) gb = gb*m(im)/mb1
178  g4 = 2.0*m(imax)/(mb + m(im))*(rm/rmax)**3
179  ELSE
180  gb = m(i3)/mb*(rb/r13)**3
181  g4 = 0.0
182  eb1 = 0.0
183  END IF
184 *
185 * Form net binary energy change (added to CHCOLL in routine CHTERM).
186  IF (iend.EQ.2) THEN
187  cm(9) = eb + eb1 - cm(9) + ecoll1
188  END IF
189 *
190 * Set relative energy production and total energy change.
191  db = cm(9)/ebch0
192  de = de + abs(energy - e0)
193 *
194 * Provide diagnostics for exchange (membership may have changed).
195  IF (name1 + name2.NE.namec(i1) + namec(i2).AND.eb.LT.0.0) THEN
196  isub = isys(5)
197  tch = t0s(isub) + timec
198  WRITE (6,15) tch, name1, name2, namec(i1), namec(i2), ecc0,
199  & ecc, a0, semi, eb0, eb
200  15 FORMAT (' EXCHANGE T NAM E0 E A0 A EB0 EB ',
201  & f9.2,4i6,2f7.3,1p,2e9.1,2e10.1)
202  END IF
203 *
204 * Print final binary if relative energy increase exceeds 0.1.
205  IF (db.GT.0.1.AND.semi.GT.0.0) THEN
206 * Scale binding energy of relative motion by initial value.
207  e1 = (energy - eb - eb1)/eb0
208  eb = eb/energy
209  eb1 = eb1/energy
210  WRITE (6,20) namec(i1), namec(i2), semi, ecc, eb, gb, g4,
211  & eb1, e1, et, db
212  20 FORMAT (' CHAIN RECOIL',' NAM =',2i6,' A =',1p,e8.1,
213  & ' E =',0p,f5.2,' EB =',f5.2,' GB =',1p,e8.1,
214  & ' G4 =',e8.1,' EB1 =',0p,f5.2,' E1 =',f5.2,
215  & ' ET =',f6.3,' DB =',f5.1)
216  END IF
217 *
218 * Include output of strong interaction from chaos case.
219  IF (iend.EQ.2.AND.pmin.LT.2.0*semi.AND.db.GT.0.1) THEN
220  CALL inclin(xx,vv,rc,vc,alpha)
221  WRITE (6,25) namec(i3), ecc, ecc1, pmin/semi, rcdot/rcp,
222  & semi, a1, rcp, gb, 180.0*alpha/3.14
223  25 FORMAT (' RECOIL: NAM E E1 PM/A RD A A1 RP GB IN ',
224  & i6,f8.4,f6.2,f5.1,f6.1,1p,4e10.2,0p,f7.1)
225  END IF
226 *
227 * Check optional print diagnostics of chain regularization.
228  IF (kz30.GT.1.AND.iend.EQ.2) THEN
229  tcr = mass**2.5/abs(2.0d0*energy)**1.5
230  tc = timec/tcr
231  ec = energy/cm(8)
232  WRITE (6,30) i1, i2, i3, i4, rb, r13, r24, de, tc, nstep1,
233  & nreg, npert, db, ec
234  30 FORMAT (/,' END CHAIN ',4i3,' RB =',1pe8.1,' R13 =',e8.1,
235  & ' R24 =',e8.1,' DE =',e8.1,' TC =',0p,f5.1,' #',
236  & i5,i4,i3,' DB =',f5.2,' EC =',f6.3)
237  END IF
238 *
239  50 RETURN
240 *
241  END