Nbody6
 All Files Functions Variables
kspred.f
Go to the documentation of this file.
1  SUBROUTINE kspred(IPAIR,I1,I,BODYIN,UI,UIDOT,XI,VI)
2 *
3 *
4 * Prediction for KS regularization.
5 * ---------------------------------
6 *
7  include 'common6.h'
8  parameter(one24=1.0/24.0d0,one120=1.0/120.0d0)
9  REAL*8 ui(4),uidot(4),xi(6),vi(6),rdot(3),a1(3,4)
10 *
11 *
12 * Add body #I to the perturber list for prediction.
13  nnb2 = list(1,i1) + 2
14  list(nnb2,i1) = i
15 *
16 * Predict coordinates & velocities of perturbers & c.m. to order FDOT.
17  DO 10 l = 2,nnb2
18  j = list(l,i1)
19  s = time - t0(j)
20  s1 = 1.5*s
21  s2 = 2.0*s
22  x(1,j) = ((fdot(1,j)*s + f(1,j))*s + x0dot(1,j))*s + x0(1,j)
23  x(2,j) = ((fdot(2,j)*s + f(2,j))*s + x0dot(2,j))*s + x0(2,j)
24  x(3,j) = ((fdot(3,j)*s + f(3,j))*s + x0dot(3,j))*s + x0(3,j)
25  xdot(1,j) = (fdot(1,j)*s1 + f(1,j))*s2 + x0dot(1,j)
26  xdot(2,j) = (fdot(2,j)*s1 + f(2,j))*s2 + x0dot(2,j)
27  xdot(3,j) = (fdot(3,j)*s1 + f(3,j))*s2 + x0dot(3,j)
28  10 CONTINUE
29 *
30 * Copy Stumpff coefficients to scalars.
31  z3 = sf(3,ipair)
32  z4 = sf(4,ipair)
33  z5 = sf(5,ipair)
34 *
35 * Predict KS motion to full accuracy including Stumpff functions.
36  dtu = dtau(ipair)
37  DO 20 k = 1,4
38  u2 = fu(k,ipair)
39  u3 = fudot(k,ipair)
40  u4 = one24*fudot2(k,ipair)
41  u5 = one120*fudot3(k,ipair)
42 *
43  ui(k) = ((((u5*z5*dtu + u4*z4)*dtu + u3)*dtu + u2)*dtu +
44  & udot(k,ipair))*dtu + u0(k,ipair)
45  uidot(k) = (((5.0*u5*z4*dtu + 4.0d0*u4*z3)*dtu + 3.0d0*u3)*dtu
46  & + 2.0*u2)*dtu + udot(k,ipair)
47  20 CONTINUE
48 *
49  r(ipair) = ui(1)**2 + ui(2)**2 + ui(3)**2 + ui(4)**2
50 *
51 * Form relative coordinates obtained from explicit KS transformation.
52  q1 = ui(1)**2 - ui(2)**2 - ui(3)**2 + ui(4)**2
53  q2 = ui(1)*ui(2) - ui(3)*ui(4)
54  q3 = ui(1)*ui(3) + ui(2)*ui(4)
55  q2 = q2 + q2
56  q3 = q3 + q3
57 *
58 * Assign global coordinates of regularized components.
59  a2 = body(i1+1)*bodyin
60  xi(1) = x(1,i) + a2*q1
61  xi(2) = x(2,i) + a2*q2
62  xi(3) = x(3,i) + a2*q3
63  xi(4) = xi(1) - q1
64  xi(5) = xi(2) - q2
65  xi(6) = xi(3) - q3
66 *
67 * Set current transformation matrix.
68  CALL matrix(ui,a1)
69 *
70 * Obtain relative velocities from KS transformation.
71  rinv = 2.0d0/r(ipair)
72  DO 30 l = 1,3
73  rdot(l) = 0.0d0
74  DO 25 k = 1,4
75  rdot(l) = rdot(l) + a1(l,k)*uidot(k)
76  25 CONTINUE
77  rdot(l) = rdot(l)*rinv
78  30 CONTINUE
79 *
80 * Set global velocities of KS components.
81  vi(1) = xdot(1,i) + a2*rdot(1)
82  vi(2) = xdot(2,i) + a2*rdot(2)
83  vi(3) = xdot(3,i) + a2*rdot(3)
84  vi(4) = vi(1) - rdot(1)
85  vi(5) = vi(2) - rdot(2)
86  vi(6) = vi(3) - rdot(3)
87 *
88  RETURN
89 *
90  END