Nbody6
 All Files Functions Variables
hiarch.f
Go to the documentation of this file.
1  SUBROUTINE hiarch(IPAIR)
2 *
3 *
4 * Hierarchical system diagnostics.
5 * --------------------------------
6 *
7  include 'common6.h'
8  REAL*8 a1(3),a2(3),xrel(3),vrel(3),ei(3),hi(3),ho(3),
9  & tk0(100),tlast(100),tterm(100)
10  INTEGER last(100),iprint(100)
11  LOGICAL first
12  SAVE first,tlast,tterm,rap,tk0,last,iprint,nl
13  DATA first /.true./
14  DATA last,iprint /200*0/
15 *
16 *
17 * Open unit #10 the first time.
18  IF (first) THEN
19  OPEN (unit=10,status='NEW',form='FORMATTED',file='HIARCH')
20  first = .false.
21 *
22 * Print cluster scaling parameters at start of the run.
23  IF (nmerg.EQ.0) THEN
24  WRITE (10,1) rbar, bodym*zmbar, body1*zmbar, tscale,
25  & nbin0, nzero
26  1 FORMAT (/,6x,'MODEL: RBAR =',f5.1,' <M> =',f6.2,
27  & ' M1 =',f6.1,' TSCALE =',f6.2,
28  & ' NB =',i4,' N0 =',i6,//)
29  END IF
30 *
31 * Define the saved variables in case of COMMON dump during merger.
32  tk0(1) = 1.0
33  rap = 1.0
34  nl = 0
35  tlast(1) = 0.0
36  tterm(1) = 0.0
37  WRITE (10,2)
38  2 FORMAT (' TIME A A1 E1 PMIN',
39  & ' P1/P0 Q PCR MB Q1 I NAME',
40  & ' KC',/)
41  WRITE (10,3)
42  3 FORMAT (' r/Rc ',
43  & ' Rc GAM IT N0 N1 KS NAME',/)
44  END IF
45 *
46 * Set c.m. and component indices.
47  i = n + ipair
48  i1 = 2*ipair - 1
49  i2 = i1 + 1
50  ttot = time + toff
51 *
52 * Skip printing termination at line #2 until after first new merger.
53  IF (iphase.EQ.7.AND.nl.EQ.0) THEN
54  go to 30
55  END IF
56 *
57 * Determine index of current merger (may have been set previously).
58  il = 0
59  DO 4 l = 1,nl
60  IF (name(i1).EQ.last(l)) THEN
61  il = l
62  END IF
63  4 CONTINUE
64 *
65 * Increase membership and set identifier = 0 temporarily for new case.
66  IF (il.EQ.0) THEN
67  nl = nl + 1
68  il = nl
69  last(il) = 0
70  END IF
71 *
72 * Decide whether new merger or termination should be printed.
73  IF (iphase.EQ.6) THEN
74  IF (name(i1).EQ.last(il)) THEN
75  tlast(il) = ttot
76 * Avoid repeated output of same system if terminated recently (< TCR).
77  IF (ttot - tterm(il).LT.tcr) THEN
78  iprint(il) = 0
79  go to 30
80  ELSE
81  iprint(il) = 1
82  END IF
83  ELSE
84 * Save identity and merge time and set print indicator for line #2.
85  last(il) = name(i1)
86  tlast(il) = ttot
87  tterm(il) = 0.0d0
88  iprint(il) = 1
89  END IF
90  ELSE
91  tterm(il) = ttot
92 * Skip output of line #2 if line #1 was not printed.
93  IF (iprint(il).EQ.0) THEN
94  go to 30
95  END IF
96  iprint(il) = 0
97  END IF
98 *
99 * Form semi-major axis and eccentricity (skip hyperbolic orbits).
100  semi = -0.5*body(i)/h(ipair)
101  IF (semi.LT.0.0) go to 30
102  ecc2 = (1.0d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
103  ecc = sqrt(ecc2)
104 * Skip highly eccentric fly-by's which are not formally stable.
105  IF (ecc.GT.0.99) go to 30
106  tk = twopi*semi*sqrt(semi/body(i))
107 *
108 * Decide between new and terminated merger.
109  IF (iphase.EQ.6) THEN
110 * Ensure that unperturbed binary is resolved.
111  IF (list(1,i1).EQ.0.OR.x(1,i1).EQ.x(1,i2)) THEN
112  CALL resolv(ipair,1)
113  END IF
114  q = body(i)/body(jcomp)
115  q1 = max(body(i2)/body(i1),body(i1)/body(i2))
116  rap = semi*(1.0d0 + ecc)
117  rij2 = 0.0
118  vij2 = 0.0
119  rdot = 0.0
120  a12 = 0.0
121  a22 = 0.0
122  a1a2 = 0.0
123  ri2 = 0.0
124  vi2 = 0.0
125  rvi = 0.0
126  DO 5 k = 1,3
127  rij2 = rij2 + (x(k,i) - x(k,jcomp))**2
128  vij2 = vij2 + (xdot(k,i) - xdot(k,jcomp))**2
129  rdot = rdot + (x(k,i) - x(k,jcomp))*
130  & (xdot(k,i) - xdot(k,jcomp))
131  k1 = k + 1
132  IF (k1.GT.3) k1 = 1
133  k2 = k1 + 1
134  IF (k2.GT.3) k2 = 1
135  a1(k) = (x(k1,i1) - x(k1,i2))*(xdot(k2,i1) - xdot(k2,i2))
136  & - (x(k2,i1) - x(k2,i2))*(xdot(k1,i1) - xdot(k1,i2))
137  a2(k) = (x(k1,jcomp) - x(k1,i))*
138  & (xdot(k2,jcomp) - xdot(k2,i))
139  & - (x(k2,jcomp) - x(k2,i))*
140  & (xdot(k1,jcomp) - xdot(k1,i))
141  a12 = a12 + a1(k)**2
142  a22 = a22 + a2(k)**2
143  a1a2 = a1a2 + a1(k)*a2(k)
144 * Form relative vectors and scalars for inner binary.
145  xrel(k) = x(k,i1) - x(k,i2)
146  vrel(k) = xdot(k,i1) - xdot(k,i2)
147  ri2 = ri2 + xrel(k)**2
148  vi2 = vi2 + vrel(k)**2
149  rvi = rvi + xrel(k)*vrel(k)
150  5 CONTINUE
151 *
152 * Evaluate orbital parameters for outer orbit (skip hyperbolic case).
153  rij = sqrt(rij2)
154  zmb = body(i) + body(jcomp)
155  semi1 = 2.0/rij - vij2/zmb
156  semi1 = 1.0/semi1
157  IF (semi1.LT.0.0) go to 30
158  ecc1 = sqrt((1.0 - rij/semi1)**2 + rdot**2/(semi1*zmb))
159  pmin = semi1*(1.0d0 - ecc1)
160  tk1 = twopi*semi1*sqrt(semi1/zmb)
161  pmin = min(pmin/rap,999.0d0)
162  tk1 = min(tk1/tk,99999.0d0)
163 * Determine inclination (8 bins of 22.5 degrees).
164  fac = a1a2/sqrt(a12*a22)
165  fac = acos(fac)
166 * II = 1 + FAC*360.0/(TWOPI*22.5)
167  ideg = 360.0*fac/twopi
168 *
169 * Construct Runge-Lenz vector (Heggie & Rasio, 1995, IAU174, Eq.(5)).
170  ei2 = 0.0
171  DO 6 k = 1,3
172  ei(k) = (vi2*xrel(k) - rvi*vrel(k))/body(i) -
173  & xrel(k)/sqrt(ri2)
174  ei2 = ei2 + ei(k)**2
175  6 CONTINUE
176 *
177 * Define unit vectors for inner eccentricity and angular momenta.
178  cosj = 0.0
179  sjsg = 0.0
180  DO 8 k = 1,3
181  ei(k) = ei(k)/sqrt(ei2)
182  hi(k) = a1(k)/sqrt(a12)
183  ho(k) = a2(k)/sqrt(a22)
184  cosj = cosj + hi(k)*ho(k)
185  sjsg = sjsg + ei(k)*ho(k)
186  8 CONTINUE
187 *
188 * Form the expressions A & Z.
189  a = cosj*sqrt(1.0 - ei2)
190  z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
191 *
192 * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
193  z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
194  emax = one6*(z + 1.0 - 4.0*a**2 + sqrt(z2))
195  emax = sqrt(emax)
196  kcm = kstar(i)
197  IF (name(i).LT.0) kcm = -10
198  newhi = newhi + 1
199 *
200  WRITE (10,10) ttot, semi, semi1, ecc1, pmin, tk1, q,
201  & pcrit/semi, body(i)/bodym, q1, ideg,
202  & name(i1), name(i2), name(jcomp), kcm
203  10 FORMAT (/,' #1',f8.1,1p,2e10.2,0p,f6.2,f7.2,f8.1,3f6.2,f5.1,
204  & i5,2i5,i6,i4)
205 *
206  WRITE (30,10) ttot, semi, semi1, ecc1, pmin, tk1, q,
207  & pcrit/semi, body(i)/bodym, q1, ideg,
208  & name(i1), name(i2), name(jcomp), kcm
209 *
210  rbig = max(radius(i1),radius(i2))
211  pmin = semi*(1.0 - ecc)
212  pmin2 = semi*(1.0 - emax)
213  IF (kz(19).GE.3) THEN
214  r1 = pmin/rbig
215  r2 = pmin2/rbig
216  ELSE
217  r1 = 0.0
218  r2 = 0.0
219  END IF
220  zi = fac
221  em = sqrt(sin(zi)**2 + ei2*cos(zi)**2)
222  WRITE (30,11) sqrt(ei2), emax, kstar(i1), kstar(i2),kstar(i),
223  & semi, a, z, pmin2, r1, r2, em
224  11 FORMAT (' INNER: E EMAX K* SEMI A Z PM2 R1 R2 EM ',
225  & 2f8.4,3i3,1p,4e10.2,0p,2f7.1,f7.3)
226  CALL flush(30)
227 * Save parameters for termination diagnostics.
228  tk0(il) = tk
229 *
230  ELSE IF (iphase.EQ.7) THEN
231  pmin = semi*(1.0d0 - ecc)
232  pmin = min(pmin/rap,999.0d0)
233  IF (tk0(il).LE.0.0d0) tk0(il) = tk
234  tk1 = min(tk/tk0(il),99999.0d0)
235  nk = (ttot - tlast(il))/tk0(il)
236  nk1 = (ttot - tlast(il))/tk
237  nk = min(nk,99999999)
238  nk1 = min(nk1,9999)
239  nk = max(nk,0)
240  ri = 0.0
241  DO 15 k = 1,3
242  ri = ri + (x(k,i) - rdens(k))**2
243  15 CONTINUE
244  rr = rc/rscale
245 *
246 * Define type index (0, 1 or KSTAR; note: #I1 is inner c.m.).
247  it = 0
248  IF (time.GT.min(tev(i1),tev(i2),tev(i))) it = 1
249  IF (kstar(i1).NE.0.AND.it.EQ.0) it = kstar(i1)
250 *
251  WRITE (10,20) ttot, sqrt(ri)/rc, semi, ecc, pmin, tk1,
252  & rr, gamma(ipair), it, nk, nk1, npairs, name(i2)
253  20 FORMAT (' #2',f8.1,1p,2e10.2,0p,f6.2,f7.2,f8.1,2f6.2,i4,i10,
254  & 3i6)
255  CALL flush(10)
256  WRITE (30,20) ttot, sqrt(ri)/rc, semi, ecc, pmin, tk1,
257  & rr, gamma(ipair), it, nk, nk1, npairs, name(i2)
258  END IF
259 *
260 * Restrict memory of previous cases by removing the oldest one.
261  30 IF (nl.GE.100) THEN
262  lt = 1
263 * Skip any current mergers (defined by TTERM = 0.0).
264  32 IF (tterm(lt).LE.0.0d0) THEN
265  lt = lt + 1
266  go to 32
267  END IF
268  DO 35 l = lt,nl-1
269  tlast(l) = tlast(l+1)
270  tterm(l) = tterm(l+1)
271  tk0(l) = tk0(l+1)
272  last(l) = last(l+1)
273  iprint(l) = iprint(l+1)
274  35 CONTINUE
275  nl = nl - 1
276  END IF
277 *
278  RETURN
279 *
280  END