Nbody6
ksres2.f
Go to the documentation of this file.
1  SUBROUTINE ksres2(J,J1,J2,RIJ2)
2 *
3 *
4 * Coordinates & velocities of KS pair.
5 * ------------------------------------
6 *
7  include 'common6.h'
8  common/slow0/ range,islow(10)
9  REAL*8 ui(4),rdot(3),v(4),a1(3,4)
10 *
11 *
12 * Skip merged ghost binary (defined by zero total mass).
13  IF (body(n+j).LE.0.0d0) go to 40
14 *
15 * Resolve components of pair #J at curent time.
16  j2 = j + j
17  j1 = j2 - 1
18  a2 = 1.0/r(j)
19  a3 = a2*(time - t0(j1))
20  IF (kslow(j).GT.1) THEN
21  imod = kslow(j)
22  a3 = a3/float(islow(imod))
23  END IF
24 *
25 * Decide appropriate order for interpolation & prediction.
26  IF (rij2.GT.625.0*r(j)**2) THEN
27 * Convert physical interval to regularized time (second order only).
28  dtu = (1.0 - 0.5d0*tdot2(j)*a2*a3)*a3
29  IF (abs(dtu).GT.dtau(j)) dtu = dtau(j)
30 *
31 * Predict regularized coordinates of distant pair to second order.
32  DO 10 k = 1,4
33  ui(k) = (fu(k,j)*dtu + udot(k,j))*dtu + u0(k,j)
34  v(k) = (3.0*fudot(k,j)*dtu + 2.0*fu(k,j))*dtu + udot(k,j)
35  10 CONTINUE
36  ELSE
37 * Expand regularized time interval to third order.
38  a4 = 3.0d0*tdot2(j)**2*a2 - tdot3(j)
39  dtu = ((one6*a4*a3 - 0.5d0*tdot2(j))*a2*a3 + 1.0)*a3
40 * Apply safety test near small pericentre or for unperturbed motion.
41  IF (dtu.GT.dtau(j)) dtu = 0.8*dtau(j)
42  IF (list(1,j1).EQ.0) dtu = 0.0
43 *
44 * Predict regularized coordinates to third order.
45  dtu1 = dtu/24.0d0
46  dtu2 = one6*dtu
47  DO 20 k = 1,4
48  ui(k) = (((fudot2(k,j)*dtu1 + fudot(k,j))*dtu +
49  & fu(k,j))*dtu + udot(k,j))*dtu + u0(k,j)
50  v(k) = ((fudot2(k,j)*dtu2 + 3.0d0*fudot(k,j))*dtu +
51  & 2.0d0*fu(k,j))*dtu + udot(k,j)
52  20 CONTINUE
53  END IF
54 *
55 * Employ KS transformation.
56  q1 = ui(1)**2 - ui(2)**2 - ui(3)**2 + ui(4)**2
57  q2 = ui(1)*ui(2) - ui(3)*ui(4)
58  q3 = ui(1)*ui(3) + ui(2)*ui(4)
59  q2 = q2 + q2
60  q3 = q3 + q3
61  j3 = n + j
62  a2 = body(j2)/body(j3)
63 *
64 * Set global coordinates of regularized components.
65  x(1,j1) = x(1,j3) + a2*q1
66  x(2,j1) = x(2,j3) + a2*q2
67  x(3,j1) = x(3,j3) + a2*q3
68  x(1,j2) = x(1,j1) - q1
69  x(2,j2) = x(2,j1) - q2
70  x(3,j2) = x(3,j1) - q3
71 *
72 * Set current transformation matrix and two-body separation.
73  CALL matrix(ui,a1)
74  ri = ui(1)**2 + ui(2)**2 + ui(3)**2 + ui(4)**2
75  rinv = 2.0/ri
76 *
77 * Obtain relative velocities from KS transformation.
78  DO 30 l = 1,3
79  rdot(l) = 0.0d0
80  DO 25 k = 1,4
81  rdot(l) = rdot(l) + a1(l,k)*v(k)*rinv
82  25 CONTINUE
83  30 CONTINUE
84 *
85 * Set global velocities of KS components.
86  xdot(1,j1) = xdot(1,j3) + a2*rdot(1)
87  xdot(2,j1) = xdot(2,j3) + a2*rdot(2)
88  xdot(3,j1) = xdot(3,j3) + a2*rdot(3)
89  xdot(1,j2) = xdot(1,j1) - rdot(1)
90  xdot(2,j2) = xdot(2,j1) - rdot(2)
91  xdot(3,j2) = xdot(3,j1) - rdot(3)
92 *
93  40 RETURN
94 *
95  END