Nbody6
 All Files Functions Variables
chrect.f
Go to the documentation of this file.
1  SUBROUTINE chrect(IPAIR,DMR)
2 *
3 *
4 * Rectification of chaotic orbit.
5 * -------------------------------
6 *
7  include 'common6.h'
8  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
9  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
10  & rp(ntmax),es(ntmax),cm(2,ntmax),iosc(ntmax),
11  & namec(ntmax)
12  REAL*8 ww(6),w(4),wg(2),qg(2),wscale(2),qscale(2)
13  CHARACTER*8 which1
14  LOGICAL sleep
15  DATA ww /2.119,3.113,8.175,3.742,4.953,9.413/
16 *
17 *
18 * Define c.m. & KS indices and search current names for chaos index.
19  i = n + ipair
20  i1 = 2*ipair - 1
21  i2 = i1 + 1
22  ic = 0
23  DO 1 k = 1,nchaos
24  IF (namec(k).EQ.name(i)) ic = k
25  1 CONTINUE
26 *
27 * Include case of chain chaos without identified NAMEC.
28  IF (ic.EQ.0.AND.nchaos.GT.0) THEN
29  WRITE (6,2) nchaos, ipair, kstar(i), name(i1), name(i2),
30  & list(1,i1), step(i1), step(i)
31  2 FORMAT (' WARNING! CHRECT NCH KS K* NAM NP DT1 DTI ',
32  & 3i4,2i6,i4,1p,2e10.2)
33 * See whether wrong component name (+ NZERO) saved as NAMEC in CHAOS2.
34  nam2 = 0
35  DO 6 l = 1,2
36  nam1 = ksave(2*l) - name(i2)
37  IF (ksave(2*l-1).LT.0.AND.nam1.EQ.name(i1)) THEN
38  nam2 = nzero + name(i2)
39  END IF
40  6 CONTINUE
41  ic = nchaos
42  namc = namec(ic)
43 * Check identification for correct value (two K*=-2 are possible).
44  DO 8 k = 1,nchaos
45  IF (namec(k).EQ.nam2) ic = k
46  8 CONTINUE
47  namec(ic) = name(i)
48  IF (nam2.EQ.namc.OR.ic.LT.nchaos) THEN
49  WRITE (6,9) ic, nchaos, nam2, namc, name(i)
50  9 FORMAT (' CHRECT RESTORE IC NCH NM2 NMC NMI ',2i4,3i8)
51  END IF
52  ELSE IF (nchaos.EQ.0) THEN
53  WRITE (6,3) nchaos, ipair, kstar(i), name(i)
54  3 FORMAT (' CHRECT RESTORE NCH KS K* NAM ',3i4,i6)
55 * Restore case of former merger with KSTARM < 0 to chaos table.
56  nchaos = 1
57  ic = 1
58  namec(nchaos) = name(i)
59  END IF
60 *
61 * Save variables for diagnostic output.
62  time0 = tosc(ic)
63  es0 = es(ic)
64  tc = 0.0
65  idis = kstar(i)
66 *
67 * Obtain current values of KS variables in case of spiral.
68  IF (kstar(i).EQ.-2) THEN
69 * Skip on call from RESET/MDOT with arbitrary phase (only updating).
70  IF (abs(tdot2(ipair)).GT.1.0d-12.AND.dmr.GE.0.0) THEN
71 * Rectify KS variables in order to obtain correct pericentre.
72  CALL ksrect(ipair)
73 * Reduce eccentric anomaly by pi for inward motion.
74  IF (tdot2(ipair).LT.0.0d0) THEN
75  CALL ksapo(ipair)
76  END IF
77 * Transform from outward motion to exact pericentre.
78  CALL ksperi(ipair)
79  END IF
80 *
81 * Form current two-body elements.
82  semi = -0.5*body(i)/h(ipair)
83  ecc2 = (1.0 - r(ipair)/semi)**2 +
84  & tdot2(ipair)**2/(body(i)*semi)
85  ecc = sqrt(ecc2)
86 *
87 * Update periastron and eccentricity (ECC modulation or mass loss).
88  qperi = semi*(1.0d0 - ecc)
89  rp(ic) = qperi
90  es(ic) = ecc
91 *
92 * Update orbital parameters after merger, mass loss or radius change.
93  IF (dmr.GE.0.0d0) THEN
94  kstar(i) = -kstar(i)
95  CALL spiral(ipair)
96  IF (iphase.LT.0.OR.kstar(i).GT.0) go to 30
97  END IF
98 *
99 * Check circularization time for spiral after radius/orbit expansion.
100  icirc = -1
101  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
102  tc = max(tc,0.01d0)
103 * Restrict lookup time to TC/2 for exit from possible merger.
104  tev(i1) = min(tev(i1),time + 0.5*tc/tstar)
105 * TEV(I2) = TEV(I1)
106  END IF
107 *
108 * Form (new) semi-major axis, eccentricity & periastron distance.
109  semi = -0.5*body(i)/h(ipair)
110  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
111  ecc = sqrt(ecc2)
112  qperi = semi*(1.0d0 - ecc)
113 *
114 * Set latest values of periastron & eccentricity and update time.
115  rp(ic) = qperi
116  es(ic) = ecc
117  tosc(ic) = time
118 *
119 * Include check for SLEEP after recent WD/NS formation.
120  sleep = .false.
121  km = max(kstar(i1),kstar(i2))
122  IF (km.GE.10) THEN
123 * Determine index of most recent degenerate object formation.
124  IF (kstar(i2).LT.10) THEN
125  j1 = i1
126  ELSE IF (kstar(i1).LT.10) THEN
127  j1 = i2
128  ELSE
129  j1 = i1
130  IF (epoch(i2).GT.epoch(i1)) j1 = i2
131  END IF
132  IF (time - tev0(j1).LT.2.0*stepx.AND.tc.GT.3000.0) THEN
133  sleep = .true.
134  END IF
135  END IF
136 *
137 * Check circularization time after chain regularization.
138  IF (iphase.EQ.8) THEN
139  icirc = -1
140  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
141  IF (tc.GT.3000.0) THEN
142  WRITE (6,4) name(i1), tc
143  4 FORMAT (' CHAIN SLEEP: NM TC ',i6,1p,e10.2)
144  sleep = .true.
145  END IF
146  END IF
147 *
148 * Re-initialize all chaos parameters after expansion or check SLEEP.
149  IF (dmr.GT.0.01.OR.kstar(i).EQ.-1.OR.sleep) THEN
150  IF (dmr.GT.0.01.AND.tc.GT.100.0.AND.km.LT.5) THEN
151  WRITE (6,5) name(i1), name(i2), kstar(i1), kstar(i2),
152  & kstar(i), tphys, radius(i1), radius(i2),
153  & qperi, semi, ecc, es0, body(i)*zmbar, tc
154  5 FORMAT (' CHRECT: NAM K* TP R* QP A E E0 M TC ',
155  & 2i6,3i4,f8.1,1p,4e10.2,0p,3f7.3,f7.1)
156  END IF
157 *
158 * Reset spiral indicator for degenerate component and long t_{circ}.
159  IF (sleep) THEN
160  kstar(i) = 0
161  nslp = nslp + 1
162  tk = semi*sqrt(semi/body(i))
163  tb = yrs*tk
164  xp = (time - time0)/(twopi*tk)
165  qps = semi*(1.0 - ecc)/max(radius(i1),radius(i2))
166  WRITE (6,10) ttot, name(i1), name(i2), kstar(i1),
167  & kstar(i2), ecc, es0, qps, semi, tc, tb, xp
168  10 FORMAT (' SLEEP SPIRAL T NM K* E E0 QP/S A TC TB DTK ',
169  & f9.2,2i6,2i4,2f8.4,1p,5e9.1)
170 * Update chaos variables at end of routine SPIRAL (argument < 0).
171  ii = -i
172  CALL spiral(ii)
173  go to 30
174  END IF
175 *
176 * Save eccentricity, binding energy & J0 and initialize EDEC & IOSC.
177  zmu = body(i1)*body(i2)/body(i)
178  cj = zmu*sqrt(body(i))
179  eb0(ic) = zmu*h(ipair)
180  zj0(ic) = cj*sqrt(qperi*(1.0 + ecc))
181  edec(ic) = 0.0
182  iosc(ic) = 1
183  xn = 0.0
184  DO 12 k = 1,2
185  ik = i1 + k - 1
186  IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5.OR.
187  & kstar(ik).EQ.6.OR.kstar(ik).EQ.9) THEN
188  CALL giant(ipair,ik,wg,qg,wscale,qscale,xn,ql)
189  w(k) = wg(1)
190  ELSE
191  ip = 3
192  IF (kstar(ik).EQ.0) ip = 1
193  w(k) = ww(ip)
194  END IF
195  12 CONTINUE
196 *
197 * Set new chaos boundary parameters (ECRIT, AR & BR).
198  CALL chaos0(qperi,ecc,eb0(ic),zj0(ic),body(i1),body(i2),
199  & radius(i1),radius(i2),w,ecrit(ic),ar(ic),br(ic),idis)
200 *
201  rcoll = radius(i1) + radius(i2)
202  IF (idis.EQ.-1.AND.kstar(i).EQ.-1) THEN
203  iosc(ic) = 2
204  WRITE (6,15) ttot, ipair, name(i1), name(i2), kstar(i1),
205  & kstar(i2), radius(i1), radius(i2), qperi,
206  & semi, es0, ecc
207  15 FORMAT (' CHAOS => SPIRAL T KS NAM K* R* QP A E0 E ',
208  & f9.2,i4,2i6,2i4,1p,4e10.2,0p,2f7.3)
209 * Activate spiral indicator and save time, pericentre & eccentricity.
210  kstar(i) = -2
211  tosc(ic) = time
212  rp(ic) = qperi
213  es(ic) = ecc
214  nsp = nsp + 1
215  IF (kz(8).GT.3) THEN
216  CALL binev(ipair)
217  END IF
218  go to 30
219  END IF
220 *
221 * Combine the two stars inelastically in case of chaos disruption.
222  IF (idis.GT.0.AND.qperi.LT.rcoll) THEN
223  r1 = max(radius(i1),radius(i2))
224  which1 = ' CHAOS '
225  IF (kstar(i).EQ.-2) which1 = ' SPIRAL '
226  WRITE (6,20) which1, ipair, name(i1), name(i2),
227  & kstar(i1), kstar(i2), r1, r(ipair), qperi,
228  & semi, ecc, es0, xn
229  20 FORMAT (' DISRUPTED',a8,' KS NM K* R* R QP A E E0 n ',
230  & i4,2i6,2i4,1p,4e10.2,0p,3f7.3)
231 * Update chaos variables at end of routine SPIRAL (argument < 0).
232  iqcoll = 1
233  IF (kstar(i).EQ.-2) iqcoll = 2
234  ii = -i
235  CALL spiral(ii)
236  kstar(i) = 0
237  CALL xvpred(i,0)
238  kspair = ipair
239  CALL cmbody(r(ipair),2)
240  dmr = -1.0
241  go to 30
242  END IF
243  ELSE IF (kstar(i).EQ.-2) THEN
244  i1 = 2*ipair - 1
245  i2 = i1 + 1
246  IF (qperi.LT.2.0*max(radius(i1),radius(i2))) THEN
247 * Determine indices for primary & secondary star (donor & accretor).
248  j1 = i1
249  IF (radius(i2).GT.radius(i1)) j1 = i2
250 * Define mass ratio and evaluate Roche radius for the primary.
251  q0 = body(j1)/(body(i) - body(j1))
252  q1 = q0**0.3333
253  q2 = q1**2
254  rl1 = 0.49*q2/(0.6*q2 + log(1.0d0 + q1))*semi
255 * Check Roche radius but skip RESET call (no second EMERGE correction).
256  IF (radius(j1).GT.rl1.AND.iphase.NE.7) THEN
257  WRITE (6,25) name(i1), name(i2), kstar(i1), kstar(i2),
258  & ecc, es0, rcoll, rl1, qperi, semi, tc
259  25 FORMAT (' DISRUPTED SPIRAL NM K* E E0 RC RL QP A',
260  & ' TC ',2i6,2i4,2f7.3,1p,5e10.2)
261 * Obtain KS variables at pericentre and enforce collision or CE.
262  CALL ksperi(ipair)
263  kstar(i) = 0
264  kspair = ipair
265  iqcoll = 2
266  CALL cmbody(qperi,2)
267  dmr = -1.0
268 * Note that same result achieved by TERM SPIRAL in case IPHASE = 7.
269  END IF
270  END IF
271  END IF
272 *
273  30 RETURN
274 *
275  END