Nbody6
 All Files Functions Variables
degen.f
Go to the documentation of this file.
1  SUBROUTINE degen(I1,I2,ICASE)
2 *
3 *
4 * Binary output for degenerate stars.
5 * -----------------------------------
6 *
7  include 'common6.h'
8  REAL*8 eb(kmax),semi(kmax),ecc(kmax),rcm(kmax),ecm(kmax)
9  CHARACTER*8 which1
10  LOGICAL first,first2,first3
11  SAVE first,first2,first3
12  DATA first,first2,first3 /.true.,.true.,.true./
13 *
14 *
15 * Skip output for start and end of merger.
16  IF (iphase.EQ.6.OR.iphase.EQ.7) go to 50
17  IF (icase.EQ.7) go to 40
18 *
19 * See whether KS binaries contain any degenerate stars.
20  nb = 0
21  DO 1 ipair = i1,i2
22  j2 = 2*ipair
23  j1 = j2 - 1
24  IF (kstar(j1).GT.9.OR.kstar(j2).GT.9) THEN
25  nb = nb + 1
26  END IF
27  1 CONTINUE
28 *
29 * Open unit #4 the first time.
30  IF (nb.GT.0.AND.first) THEN
31  OPEN (unit=4,status='NEW',form='FORMATTED',file='DEGEN')
32  first = .false.
33 *
34 * Print cluster scaling parameters at start of the run.
35  WRITE (4,2) rbar, bodym*zmbar, body1*zmbar, tscale,
36  & nbin0, nzero
37  2 FORMAT (/,6x,'MODEL: RBAR =',f5.1,' <M> =',f6.2,
38  & ' M1 =',f6.1,' TSCALE =',f6.2,
39  & ' NB =',i4,' N0 =',i6,//)
40 *
41  WRITE (4,3)
42  3 FORMAT (' # TPHYS A E Rp/R* P r',
43  & ' M1 M2 K* NAME',/)
44  END IF
45 *
46 * Form binding energy and central distance for degenerate stars.
47  jpair = 0
48  iprint = 0
49  DO 20 ipair = i1,i2
50  j2 = 2*ipair
51  j1 = j2 - 1
52  IF (kstar(j1).LT.10.AND.kstar(j2).LT.10) go to 20
53  jpair = jpair + 1
54  icm = n + ipair
55 * Avoid division by zero for merged or synchronous ghost binary.
56  IF (body(j1).GT.0.0) THEN
57  eb(jpair) = body(j1)*body(j2)*h(ipair)/
58  & (body(j1) + body(j2))
59  semi(jpair) = -0.5d0*body(icm)/h(ipair)
60  ecc2 = (1.d0 - r(ipair)/semi(jpair))**2 +
61  & tdot2(ipair)**2/(body(icm)*semi(jpair))
62  ecc(jpair) = sqrt(ecc2)
63 * Set zero eccentricity after common envelope stage (still large R).
64  IF (r(ipair).GT.2.0*semi(jpair)) THEN
65  ecc(jpair) = 0.d0
66  END IF
67  eb(jpair) = max(eb(jpair),-9.99999d0)
68  ELSE
69  eb(jpair) = 0.d0
70  semi(jpair) = r(ipair)
71  ecc(jpair) = 0.d0
72  END IF
73  rcm(jpair) = sqrt((x(1,icm) - rdens(1))**2 +
74  & (x(2,icm) - rdens(2))**2 +
75  & (x(3,icm) - rdens(3))**2)
76 * Obtain binding energy (per unit mass) of c.m. motion.
77  vj2 = xdot(1,icm)**2 + xdot(2,icm)**2 +
78  & xdot(3,icm)**2
79 * POTJ = 0.0
80 * DO 5 J = IFIRST,NTOT
81 * IF (J.EQ.ICM) GO TO 5
82 * RIJ2 = (X(1,ICM) - X(1,J))**2 + (X(2,ICM) - X(2,J))**2 +
83 * & (X(3,ICM) - X(3,J))**2
84 * POTJ = POTJ + BODY(J)/SQRT(RIJ2)
85 * 5 CONTINUE
86 * IF (TIME.EQ.TADJ) THEN
87 * POTJ = -PHI(ICM)
88 * ELSE
89  potj = 0.d0
90 * END IF
91  ecm(jpair) = 0.5d0*vj2 - potj
92 * Check for external tidal field (NB! already included on GRAPE).
93 * IF (KZ(14).GT.0) THEN
94 * CALL XTRNLV(ICM,ICM)
95 * ECM(JPAIR) = ECM(JPAIR) + HT/(BODY(ICM) + 1.0D-20)
96 * END IF
97  tphys = ttot*tstar
98  tk = semi(jpair)*sqrt(abs(semi(jpair))/(body(icm) + 1.0d-20))
99  tk = days*tk
100  tk = min(tk,999999.9d0)
101  rbig = max(radius(j1),radius(j2))
102  ratio = semi(jpair)*(1.d0 - ecc(jpair))/rbig
103  ratio = min(ratio,99.9d0)
104  IF (semi(jpair).LT.0.0.AND.ratio.GT.5.0) go to 20
105  semi(jpair) = semi(jpair)*rbar*au
106  IF (iprint.EQ.0.AND.icase.EQ.0) THEN
107  WRITE (4,*)
108  END IF
109  WRITE (4,10) icase, tphys, semi(jpair), ecc(jpair), ratio,
110  & tk, rcm(jpair), body(j1)*zmbar, body(j2)*zmbar,
111  & kstar(j1), kstar(j2), kstar(icm),
112  & name(j1), name(j2)
113  10 FORMAT (i2,f8.1,f8.2,f7.3,f6.1,f9.1,f6.2,2f5.1,3i4,2i6)
114  iprint = iprint + 1
115  20 CONTINUE
116 *
117 * Close file at end of main output.
118  IF (iprint.GT.0.AND.icase.EQ.0) THEN
119  CALL flush(4)
120  END IF
121 *
122 * Search for neutron stars at main output.
123  IF (icase.EQ.0) THEN
124  DO 30 j = 1,n
125  IF (kstar(j).EQ.13) THEN
126  IF (first2) THEN
127  OPEN (unit=33,status='NEW',form='FORMATTED',
128  & file='NS')
129  first2 = .false.
130  END IF
131  IF (j.LT.ifirst) THEN
132  which1 = ' BINARY '
133  i = kvec(j) + n
134  vi2 = xdot(1,i)**2 + xdot(2,i)**2 +
135  & xdot(3,i)**2
136  ELSE
137  which1 = ' SINGLE '
138  vi2 = xdot(1,j)**2 + xdot(2,j)**2 +
139  & xdot(3,j)**2
140  END IF
141  WRITE (33,25) which1, j, name(j), ifirst, kstar(j),
142  & tphys, sqrt(vi2)*vstar
143  25 FORMAT (1x,a8,'NS',' J NAM I* K* TPH V ',
144  & 2i6,i5,i4,2f7.1)
145  CALL flush(33)
146  ELSE IF (kstar(j).GT.13) THEN
147  IF (j.LT.ifirst) THEN
148  jcm = n + kvec(j)
149  IF (name(jcm).LT.0) go to 30
150  END IF
151  IF (first3) THEN
152  OPEN (unit=34,status='NEW',form='FORMATTED',
153  & file='BH')
154  first3 = .false.
155  END IF
156  vi2 = xdot(1,j)**2 + xdot(2,j)**2 +
157  & xdot(3,j)**2
158  vi = sqrt(vi2)
159  which1 = ' AIC '
160  IF (kstar(j).EQ.14) which1 = ' BH '
161  IF (kstar(j).EQ.15) THEN
162 * Ensure that errant massless remnant will escape.
163  ri = sqrt(x(1,j)**2 + x(2,j)**2 + x(3,j)**2)
164  DO 26 l = 1,3
165  x0(l,j) = 1000.0*rscale*x(l,j)/ri
166  x(l,j) = x0(l,j)
167  x0dot(l,j) = sqrt(0.004*zmass/rscale)*
168  & xdot(l,j)/vi
169  xdot(l,j) = x0dot(l,j)
170  26 CONTINUE
171  ELSE
172  vi = vi*vstar
173  WRITE (34,28) which1,j,name(j),kstar(j),tphys,vi
174  28 FORMAT (1x,a8,'J NAM K* TPH V ',2i6,i4,2f7.1)
175  CALL flush(34)
176  END IF
177  END IF
178  30 CONTINUE
179  END IF
180 *
181 * Print single neutron stars at escape time.
182  40 IF (icase.EQ.7) THEN
183  j = i1
184  vi2 = xdot(1,j)**2 + xdot(2,j)**2 + xdot(3,j)**2
185  vi = sqrt(vi2)*vstar
186  which1 = ' ESCAPE '
187  WRITE (33,25) which1, j, name(j), ifirst, kstar(j), tphys, vi
188  END IF
189 *
190 * Include counter for doubly degenerate binary.
191  IF (icase.EQ.3.OR.icase.EQ.4) THEN
192  ipair = i1
193  j1 = 2*ipair - 1
194  j2 = j1 + 1
195  IF (kstar(j1).GE.10.AND.kstar(j2).GE.10) THEN
196  ndd = ndd + 1
197  a = -0.5d0*su*body(n+ipair)/h(ipair)
198  WRITE (6,48) ipair, name(j1), name(j2), kstar(j1),
199  & kstar(j2), kstar(n+ipair), r(ipair), a, tk
200  48 FORMAT (' NEW DD KS NM K* R A P ',i4,2i6,3i4,1p,3e10.2)
201  END IF
202  END IF
203 *
204  50 RETURN
205 *
206  END