Nbody6
 All Files Functions Variables
resolv.f
Go to the documentation of this file.
1  SUBROUTINE resolv(IPAIR,KDUM)
2 *
3 *
4 * Transformation of KS variables.
5 * -------------------------------
6 *
7  include 'common6.h'
8  common/slow0/ range,islow(10)
9  REAL*8 q(3),rdot(3),ui(4),v(4),a1(3,4)
10 *
11 *
12 * Note: KDUM = 1 for prediction of U & UDOT and = 2 for termination.
13  i1 = 2*ipair - 1
14  i2 = 2*ipair
15 *
16 * Ensure appropriate prediction (second pair in binary collision).
17  IF (list(1,i1).EQ.0.OR.t0(i1).EQ.time) THEN
18 * Set current values of regularized coordinates & velocities.
19  DO 1 k = 1,4
20  ui(k) = u0(k,ipair)
21  v(k) = udot(k,ipair)
22  1 CONTINUE
23 *
24 * Copy binding energy for routine ENERGY.
25  ht = h(ipair)
26  go to 4
27  END IF
28 *
29 * Predict U & UDOT to order FUDOT3 using third-order time inversion.
30  a2 = 1.0/r(ipair)
31  a3 = a2*(time - t0(i1))
32 *
33 * See whether the time interval should be modified by KSLOW procedure.
34  IF (kslow(ipair).GT.1) THEN
35  imod = kslow(ipair)
36  a3 = a3/float(islow(imod))
37  END IF
38 *
39 * Expand regularized interval to third order.
40  a4 = 3.0d0*tdot2(ipair)**2*a2 - tdot3(ipair)
41  dtu = ((one6*a4*a3 - 0.5d0*tdot2(ipair))*a2*a3 + 1.0)*a3
42 * Apply safety test near small pericentre or for unperturbed motion.
43  IF (dtu.GT.dtau(ipair)) dtu = 0.8*dtau(ipair)
44  IF (dtu.LT.0.0) dtu = 0.0
45 *
46  dtu1 = 0.2d0*dtu
47  dtu2 = dtu/24.0d0
48  DO 3 k = 1,4
49  ui(k) = ((((fudot3(k,ipair)*dtu1 + fudot2(k,ipair))*dtu2 +
50  & fudot(k,ipair))*dtu + fu(k,ipair))*dtu +
51  & udot(k,ipair))*dtu + u0(k,ipair)
52  v(k) = (((fudot3(k,ipair)*dtu2 + one6*fudot2(k,ipair))*dtu +
53  & 3.0d0*fudot(k,ipair))*dtu + 2.0d0*fu(k,ipair))*dtu +
54  & udot(k,ipair)
55  3 CONTINUE
56 *
57 * Predict current binding energy per unit mass for routine ADJUST.
58  ht = (((hdot4(ipair)*dtu2 + one6*hdot3(ipair))*dtu +
59  & 0.5d0*hdot2(ipair))*dtu + hdot(ipair))*dtu + h(ipair)
60 *
61 * Form current transformation matrix and two-body separation.
62  4 CALL matrix(ui,a1)
63  ri = ui(1)**2 + ui(2)**2 + ui(3)**2 + ui(4)**2
64 *
65 * Obtain relative coordinates & velocities.
66  DO 7 j = 1,3
67  q(j) = 0.0d0
68  rdot(j) = 0.0d0
69  DO 6 k = 1,4
70  q(j) = q(j) + a1(j,k)*ui(k)
71  rdot(j) = rdot(j) + 2.0d0*a1(j,k)*v(k)/ri
72  6 CONTINUE
73  7 CONTINUE
74 *
75  i = n + ipair
76  IF (kdum.EQ.1) go to 12
77 *
78 * Predict improved X & X0DOT for c.m. & components at termination.
79  dt = time - t0(i)
80  a0 = 0.05*dt
81  a2 = 0.25*dt
82  a3 = t0(i) - t0r(i)
83  a4 = 0.5*a3
84  t0(i) = time
85 *
86  DO 10 k = 1,3
87  fk = ((one6*d3r(k,i)*a3 + 0.5*d2r(k,i))*a3 + d1r(k,i))*a3
88  & + d0r(k,i) + d0(k,i)
89  f1dotk = (d3r(k,i)*a4 + d2r(k,i))*a3 + d1r(k,i) + d1(k,i)
90  f2dotk = 0.5*(d3r(k,i)*a3 + d2r(k,i) + d2(k,i))
91  f3dotk = one6*(d3r(k,i) + d3(k,i))
92  x(k,i) = ((((f3dotk*a0 + one12*f2dotk)*dt +
93  & one6*f1dotk)*dt + 0.5*fk)*dt + x0dot(k,i))*dt +
94  & x0(k,i)
95 * Note: Update X0 for call from KSTERM after RESOLV call in MDOT.
96  x0(k,i) = x(k,i)
97  x0dot(k,i) = (((f3dotk*a2 + one3*f2dotk)*dt +
98  & 0.5*f1dotk)*dt + fk)*dt + x0dot(k,i)
99  xdot(k,i) = x0dot(k,i)
100  x0dot(k,i1) = x0dot(k,i) + body(i2)*rdot(k)/body(i)
101  x0dot(k,i2) = x0dot(k,i) - body(i1)*rdot(k)/body(i)
102  10 CONTINUE
103 *
104 * Set coordinates & velocities (XDOT for KDUM = 1) of KS components.
105  12 DO 15 k = 1,3
106  x(k,i1) = x(k,i) + body(i2)*q(k)/body(i)
107  x(k,i2) = x(k,i) - body(i1)*q(k)/body(i)
108  xdot(k,i1) = xdot(k,i) + body(i2)*rdot(k)/body(i)
109  xdot(k,i2) = xdot(k,i) - body(i1)*rdot(k)/body(i)
110  15 CONTINUE
111 *
112 * Check updating of KS variables (KDUM = 3 in routine MDOT).
113  IF (kdum.LT.3) go to 20
114 *
115  DO 18 k = 1,4
116  u(k,ipair) = ui(k)
117  u0(k,ipair) = ui(k)
118  udot(k,ipair) = v(k)
119  18 CONTINUE
120  r(ipair) = ri
121  h(ipair) = ht
122  t0(i1) = time
123  t0(i) = time
124 *
125  20 RETURN
126 *
127  END