Nbody6
nbtide.f
Go to the documentation of this file.
1  SUBROUTINE nbtide(I,J,RJMIN2)
2 *
3 *
4 * Close two-body interaction.
5 * ---------------------------
6 *
7  include 'common6.h'
8  REAL*8 xi(3),xj(3),vi(3),vj(3),vcm(3),vr(3)
9 *
10 *
11 * Copy coordinates of close bodies and predict velocity for body #J.
12  rdot = 0.0
13  vrel2 = 0.0
14  dt = time - t0(j)
15  DO 10 k = 1,3
16  xi(k) = x(k,i)
17  xj(k) = x(k,j)
18  vi(k) = x0dot(k,i)
19  vj(k) = (3.0*fdot(k,j)*dt + 2.0*f(k,j))*dt + x0dot(k,j)
20  rdot = rdot + (xi(k) - xj(k))*(vi(k) - vj(k))
21  vrel2 = vrel2 + (vi(k) - vj(k))**2
22  10 CONTINUE
23 *
24 * Only consider approaching bodies.
25  IF (rdot.LE.0.0) go to 100
26 *
27 * Predict coordinates & velocities at beginning of step for body #I.
28  rdot0 = 0.0
29  dti = time - t0(i) - step(i)
30  dtj = time - t0(j) - step(i)
31  DO 20 k = 1,3
32  xi(k) = ((fdot(k,i)*dti + f(k,i))*dti + x0dot(k,i))*dti +
33  & x0(k,i)
34  xj(k) = ((fdot(k,j)*dtj + f(k,j))*dtj + x0dot(k,j))*dtj +
35  & x0(k,j)
36  vi(k) = (3.0*fdot(k,i)*dti + 2.0*f(k,i))*dti + x0dot(k,i)
37  vj(k) = (3.0*fdot(k,j)*dtj + 2.0*f(k,j))*dtj + x0dot(k,j)
38  rdot0 = rdot0 + (xi(k) - xj(k))*(vi(k) - vj(k))
39  20 CONTINUE
40 *
41 * Check pericentre condition (old radial velocity < 0).
42  IF (rdot0.GT.0.0) go to 100
43 *
44  zm = body(i) + body(j)
45  erel = 0.5*vrel2 - zm/sqrt(rjmin2)
46 *
47 * Specify the energy loss (experimental).
48  dh = 0.1*vrel2/2.0
49  hnew = erel - dh
50  semi0 = -0.5*zm/erel
51  semi = -0.5*zm/hnew
52  WRITE (6,50) name(i), name(j), semi0, semi, sqrt(rjmin2), dh
53  50 FORMAT (5x,'NBTIDE: NAMES A0 A RIJ DH ',2i5,2f10.5,f8.4,f8.3)
54 *
55 * Skip energy loss treatment unless final orbit is significantly bound.
56  IF (semi.LT.0.0.OR.semi.GT.4.0*rmin) go to 100
57 *
58 * Predict coordinates & velocity for body #J to highest order.
59  CALL xvpred(j,0)
60 *
61 * Set current velocities and form c.m. velocity.
62  DO 30 k = 1,3
63  vi(k) = xdot(k,i)
64  vj(k) = xdot(k,j)
65  vcm(k) = (body(i)*xdot(k,i) + body(j)*xdot(k,j))/zm
66  30 CONTINUE
67 *
68 * Introduce velocity change for each component and initialize X0DOT.
69  fac = sqrt((vrel2 - 2.0d0*dh)/vrel2)
70  DO 40 k = 1,3
71  vr(k) = fac*(vi(k) - vj(k))
72  xdot(k,i) = vcm(k) + body(j)*vr(k)/zm
73  xdot(k,j) = vcm(k) - body(i)*vr(k)/zm
74  x0dot(k,i) = vi(k)
75  x0dot(k,j) = vj(k)
76  40 CONTINUE
77 *
78  de = body(i)*body(j)*dh/zm
79 * Increase event counter and update total energy loss.
80  ndiss = ndiss + 1
81  ecoll = ecoll + de
82 *
83 * Set components and phase indicator for new KS regularization.
84  icomp = min(i,j)
85  jcomp = max(i,j)
86  iphase = 1
87 *
88  WRITE (6,60) name(icomp), name(jcomp), semi, be(3) + de , de
89  60 FORMAT (' TIDAL CAPTURE NM A E DE ',2i5,f8.4,2f11.6)
90 *
91  100 RETURN
92 *
93  END