Nbody6
 All Files Functions Variables
brake.f
Go to the documentation of this file.
1  SUBROUTINE brake(IPAIR,DSEP,ECC1)
2 *
3 *
4 * Orbital changes (GR, mass loss and/or tides).
5 * ---------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 m1,m2
9  SAVE ndiag
10  DATA ndiag /0/
11 *
12 *
13 * Set indices of KS components & c.m.
14  i1 = 2*ipair - 1
15  i2 = i1 + 1
16  i = n + ipair
17  dsep = dsep/su
18 *
19 * Set mass and radius of each star in solar units.
20  m1 = body(i1)*zmbar
21  m2 = body(i2)*zmbar
22  r1 = radius(i1)*su
23  r2 = radius(i2)*su
24 *
25 * Obtain current semi-major axis & eccentricity and define new semi.
26  semi = -0.5d0*body(i)/h(ipair)
27  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
28  ecc = sqrt(ecc2)
29  semi1 = semi - dsep
30 *
31 * Include circularization check (otherwise unperturbed KS missed).
32  IF (kz(27).GT.0.AND.ecc.GT.0.01) THEN
33  qperi = semi*(1.0 - ecc)
34  icirc = -1
35  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
36  IF (tc.LT.100.0) THEN
37  kx = max(kstar(i1),kstar(i2))
38  WRITE (6,3) icirc, kx, ecc, qperi*su, tc, max(r1,r2)
39  3 FORMAT (' TCIRC CHECK IC K* E QP TC R* ',2i5,1p,3e10.2)
40  CALL flush(6)
41  END IF
42  IF (icirc.GT.0.AND.tc.LT.100.0) THEN
43  CALL kstide(ipair,qperi)
44  go to 50
45  END IF
46  END IF
47 *
48 * Check for circularized case.
49  IF(kstar(i).LT.10.AND.ecc1.LT.0.002)THEN
50  ecc1 = 0.001d0
51  kstar(i) = 10
52  ENDIF
53 *
54 * Check collision/coalescence (IQCOLL=-2 skips ECOLL updating in COAL).
55  rp = semi1*su*(1.d0 - ecc1)
56 *
57 * Include GR coalescence criterion for compact objects.
58  ksx = max(kstar(i1),kstar(i2))
59  IF (ksx.GE.13.AND.kz(28).GT.0) THEN
60  rcoal = 6.0*body(i)/clight**2
61  ELSE
62  rcoal = r1 + r2
63  END IF
64  IF (kz(28).GT.1.AND.dsep.GT.0.0) THEN
65  WRITE (6,5) ttot, name(i1), ksx, ecc, semi*su, dsep/semi
66  5 FORMAT (' BRAKE T NAM KX* E A DA/A ',
67  & f8.2,i6,i4,f9.5,1p,2e10.2)
68  END IF
69 *
70 * Check collision condition for stars or degenerate objects.
71  IF (rp.LT.rcoal) THEN
72  CALL ksperi(ipair)
73  iqcoll = -2
74  kspair = ipair
75  CALL cmbody(r(ipair),2)
76  goto 50
77  ENDIF
78  IF(dsep.EQ.0.d0) goto 50
79 *
80 * Include safety test on new semi-major axis.
81  rchck = min(radius(i1),radius(i2))
82  IF(semi1.LT.rchck) semi1 = rchck
83 *
84 * Transform to pericentre (R > A & unperturbed).
85 * IF(R(IPAIR).GT.SEMI.AND.LIST(1,I1).EQ.0)THEN
86 * CALL KSAPO(IPAIR)
87 * ENDIF
88 *
89 * Form square of regularized velocity.
90  v20 = 0.0
91  DO 10 k = 1,4
92  v20 = v20 + udot(k,ipair)**2
93  10 CONTINUE
94 *
95 * Update binding energy and collision energy.
96  hi = h(ipair)
97  h(ipair) = -0.5*body(i)/semi1
98  zmu = body(i1)*body(i2)/body(i)
99  ecoll = ecoll + zmu*(hi - h(ipair))
100  egrav = egrav + zmu*(hi - h(ipair))
101 *
102 * Distinguish between update of eccentric binaries and standard case.
103  IF(abs(ecc-ecc1).GT.tiny)THEN
104 * Change KS variables at original ecc and modify ECC at H = const.
105  CALL expand(ipair,semi)
106  CALL ksperi(ipair)
107  CALL deform(ipair,ecc,ecc1)
108  ELSE
109 * Specify KS coordinate & velocity scaling factors at arbitrary point.
110  c2 = sqrt(semi1/semi)
111 * V2 = 0.5*(BODY(I) + H(IPAIR)*SEMI1*(1.0 - ECC))
112  v2 = 0.5*(body(i) + h(ipair)*r(ipair)*(semi1/semi))
113  c1 = sqrt(v2/v20)
114 *
115 * Re-scale KS variables to new energy with constant eccentricity.
116  r(ipair) = 0.0d0
117 * TDOT2(IPAIR) = 0.0D0
118  DO 20 k = 1,4
119  u(k,ipair) = c2*u(k,ipair)
120  udot(k,ipair) = c1*udot(k,ipair)
121  u0(k,ipair) = u(k,ipair)
122  r(ipair) = r(ipair) + u(k,ipair)**2
123 * TDOT2(IPAIR) = TDOT2(IPAIR) + 2.0*U(K,IPAIR)*UDOT(K,IPAIR)
124  20 CONTINUE
125  END IF
126 *
127 * Transform back to apocentre for standard unperturbed motion.
128 * IF (LIST(1,I1).EQ.0) THEN
129 * CALL KSAPO(IPAIR)
130 * END IF
131 *
132 * Include new initialization for perturbed orbit.
133  IF (list(1,i1).GT.0) THEN
134  imod = kslow(ipair)
135  CALL kspoly(ipair,imod)
136  ENDIF
137 *
138 * Include some diagnostic output.
139  IF(kstar(i).EQ.13)THEN
140  ndiag = ndiag + 1
141  IF(ndiag.LT.100.OR.mod(ndiag,100).EQ.100)THEN
142  rcoll = radius(i1) + radius(i2)
143  WRITE (6,25) ttot, ipair, m1, m2, r1, r2, r(ipair),
144  & rcoll
145  25 FORMAT (' BRAKE T KS M12 R12 R RCOLL ',
146  & f10.4,i4,2f6.2,2f7.3,1p,2e10.2)
147  ENDIF
148  ENDIF
149 *
150  50 RETURN
151 *
152  END