Nbody6
 All Files Functions Variables
unpert.f
Go to the documentation of this file.
1  SUBROUTINE unpert(IPAIR)
2 *
3 *
4 * Unperturbed two-body motion.
5 * ----------------------------
6 *
7  include 'common6.h'
8  REAL*8 ui(4),uidot(4)
9  SAVE tcall
10  DATA tcall /0.0d0/
11 *
12 *
13 * Set first component & c.m. index and semi-major axis.
14  i1 = 2*ipair - 1
15  i = n + ipair
16  semi = -0.5d0*body(i)/h(ipair)
17 *
18 * Add elapsed unperturbed Kepler periods and update the time T0(I1).
19  tk = twopi*semi*sqrt(semi/body(i))
20  IF (time - t0(i1).LT.2.0e+09*tk) THEN
21  k = nint((time - t0(i1))/tk)
22  ELSE
23  k = 0
24  nprect = nprect + 1
25  END IF
26  t0(i1) = time
27 *
28 * Reset unperturbed counter if > 2*10**9 and update the frequency.
29  IF (nksper.GT.2000000000.OR.nksper.LT.0) THEN
30  nksper = 0
31  nprect = nprect + 1
32  END IF
33  nksper = nksper + k
34 *
35 * Include case of tidal dissipation or partial reflection (suppressed).
36  IF (tdot2(ipair).GE.0.0d0) THEN
37 * IF (KZ(25).GT.0.AND.ABS(TDOT2(IPAIR).GT.1.0E-10) THEN
38 * DT = 0.0
39 * GO TO 10
40 * END IF
41 * Ensure perturbation check at least once every c.m. step (KSTAR < 11).
42  IF (kz(27).GT.0) THEN
43  IF (kstar(i).LT.11) THEN
44  IF (step(i1).LT.tk) THEN
45  IF (step(i1).LT.0.0001*step(i)) go to 9
46  END IF
47  IF (time - t0(i).GT.2.0*step(i1)) THEN
48  dt = min(3.0d0*step(i1),step(i))
49  kpert = 0
50  jclose = 0
51  go to 20
52  END IF
53  END IF
54  END IF
55  END IF
56 *
57 * Evaluate interval for unperturbed motion (GAMMA < GMIN).
58  9 kpert = 1
59  CALL tpert(ipair,gmin,dt)
60 *
61 * Restore KS indicator and re-initialize if interval < period.
62  IF (dt.LT.tk) THEN
63 * Form perturber list and restart KS motion if required.
64  CALL kslist(ipair)
65  IF (list(1,i1).GT.0) THEN
66 * Transform to apocentre variables in case of tidal dissipation.
67  IF (r(ipair).LT.semi.AND.semi.LT.5.0*rmin) THEN
68  np = list(1,i1)
69  list(1,i1) = 0
70 * Do not allow backwards integration on switch from unperturbed motion.
71  CALL ksperi(ipair)
72  list(1,i1) = np
73  CALL ksapo(ipair)
74  END IF
75  kslow(ipair) = 1
76  CALL resolv(ipair,1)
77  CALL kspoly(ipair,1)
78 * Include differential correction (experimental 10/02).
79  jlist(1) = i
80  nnb = list(1,i1)
81  DO 12 l = 1,nnb
82  jpert(l) = list(l+1,i1)
83  12 CONTINUE
84  CALL nbpot(1,nnb,pot1)
85  jlist(1) = i1
86  jlist(2) = i1 + 1
87  CALL nbpot(2,nnb,pot2)
88 * WRITE (6,13) NAME(I1), POT1-POT2, GAMMA(IPAIR)
89 * 13 FORMAT (' CORRECT NAM DPHI G ',I6,1P,2E10.2)
90  emerge = emerge + (pot2 - pot1)
91  be(3) = be(3) + (pot2 - pot1)
92  go to 30
93  END IF
94  ir = 1
95  go to 28
96  END IF
97 *
98 * Check for tidal dissipation (peri & apocentre included; skip merger).
99  20 IF (kz(27).GT.0.AND.
100  & (semi.LT.rmin.OR.tdot2(ipair).GT.0.0d0).AND.
101  & (kstar(i).EQ.0.OR.kstar(i).EQ.-1).AND.
102  & (name(i).GT.0)) THEN
103 *
104 * Compare pericentre and effective capture distance.
105  ecc2 = (1.0 - r(ipair)/semi)**2 +
106  & tdot2(ipair)**2/(body(i)*semi)
107  ecc = sqrt(ecc2)
108  rp = semi*(1.0d0 - ecc)
109  rt = 4.0*max(radius(i1),radius(i1+1))
110 *
111 * Specify circularization index (skip SPIRAL but include CHAOS).
112  icirc = 0
113  IF (kz(27).EQ.1.AND.ttot.GT.tcall) THEN
114  IF (rp.LT.rt) icirc = 1
115  tcall = ttot + 0.01
116  ELSE IF (rp.LT.2.5*rt.AND.kstar(i).EQ.0.AND.
117  & ttot.GT.tcall) THEN
118 * Perform occasional checks of TC (NB! small periods).
119  CALL tcirc(rp,ecc,i1,i1+1,icirc,tc)
120  tcall = ttot + 0.01
121  ELSE IF (kstar(i).EQ.-1) THEN
122  icirc = 1
123  END IF
124 *
125 * See whether dissipation is active (A*(1 - E) < 4*MAX(R) & E > 0.002).
126 * IF (RP.LT.0.99*RT.AND.ECC.GT.0.002) THEN
127  IF (icirc.GT.0.AND.ecc.GT.0.002) THEN
128 * Transform to pericentre if R > A (KSAPO works both ways with PI/2).
129  IF (r(ipair).GT.semi) THEN
130  CALL ksapo(ipair)
131  END IF
132 * Implement energy loss at pericentre (exit on collision).
133  CALL kstide(ipair,rp)
134  IF (iphase.LT.0) go to 30
135  semi = -0.5d0*body(i)/h(ipair)
136  tk = twopi*semi*sqrt(semi/body(i))
137 * Assign one Kepler period for active state and impose TDOT2 > 0.
138  IF (r(ipair).LT.0.99*rt.AND.kstar(i).NE.10) THEN
139  step(i1) = tk
140  tdot2(ipair) = 1.0e-20
141  ELSE
142 * Define apocentre and set 1/2 period for inactive or synchronous case.
143  CALL ksapo(ipair)
144  step(i1) = 0.5*tk
145  ecc2 = (1.0 - r(ipair)/semi)**2 +
146  & tdot2(ipair)**2/(body(i)*semi)
147  ecc = sqrt(ecc2)
148 * Note next unperturbed check at apocentre since T'' < 0 in KSAPO.
149  IF (abs(rt - rp)/rt.GT.0.1.AND.kz(27).LE.1) THEN
150  WRITE(6,25) ecc, semi, r(ipair), rp, radius(i1)
151  25 FORMAT (' INACTIVE PHASE E A R RP R* ',
152  & f7.3,1p,4e10.2)
153 * Include safety condition for extremely small SEMI (i.e. tiny period).
154  IF (semi.LT.0.001*rmin) step(i1) = step(i)
155  END IF
156  IF (ecc.LT.0.002.AND.semi.LT.0.01*rmin) THEN
157  kstar(i) = 10
158  tev(i) = time
159  END IF
160  END IF
161  go to 30
162  END IF
163  END IF
164 *
165 * Check for Roche overflow from synchronous orbit (but KSTAR < 13).
166  ir = 0
167  IF (kstar(i).EQ.13.AND.min(tev(i1),tev(i1+1)).LT.time + dt) THEN
168  ir = -1
169  END IF
170  IF (kstar(i).GT.13) ir = 1
171  IF (kz(34).GT.0.AND.(kstar(i).GE.10.AND.kstar(i).LE.12)) THEN
172  tm = min(tev(i1),tev(i1+1),tev(i))
173  IF (tm.LT.time) THEN
174  CALL trflow(ipair,dtr)
175 * Exit if ROCHE is indicated (change 16/08/2006).
176 * IF (DTR.LT.STEP(I1)) THEN
177 * CALL ROCHE(IPAIR)
178 * Exit on coalescence (KPERT > 0 would cause trouble).
179 * IF (IPHASE.LT.0) GO TO 30
180 * ELSE
181 * IR = 1
182 * END IF
183  ir = 1
184  IF (dtr.LT.step(i1)) go to 30
185 *
186  IF (dtr.GT.tk) THEN
187  dt = min(dtr,dt)
188  kpert = 2
189  END IF
190  ELSE
191  ir = 1
192  END IF
193  END IF
194 *
195 * Perform general two-body collision test.
196  IF (kz(19).GE.3.AND.name(i).GT.0) THEN
197  ri = 0.0
198  DO 26 k = 1,4
199  ui(k) = u0(k,ipair)
200  uidot(k) = udot(k,ipair)
201  ri = ri + ui(k)**2
202  26 CONTINUE
203  CALL peri(ui,uidot,ri,body(i1),body(i1+1),qperi)
204  IF (qperi.LT.0.75*(radius(i1) + radius(i1+1))) THEN
205 * Obtain KS variables at pericentre before coalescence to one body.
206  CALL ksperi(ipair)
207  kspair = ipair
208 * Set indicator for skipping ECOLL updating in COAL.
209  iqcoll = -2
210  CALL cmbody(qperi,2)
211  IF (iphase.LT.0) go to 30
212  END IF
213  END IF
214 *
215 * Specify a conservative number of unperturbed orbits from TPERT.
216  IF (kpert.GT.0.AND.dt.GT.0.0) THEN
217 * Ensure that next look-up time is not exceeded (option #27).
218  IF (kz(27).GT.0.AND.kpert.GT.1) THEN
219  tm = min(tev(i1),tev(i1+1))
220  IF (tm - t0(i1).GT.0.0.AND.tm - t0(i1).LT.0.5*dt) THEN
221  nwarn = nwarn + 1
222  IF (nwarn.LT.1000) THEN
223  WRITE (25,27) ipair, kstar(i1), kstar(i1+1),
224  & kstar(i), tm-t0(i1), step(i1)
225  27 FORMAT (' UNPERT WARNING! KS K* TEV-T0 STEP1 ',
226  & 4i4,1p,2e10.2)
227  CALL flush(25)
228  END IF
229  dt = min(2.0*(tm - t0(i1)),dt)
230  dt = max(dt,tk,stepx)
231  END IF
232  END IF
233 * Adopt c.m. step instead if integer argument exceeds 10**9.
234  IF (dt.LT.2.0e+09*tk) THEN
235  k = 1 + int(0.5d0*dt/tk)
236 * Restrict Kepler period to c.m. step (case of very wide orbit).
237  step(i1) = float(k)*min(tk,step(i))
238  ELSE
239  step(i1) = step(i)
240  END IF
241 * Include optional treatment for spiralling of chaotic binary orbit.
242  IF (kz(27).GT.1.AND.kstar(i).EQ.-2) THEN
243 * Ensure pericentre position after possible perturbed motion.
244  IF (r(ipair).GT.semi) THEN
245  CALL ksrect(ipair)
246 * Reduce eccentric anomaly by pi for inward motion.
247  IF (tdot2(ipair).LT.0.0d0) THEN
248  CALL ksapo(ipair)
249  END IF
250 * Transform from outward motion (anomaly < pi) to exact pericentre.
251  CALL ksperi(ipair)
252  END IF
253  CALL spiral(ipair)
254  IF (iphase.LT.0) go to 30
255  END IF
256  END IF
257 *
258 * Check merger condition before continuing (skip small Roche steps).
259  28 IF (kz(15).GT.0.AND.ir.GE.0.AND.time+step(i1).GT.tblock) THEN
260  IF (step(i).LT.dtmin) THEN
261  CALL impact(i)
262  ELSE IF (jclose.GT.0.AND.step(i).LT.10.0*dtmin) THEN
263  CALL histab(ipair,jclose,pmin,rstab)
264  IF (rstab.LT.pmin) THEN
265  CALL impact(i)
266  END IF
267 * Include a more generous test for massive quadruples.
268  ELSE IF (jclose.GT.n) THEN
269  fac = 2.0*(body(i) + body(jclose))/bodym
270  IF (step(i).LT.fac*dtmin) THEN
271  CALL histab(ipair,jclose,pmin,rstab)
272  IF (rstab.LT.pmin) THEN
273  CALL impact(i)
274  END IF
275  END IF
276  END IF
277  END IF
278 *
279 * Check collision criterion for special case.
280  IF (kz(27).EQ.-1.AND.kz(13).LT.0) THEN
281  ecc2 = (1.0 - r(ipair)/semi)**2 +
282  & tdot2(ipair)**2/(body(i)*semi)
283  ecc = sqrt(ecc2)
284  qperi = semi*(1.0 - ecc)
285  rfac = 2.0
286  i2 = i1 + 1
287  IF (qperi.LT.rfac*max(radius(i1),radius(i2))) THEN
288  j1 = i1
289  IF (radius(i2).GT.radius(i1)) j1 = i2
290  fac = 0.5*body(i)/body(j1)
291 * Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
292  rcoll = 1.7*fac**0.3333*radius(j1)
293  IF (qperi.LT.rcoll) THEN
294  CALL touch(ipair,i1,i2,rcoll)
295  END IF
296  END IF
297  END IF
298 *
299  30 RETURN
300 *
301  END