Nbody6
 All Files Functions Variables
bindat.f
Go to the documentation of this file.
1  SUBROUTINE bindat
2 *
3 *
4 * Binary data bank.
5 * -----------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflag(mmax)
11  REAL*4 eb(kmax),ecc(kmax),rcm(kmax),ecm(kmax),pb(kmax),as(30)
12  REAL*8 xx(3,3),vv(3,3)
13  LOGICAL first
14  SAVE first
15  DATA first /.true./
16 *
17 *
18 * Form binding energy and central distance for each KS pair.
19  zmbin = 0.0
20  DO 10 jpair = 1,npairs
21  j2 = 2*jpair
22  j1 = j2 - 1
23  icm = n + jpair
24  zmbin = zmbin + body(icm)
25  bodycm = body(icm)
26 * Determine merger & ghost index for negative c.m. name (skip ghost).
27  IF (name(icm).LT.0.AND.body(icm).GT.0.0) THEN
28  CALL findj(j1,j,imerge)
29 * Employ actual masses and two-body distance for energy & eccentricity.
30  bodycm = cm(1,imerge) + cm(2,imerge)
31  eb(jpair) = cm(1,imerge)*cm(2,imerge)*hm(imerge)/bodycm
32  semi = -0.5*bodycm/hm(imerge)
33  rj = sqrt(xrel(1,imerge)**2 + xrel(2,imerge)**2 +
34  & xrel(3,imerge)**2)
35 * Assume that merged binary is near apo or peri (hence ignore TDOT2).
36  ecc2 = (1.0 - rj/semi)**2
37 * Include separate diagnostics for the hierarchy (inner comps J1 & J).
38  semi1 = -0.5*body(icm)/h(jpair)
39  ecc1 = (1.0 - r(jpair)/semi1)**2 +
40  & tdot2(jpair)**2/(body(icm)*semi1)
41  e0 = sqrt(ecc2)
42  e1 = sqrt(ecc1)
43  IF (j.LT.0) j = j1
44  rm = semi*(1.0 - e0)/max(radius(j1),radius(j),1.0d-20)
45  rm = min(rm,99.9d0)
46  p0 = days*semi*sqrt(abs(semi)/bodycm)
47  p1 = days*semi1*sqrt(abs(semi1)/body(icm))
48  DO 2 k = 1,3
49  xx(k,1) = xrel(k,imerge)
50  xx(k,2) = 0.0
51  xx(k,3) = x(k,j2)
52  vv(k,1) = vrel(k,imerge)
53  vv(k,2) = 0.0
54  vv(k,3) = xdot(k,j2)
55  2 CONTINUE
56  CALL inclin(xx,vv,x(1,icm),xdot(1,icm),alph)
57  pcr = stability(cm(1,imerge),cm(2,imerge),body(icm),e0,e1,
58  & alph)*semi
59  pm = semi1*(1.0 - e1)/pcr
60  WRITE (84,3) ttot, name(j1), name(j), kstar(j1), kstar(j),
61  & kstarm(imerge), e0, e1, pm, rm, p0, p1, semi1
62  3 FORMAT (' BINDAT: T NM K* E0 E1 PM/PC PM0/R* P0 P1 A1',
63  & f8.1,2i5,3i4,2f7.3,2f6.1,1p,3e9.1)
64  CALL flush(84)
65  ELSE IF (body(j1).GT.0.0d0) THEN
66 * Form binding energy and eccentricity for standard case.
67  eb(jpair) = body(j1)*body(j2)*h(jpair)/
68  & (body(j1) + body(j2))
69  semi = -0.5*body(icm)/h(jpair)
70  ecc2 = (1.0 - r(jpair)/semi)**2 +
71  & tdot2(jpair)**2/(body(icm)*semi)
72  ELSE
73  im = 0
74 * Search merger table to identify corresponding index of c.m. name.
75  DO 5 k = 1,nmerge
76  IF (nameg(k).EQ.name(icm)) THEN
77  im = k
78  END IF
79  5 CONTINUE
80  bodyj1 = cm(3,im)
81  bodyj2 = cm(4,im)
82  bodycm = bodyj1 + bodyj2
83  bodycm = max(bodycm,1.0d-10)
84  eb(jpair) = bodyj1*bodyj2*h(jpair)/bodycm
85  semi = -0.5*bodycm/h(jpair)
86  ecc2 = (1.0 - semi/r(jpair))**2
87  END IF
88  ecc(jpair) = sqrt(ecc2)
89  eb(jpair) = max(eb(jpair),-9.99999)
90  pb(jpair) = days*semi*sqrt(abs(semi)/bodycm)
91  pb(jpair) = min(pb(jpair),99999.9)
92  IF (semi.LT.0.0) pb(jpair) = 0.0
93 * Obtain binding energy (per unit mass) of c.m. motion.
94  vj2 = xdot(1,icm)**2 + xdot(2,icm)**2 + xdot(3,icm)**2
95  IF (body(icm).EQ.0.0d0) vj2 = 0.0
96  potj = 0.0
97  DO 9 j = ifirst,ntot
98  IF (j.EQ.icm) go to 9
99  rij2 = (x(1,icm) - x(1,j))**2 + (x(2,icm) - x(2,j))**2 +
100  & (x(3,icm) - x(3,j))**2
101  potj = potj + body(j)/sqrt(rij2)
102  9 CONTINUE
103  ecm(jpair) = 0.5*vj2 - potj
104 * Check for external tidal field (note that HT includes mass).
105  IF (kz(14).GT.0) THEN
106  CALL xtrnlv(icm,icm)
107  ecm(jpair) = ecm(jpair) + ht/(body(icm) + 1.0e-20)
108  END IF
109  rcm(jpair) = sqrt((x(1,icm) - rdens(1))**2 +
110  & (x(2,icm) - rdens(2))**2 +
111  & (x(3,icm) - rdens(3))**2)
112  rcm(jpair) = min(rcm(jpair),99.9)
113  10 CONTINUE
114 *
115 * Copy relevant binary diagnostics to single precision.
116  as(1) = time + toff
117  as(2) = rscale
118  as(3) = rtide
119  as(4) = rc
120  as(5) = tphys
121  as(6) = -1.5*(tidal(1)*zmass**2)**0.3333
122  as(7) = 0.0
123  DO 20 k = 1,10
124  as(k+7) = e(k)
125  20 CONTINUE
126  as(18) = sbcoll
127  as(19) = bbcoll
128  as(20) = zkin
129  as(21) = pot
130  as(22) = ebin0
131  as(23) = ebin
132  as(24) = esub
133  as(25) = emerge
134  as(26) = be(3)
135  as(27) = zmass
136  as(28) = zmbin
137  as(29) = chcoll
138  as(30) = ecoll
139 *
140 * Write formatted data bank on unit 9.
141  IF (first) THEN
142  OPEN (unit=9,status='NEW',form='FORMATTED',file='OUT9')
143  first = .false.
144  END IF
145 *
146  WRITE (9,30) npairs, model, nrun, n, nc, nmerge, (as(k),k=1,7)
147  30 FORMAT (3i4,i6,2i4,2x,f7.1,2f7.2,f7.3,f8.1,2f9.4)
148  WRITE (9,35) (as(k),k=8,17)
149  35 FORMAT (10f11.6)
150  WRITE (9,40) (as(k),k=18,30)
151  40 FORMAT (13f10.5)
152 *
153  DO 50 jpair = 1,npairs
154  j1 = 2*jpair - 1
155  j2 = 2*jpair
156  kcm = kstar(n+jpair)
157  IF (name(n+jpair).LT.0) THEN
158  kcm = -10
159  END IF
160  WRITE (9,45) eb(jpair), ecc(jpair), ecm(jpair), rcm(jpair),
161  & body(j1)*zmbar, body(j2)*zmbar, pb(jpair),
162  & name(j1), name(j2), kstar(j1), kstar(j2), kcm
163  45 FORMAT (f8.5,f7.3,f7.2,f6.2,2f5.1,f8.1,2i6,3i4)
164  50 CONTINUE
165  CALL flush(9)
166 *
167  RETURN
168 *
169  END