Nbody6
brake3.f
Go to the documentation of this file.
1  SUBROUTINE brake3(IPAIR,ITERM)
2 *
3 *
4 * Gravitational radiation of hierarchical BH/NS binary.
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),iflagm(mmax)
11  common/slow0/ range,islow(10)
12  parameter(aursun=214.95d0)
13  REAL*8 m1,m2,jorb
14 *
15 *
16 * Define scalars and determine merger index.
17  i = n + ipair
18  i1 = 2*ipair - 1
19  iterm = 0
20  im = 0
21  DO 1 k = 1,nmerge
22  IF (namem(k).EQ.name(i)) im = k
23  1 CONTINUE
24 *
25 * Find location of secondary (i.e. ghost) and form mass factors & SEMI.
26  CALL findj(i,i2,im)
27  IF (i2.LT.0) go to 30
28  IF(kstar(i2).LT.10)THEN
29  iterm = 1
30  go to 30
31  ENDIF
32  m1 = cm(1,im)*zmbar
33  m2 = cm(2,im)*zmbar
34  zmb = cm(1,im) + cm(2,im)
35  zmu = cm(1,im)*cm(2,im)/zmb
36  semi = -0.5*zmb/hm(im)
37 *
38 * Obtain current separation, square velocity and t'' = 2*U*U'.
39  rb = 0.0
40  v20 = 0.0
41  td2 = 0.0
42  DO 5 k = 1,4
43  rb = rb + um(k,im)**2
44  v20 = v20 + umdot(k,im)**2
45  td2 = td2 + 2.0*um(k,im)*umdot(k,im)
46  5 CONTINUE
47 *
48 * Form inner eccentricity and angular momentum.
49  sep = semi*su
50  ecc2 = (1.0 - rb/semi)**2 + td2**2/(zmb*semi)
51  ecc0 = sqrt(ecc2)
52  oorb = twopi*sqrt((m1 + m2)/(sep/aursun)**3)
53  jorb = m1*m2/(m1 + m2)*sqrt(1.d0-ecc2)*sep**2*oorb
54 *
55 * Obtain angular momentum and eccentricity derivatives (SEMI < 10 R_s).
57 *
58 * Allow 2% angular momentum change (expressions adapted from MDOT).
59  IF (abs(djgr).GT.0.d0) THEN
60  dtgr = 0.02d0*jorb/abs(djgr)
61  IF(delet.GT.tiny.AND.ecc0.GT.0.0011d0)THEN
62  dtgr = min(dtgr,0.05d0*ecc0/delet)
63  ENDIF
64  dtgr = max(dtgr,100.d0)
65 * Note that time interval is in units of yr (cf. grrad.f).
66  djorb = djgr*dtgr
67 * Terminate on large change of JORB which occurs at GR coalescence.
68  IF (djorb.GT.0.1*jorb) THEN
69  iterm = 1
70  go to 30
71  END IF
72  jorb = jorb - djorb
73  jorb = max(jorb,1.d0)
74  ecc = max(ecc0 - delet*dtgr,0.001d0)
75 * Convert to N-body units and update both TEV to same time (NB!).
76  dtgr = dtgr/(1.0d+06*tstar)
77  tev0(i1) = time
78  tev0(i2) = time
79  tev(i1) = time + dtgr
80  tev(i2) = time + dtgr
81  ELSE
82  dt = min(tev(i1)-tev0(i1),tev(i2)-tev0(i2))
83  dt = max(dt,0.1d0)
84  tev0(i1) = time
85  tev0(i2) = time
86  tev(i1) = time + dt
87  tev(i2) = time + dt
88  go to 30
89  END IF
90 *
91 * Obtain new semi-major axis (in N-body units) from angular momentum.
92  semi1 = (m1 + m2)*jorb*jorb/
93  & ((m1*m2*twopi)**2*aursun**3*(1.d0-ecc**2))
94  semi1 = semi1/su
95 *
96 * Update inner binding energy.
97  hi = hm(im)
98  hm(im) = -0.5*zmb/semi1
99 *
100 * Correct EMERGE & ECOLL (consistency; no net effect).
101  decorr = zmu*(hi - hm(im))
102  emerge = emerge - decorr
103  ecoll = ecoll + decorr
104  egrav = egrav + decorr
105 *
106 * Specify KS coordinate & velocity scaling factors at general point.
107  c2 = sqrt(semi1/semi)
108  v2 = 0.5*(zmb + hm(im)*rb*(semi1/semi))
109  c1 = sqrt(v2/v20)
110 *
111 * Re-scale KS variables to new energy with constant eccentricity.
112  rb = 0.0d0
113  DO 10 k = 1,4
114  um(k,im) = c2*um(k,im)
115  umdot(k,im) = c1*umdot(k,im)
116  rb = rb + um(k,im)**2
117  10 CONTINUE
118 *
119 * Specify new peri- or apo-centre at turning-point.
120  IF (rb.LT.semi1) THEN
121  rnew = semi1*(1.0 - ecc)
122  efac = (1.0 - ecc)/(1.0 - ecc0)
123  ELSE
124  rnew = semi1*(1.0 + ecc)
125  efac = (1.0 + ecc)/(1.0 + ecc0)
126  END IF
127  c1 = sqrt(efac)
128  v2 = 0.5*(zmb + hm(im)*rnew)
129  IF (v2.LE.0.d0) THEN
130  c2 = 1.0d-06
131  ELSE
132  c2 = sqrt(v2/v20)
133  END IF
134 *
135 * Re-scale KS variables at constant energy to new eccentricity (ECC).
136  rb = 0.0d0
137  DO 20 k = 1,4
138  um(k,im) = c1*um(k,im)
139  umdot(k,im) = c2*umdot(k,im)
140  rb = rb + um(k,im)**2
141  20 CONTINUE
142 *
143 * Rectify the hierarchical KS variables.
144  CALL hirect(im)
145 *
146 * Include GR coalescence criterion for compact objects.
147  ksx = max(kstar(i1),kstar(i2))
148  IF (ksx.GE.13.AND.kz(28).GT.0) THEN
149  rcoal = 6.0*zmb/clight**2
150  ELSE
151  rcoal = 0.0
152  END IF
153 *
154  WRITE (6,25) name(i1), toff+time, ecc, semi1, dtgr, rcoal
155  25 FORMAT (' BRAKE3 NM T E A DTGR RCOAL',i8,f9.2,f8.4,1p,3e10.2)
156 *
157 * Check termination criterion for coalescence.
158  IF (rb.LT.rcoal) THEN
159  iterm = 1
160  END IF
161  30 CONTINUE
162 *
163  RETURN
164 *
165  END