Nbody6
 All Files Functions Variables
evolve.f
Go to the documentation of this file.
1  SUBROUTINE evolve(IPAIR,LEVEL)
2 *
3 *
4 * Binary diagnostics.
5 * -------------------
6 *
7  include 'common6.h'
8  common/gamdot/ dgam
9  REAL*8 xcm(3),vcm(3),work(10)
10  INTEGER iflag(kmax),iout(lmax)
11 *
12 *
13  zkt = -be(1)/float(nzero)
14  m = 0
15  IF (level.EQ.0) go to 500
16 * NBSTAT = NBSTAT + 1
17 * Prepare diagnostic output for regularized binary.
18 *
19  DO 10 ip = 1,npairs
20  iflag(ip) = 0
21  10 CONTINUE
22  jpair = ipair
23 *
24  20 CALL resolv(jpair,1)
25  icm = n + jpair
26  j2 = 2*jpair
27  j1 = j2 - 1
28  semi = -0.5*body(icm)/h(jpair)
29  zn = sqrt(body(icm)/abs(semi)**3)
30  rdot = (x(1,j1) - x(1,j2))*(xdot(1,j1) - xdot(1,j2))
31  & + (x(2,j1) - x(2,j2))*(xdot(2,j1) - xdot(2,j2))
32  & + (x(3,j1) - x(3,j2))*(xdot(3,j1) - xdot(3,j2))
33  ecc = (1.0 - r(jpair)/semi)**2 + rdot**2/(semi*body(icm))
34  ecc = sqrt(ecc)
35  apo = abs(semi)*(1.0 + ecc)
36  ri = sqrt((x(1,icm) - rdens(1))**2 + (x(2,icm) - rdens(2))**2 +
37  & (x(3,icm) - rdens(3))**2)
38  zk = 0.5*body(icm)*(xdot(1,icm)**2+xdot(2,icm)**2+xdot(3,icm)**2)
39  zk = zk/zkt
40  eb = body(j1)*body(j2)*h(jpair)/body(icm)
41  eb = -eb/zkt
42  tk = twopi/zn
43  nk = tcr0/tk
44 *
45 * Print time interval since last binary output and set new time.
46  dt = 1.0e+06
47  k = kz(4)
48  DO 40 l = 1,k
49  IF (time - tlastb(l).LT.dt) dt = time - tlastb(l)
50  40 CONTINUE
51  tlastb(level) = time
52 *
53  IF (level.LT.0) dgam = 1.
54  gs = sign(gamma(jpair),dgam)
55  IF (name(icm).GT.0) THEN
56  WRITE (4,50) level, time, name(j1), name(j2), eb, semi, ecc,
57  & nk, r(jpair), gs, list(1,j1), ri, zk, dt, tk, m
58  ELSE
59  WRITE (4,50) level, time, nzero + name(j1), name(j2), eb,
60  & semi, ecc, nk, r(jpair), gs, list(1,j1), ri, zk,
61  & dt, tk, m, name(j1), -(name(icm) + nzero)
62  END IF
63 *
64  50 FORMAT (i3,f8.2,i6,i5,f7.1,f9.5,f6.2,i7,f9.5,f6.2,
65  & i5,f7.2,f6.1,2f10.6,5i5)
66  m = m + 1
67  iflag(jpair) = -1
68 *
69  nnb = list(1,j1)
70  DO 60 l = 1,nnb
71  iout(l) = 0
72  60 CONTINUE
73  rjmin2 = 1.e+06
74 *
75  nnb1 = nnb + 1
76  DO 150 l = 2,nnb1
77  j = list(l,j1)
78  rij2 = (x(1,icm) - x(1,j))**2 + (x(2,icm) - x(2,j))**2
79  & + (x(3,icm) - x(3,j))**2
80  vij2 = (xdot(1,icm) - xdot(1,j))**2 +
81  & (xdot(2,icm) - xdot(2,j))**2 +
82  & (xdot(3,icm) - xdot(3,j))**2
83  bb = body(icm) + body(j)
84  bbi = 1.0/bb
85  rij = sqrt(rij2)
86  eb = 0.5*vij2 - bb/rij
87  semij = -0.5*bb/eb
88  rdot = (x(1,icm) - x(1,j))*(xdot(1,icm) - xdot(1,j)) +
89  & (x(2,icm) - x(2,j))*(xdot(2,icm) - xdot(2,j)) +
90  & (x(3,icm) - x(3,j))*(xdot(3,icm) - xdot(3,j))
91  eccj = (1.0 - rij/semij)**2 + rdot**2/(semij*bb)
92  eccj = sqrt(eccj)
93  p = semij*(1.0 - eccj)
94  rdot = rdot/rij
95  eb = -body(j)*body(icm)*eb*bbi/zkt
96  geff = 2.0*body(j)*(abs(semi)/rij)**3/body(icm)
97  IF (geff.GT.99.9) geff = 99.9
98 * Exclude apocentres > A with A = 2/N for hard binary.
99  apo = min(apo,2.0d0/float(n))
100 *
101 * Print small impact parameter & significant effective perturbation.
102  IF (p.LT.3.*apo.AND.geff.GT.gprint(1)) THEN
103  WRITE (4,100) name(j), eb, semij, eccj, rij, geff, p,rdot
104  100 FORMAT (16x,i6,f7.1,f9.5,f6.2,f16.5,f6.2,f9.5,f7.2)
105  iout(l-1) = 1
106  IF (j.GT.n.AND.iflag(j-n).EQ.0) iflag(j-n) = 1
107  END IF
108 *
109 * Save parameters for the nearest perturber.
110  IF (rij2.LT.rjmin2) THEN
111  rjmin2 = rij2
112  lmin = l
113  work(1) = eb
114  work(2) = semij
115  work(3) = eccj
116  work(4) = rij
117  work(5) = geff
118  work(6) = p
119  work(7) = rdot
120  END IF
121  150 CONTINUE
122 *
123  IF (iout(lmin-1).EQ.0) THEN
124  jmin = list(lmin,j1)
125  WRITE (4,100) name(jmin), (work(k),k=1,7)
126  IF (jmin.GT.n.AND.iflag(jmin-n).EQ.0) iflag(jmin-n) = 1
127  END IF
128 *
129  DO 190 ip = 1,npairs
130  IF(iflag(ip).GT.0) THEN
131  jpair = ip
132 * Only include full output for active binaries.
133  IF (gamma(jpair).GT.gmax) go to 20
134  END IF
135  190 CONTINUE
136 *
137 * Set output times to facilitate DTCRIT estimate in routine KSINT.
138  kk = kz(4)
139  DO 200 k = 1,kk
140  tlastb(k) = time
141  200 CONTINUE
142  go to 900
143 *
144 * Search for binaries of single particles & c.m. bodies.
145  500 simax = 0.01*tcr0
146  IF (ipair.EQ.0) THEN
147  ilow = ifirst
148  ihigh = ntot
149  ELSE
150  ilow = ipair
151  ihigh = ipair
152  END IF
153  nam = ipair
154  IF (nam.GT.0) nam = name(nam)
155 *
156  DO 600 i = ilow,ihigh
157  IF (step(i).GT.simax) go to 600
158  j2 = 0
159  rij = 1.0e+06
160  nnb1 = list(1,i) + 1
161 *
162  DO 250 l = 2,nnb1
163  j = list(l,i)
164  IF (step(j).GT.simax) go to 250
165  rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
166  & (x(3,i) - x(3,j))**2
167  IF (rij2.GT.rij) go to 250
168  rij = rij2
169  j2 = j
170  250 CONTINUE
171 *
172  IF (j2.LT.i.AND.ipair.EQ.0) go to 600
173  rij = sqrt(rij)
174  vij2 = (xdot(1,i) - xdot(1,j2))**2 +
175  & (xdot(2,i) - xdot(2,j2))**2 +
176  & (xdot(3,i) - xdot(3,j2))**2
177  bodyij = body(i) + body(j2)
178  bbi = 1.0/bodyij
179  erel = 0.5*vij2 - bodyij/rij
180  IF (erel.GT.-0.1*eclose) go to 600
181  semi = -0.5*bodyij/erel
182  zn = sqrt(bodyij/semi**3)
183  rdot = (x(1,i) - x(1,j2))*(xdot(1,i) - xdot(1,j2)) +
184  & (x(2,i) - x(2,j2))*(xdot(2,i) - xdot(2,j2)) +
185  & (x(3,i) - x(3,j2))*(xdot(3,i) - xdot(3,j2))
186  ecc = (1.0 - rij/semi)**2 + rdot**2/(semi*bodyij)
187  ecc = sqrt(ecc)
188 *
189  DO 520 k = 1,3
190  xcm(k) = (body(i)*x(k,i) + body(j2)*x(k,j2))*bbi
191  vcm(k) = (body(i)*xdot(k,i) + body(j2)*xdot(k,j2))*bbi
192  520 CONTINUE
193 *
194  ri = sqrt((xcm(1) - rdens(1))**2 + (xcm(2) - rdens(2))**2 +
195  & (xcm(3) - rdens(3))**2)
196  zk = 0.5*bodyij*(vcm(1)**2 + vcm(2)**2 + vcm(3)**2)
197  zk = zk/zkt
198  eb = body(i)*body(j2)*erel*bbi
199  eb = -eb/zkt
200  tk = twopi/zn
201  nk = tcr0/tk
202  gb = 0.0
203 *
204 * Estimate the relative perturbation at apocentre.
205  dt = t0(i) - t0(j2)
206  DO 540 k = 1,3
207  gb = gb + (f(k,i) - (f(k,j2) + 3.0*fdot(k,j2)*dt))**2
208  540 CONTINUE
209  gb = 2.0*sqrt(gb)*rij**2*bbi
210  gb = gb*(semi*(1.0 + ecc)/rij)**3
211  IF (gb.GT.99.9) gb = 99.9
212 *
213 * Set time interval since last binary output.
214  dt = 1.0e+06
215  k = kz(4)
216  DO 550 l = 1,k
217  IF (time - tlastb(l).LT.dt) dt = time - tlastb(l)
218  550 CONTINUE
219 *
220 * Only output soft binaries for dominant motion at apocentre.
221  IF (eb.LT.1.0.AND.gb.GT.0.5) go to 600
222 *
223  IF (i.LE.n.AND.j2.LE.n) THEN
224  WRITE (4,50) level, time, name(i), name(j2), eb, semi,
225  & ecc, nk, rij, gb, list(1,i), ri, zk, dt, tk,
226  & nam
227  ELSE IF (i.LE.n.AND.j2.GT.n) THEN
228  jj2 = 2*(j2-n)
229  jj1 = jj2-1
230  WRITE (4,50) level, time, name(i), name(j2), eb, semi,
231  & ecc, nk, rij, gb, list(1,i), ri, zk, dt, tk,
232  & nam, name(jj1), name(jj2)
233  ELSE IF (i.GT.n.AND.j2.LE.n) THEN
234  ii2 = 2*(i-n)
235  ii1 = ii2-1
236  WRITE (4,50) level, time, name(i), name(j2), eb, semi,
237  & ecc, nk, rij, gb, list(1,i), ri, zk, dt, tk,
238  & nam, name(ii1), name(ii2)
239  ELSE
240  ii2 = 2*(i-n)
241  ii1 = ii2-1
242  jj2 = 2*(j2-n)
243  jj1 = jj2-1
244  WRITE (4,50) level, time, name(i), name(j2), eb, semi,
245  & ecc, nk, rij, gb, list(1,i), ri, zk, dt, tk,
246  & nam, name(ii1), name(ii2), name(jj1),
247  & name(jj2)
248  END IF
249  600 CONTINUE
250 *
251  IF (ipair.EQ.0) THEN
252  tlasts = time
253  ELSE
254  tlastt = time
255  END IF
256 *
257  900 RETURN
258 *
259  END