Nbody6
 All Files Functions Variables
brake2.f
Go to the documentation of this file.
1  SUBROUTINE brake2(IPAIR,ITERM)
2 *
3 *
4 * Gravitational radiation of hierarchical 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  REAL*8 m0,aj,r1,lum1,mc,rcc,menv,renv,k2
13  REAL*8 tscls(20),lums(10),gb(10),tm,tn
14  REAL*8 m1,m2,rl,q
15  EXTERNAL rl,rtmsf
16 *
17 *
18  i = n + ipair
19  i1 = 2*ipair - 1
20  im = 0
21  DO 1 k = 1,nmerge
22  IF (namem(k).EQ.name(i)) im = k
23  1 CONTINUE
24 * Find location of the secondary (i.e. ghost).
25  CALL findj(i,i2,im)
26  m1 = cm(1,im)*zmbar
27  m2 = cm(2,im)*zmbar
28 *
29 * Quit if active ROCHE or mass-losing star.
30  iterm = 1
31  IF(kstarm(im).EQ.11.OR.kstarm(im).EQ.13)THEN
32  go to 30
33  END IF
34  IF (m1.GE.10.0.OR.(kstar(i1).GT.1.AND.kstar(i1).LT.10))THEN
35  go to 30
36  END IF
37  IF (m2.GE.10.0.OR.(kstar(i2).GT.1.AND.kstar(i2).LT.10))THEN
38  go to 30
39  END IF
40 *
41 * Skip braking if semi-major axis exceeds 10 solar radii.
42  iterm = 0
43  zmb = cm(1,im) + cm(2,im)
44  semi = -0.5*zmb/hm(im)
45  IF (semi*su.GT.10.0)THEN
46  tgr = 1.0d+06
47  tbr = 1.0d+06
48  semi1 = semi
49  go to 20
50  END IF
51 *
52 * Obtain current separation, square velocity and t'' = 2*U*U'.
53  rb = 0.0
54  v20 = 0.0
55  td2 = 0.0
56  DO 5 k = 1,4
57  rb = rb + um(k,im)**2
58  v20 = v20 + umdot(k,im)**2
59  td2 = td2 + 2.0*um(k,im)*umdot(k,im)
60  5 CONTINUE
61 *
62 * Evaluate inner eccentricity and outer period.
63  ecc2 = (1.0 - rb/semi)**2 + td2**2/(zmb*semi)
64  ecc = sqrt(ecc2)
65  semi0 = -0.5*body(i)/h(ipair)
66 *
67 * Identify the most magnetically active star (#J2).
68  IF (radius(i1).GE.radius(i2)) THEN
69  IF ((kstar(i1).EQ.1.AND.m1.GT.1.25).OR.
70  & (kstar(i1).EQ.0.AND.m1.LT.0.35)) THEN
71  j2 = i2
72  ELSE
73  j2 = i1
74  END IF
75  ELSE
76  IF ((kstar(i2).EQ.1.AND.m2.GT.1.25).OR.
77  & (kstar(i2).EQ.0.AND.m2.LT.0.35)) THEN
78  j2 = i1
79  ELSE
80  j2 = i2
81  END IF
82  END IF
83 *
84 * Set mass and radius of star #J2 in solar units.
85  IF (j2.EQ.i1) m2 = m1
86  r2 = radius(j2)*su
87 *
88 * Define binary period, solar radius & GMS in cgs units.
89  tb = 3.147d+07*yrs*semi*sqrt(semi/zmb)
90  rsun = 6.96d+10
91  gms = gm*zmbar*zmb
92 *
93 * Form time derivative of angular momentum (Regos & Tout 1995).
94  zjdot = 3.8d-30*1.989d+33*m2*rsun**4*r2**3*(twopi/tb)**3
95 *
96 * Determine new semi-major axis from angular momentum relation.
97  zmu = cm(1,im)*cm(2,im)/zmb
98  acm = 3.08d+18*rbar*semi
99  adot = 2.0*sqrt(acm/gms)*zjdot/(1.989d+33*zmbar*zmu)
100 *
101 * Define old primary and evaluate time scale for spin-down (yrs).
102  tbr = acm/(3.147d+07*adot)
103 * Evaluate timescale for 10% change in semi-major axis (N-body units).
104  tbr = 0.1d0*tbr/(1.0d+06*tstar)
105 *
106 * Include gravitational radiation (64*GM**3/5*c**5; Ap.J. 418, 147).
107 * FAC = 64.0*((6.67D-08*1.989D+33)/3.0D+10)**3/(5.0*9.0D+20)
108  agdot = 1.23d+27*cm(1,im)*cm(2,im)*zmb*(zmbar/acm)**3
109  tgr = acm/(3.147d+07*agdot)
110  tgr = 0.1d0*tgr/(1.0d+06*tstar)
111 *
112 * Suppress magnetic braking for massive MS/low-mass or evolved star.
113  IF (((m2.GT.1.25.OR.m2.LT.0.35).AND.kstar(j2).LE.2).OR.
114  & kstar(j2).GE.10) THEN
115  adot = 0.0
116  tbr = 1.0d+06
117  END IF
118 *
119 * Specify time interval.
120  dt = tev(i) - tev0(i)
121  tk = twopi*abs(semi0)*sqrt(abs(semi0)/body(i))
122 *
123 * Combine effects but include possible cutoff of magnetic braking.
124  adot = adot + agdot
125 *
126 * Convert from cgs to scaled units and update semi-major axis.
127  adot = adot/(1.0d+05*vstar)
128  IF (adot*dt.GT.0.1*semi) THEN
129  semi1 = 0.9*semi
130  ELSE
131  semi1 = semi - adot*dt
132  END IF
133 *
134  WRITE (6,71) semi, semi1, adot, dt
135  71 FORMAT (' BRAKE2 !! A A1 AD DT ',1p,4e10.2)
136 * Include safety test on new semi-major axis.
137  IF (abs(semi1).LT.radius(j2).OR.semi1.LT.0.0) THEN
138  semi1 = radius(j2)
139  END IF
140 *
141 * Update binding energy.
142  hi = hm(im)
143  hm(im) = -0.5*zmb/semi1
144 *
145 * Correct EMERGE & ECOLL (consistency; no net effect).
146  decorr = zmu*(hi - hm(im))
147  emerge = emerge - decorr
148  ecoll = ecoll + decorr
149  egrav = egrav + decorr
150 *
151 * Specify KS coordinate & velocity scaling factors at general point.
152  c2 = sqrt(semi1/semi)
153  v2 = 0.5*(zmb + hm(im)*rb*(semi1/semi))
154  c1 = sqrt(v2/v20)
155 *
156 * Re-scale KS variables to new energy with constant eccentricity.
157  rb = 0.0d0
158  DO 10 k = 1,4
159  um(k,im) = c2*um(k,im)
160  umdot(k,im) = c1*umdot(k,im)
161  rb = rb + um(k,im)**2
162  10 CONTINUE
163 *
164 * Rectify the hierarchical KS variables.
165  CALL hirect(im)
166 *
167 * Determine indices for primary & secondary star (donor & accretor).
168  20 j1 = i1
169  j2 = i2
170  q = cm(1,im)/cm(2,im)
171  rl1 = rl(q)*semi1
172 * Evaluate Roche radius for the second star.
173  q = 1.0/q
174  rl2 = rl(q)*semi1
175 *
176 * Update radius for MS stars to current time.
177  tm0 = 100.d0
178  rtms = 0.d0
179  IF(kstar(i1).LE.1)THEN
180  m0 = m1
181  mc = 0.d0
182  aj = time*tstar - epoch(i1)
183  kw = kstar(i1)
184  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
185  CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
186  & r1,lum1,kw,mc,rcc,menv,renv,k2)
187  IF(kw.GT.1)THEN
188  tm0 = 0.d0
189  tev(i1) = min(tev(i1),time)
190  ELSE
191  tm = (tm + epoch(i1))/tstar - time
192  tm0 = min(tm0,tm)
193  radius(i1) = r1/su
194  tev0(i1) = time
195  rtms = rtmsf(m1)/su
196  dtms = tm
197  ENDIF
198  ENDIF
199  IF(kstar(i2).LE.1)THEN
200  m0 = m2
201  mc = 0.d0
202  aj = time*tstar - epoch(i2)
203  kw = kstar(i2)
204  CALL star(kw,m0,m2,tm,tn,tscls,lums,gb,zpars)
205  CALL hrdiag(m0,aj,m2,tm,tn,tscls,lums,gb,zpars,
206  & r1,lum1,kw,mc,rcc,menv,renv,k2)
207  IF(kw.GT.1)THEN
208  tm0 = 0.d0
209  tev(i2) = min(tev(i2),time)
210  ELSE
211  tm = (tm + epoch(i2))/tstar - time
212  tm0 = min(tm0,tm)
213  radius(i2) = r1/su
214  tev0(i2) = time
215  rtms2 = rtmsf(m2)/su
216  IF(rtms2.GT.rtms)THEN
217  rtms = rtms2
218  dtms = tm
219  ENDIF
220  ENDIF
221  ENDIF
222  IF(kstar(i1).GE.10.AND.tm0.GT.0.d0) tev0(i1) = time
223  IF(kstar(i2).GE.10.AND.tm0.GT.0.d0) tev0(i2) = time
224 *
225 * Compare scaled Roche radii when choosing the primary.
226  IF (radius(j1)/rl1.LT.radius(j2)/rl2) THEN
227  j1 = i2
228  j2 = i1
229  rl1 = rl2
230  END IF
231 *
232 * Check for Roche mass transfer at every stage (AM CVn formation).
233  IF (radius(j1).GT.rl1.AND.kz(34).GT.0) THEN
234  WRITE (6,25) name(i1), nameg(im), kstar(i1), kstar(i2),
235  & semi1*su, rl1*su, radius(j1)*su
236  25 FORMAT (' BRAKE2 TERM NM K* A RL R* ',2i6,2i4,3f7.3)
237  iterm = 1
238 *
239 * Check for termination of MS stage.
240  ELSEIF (tm0.EQ.0.d0) THEN
241  tev(i) = time + step(i1)
242  iterm = 1
243 *
244 * Increase TEV for standard case.
245  ELSE
246  dt = min(tgr,tbr,tm0)
247  IF(rtms.GT.rl1)THEN
248 * Estimate time of Roche-lobe filling for MS star.
249  dt1 = dtms*(rl1 - radius(j1))/(rtms - radius(j1))
250  dt = min(dt1,dt)
251  ENDIF
252  tev0(i) = tev(i)
253  tev(i) = tev(i) + dt
254  tev(i1) = tev(i) + step(i1)
255  tev(i2) = tev(i1)
256  END IF
257 *
258  30 RETURN
259 *
260  END