Nbody6
 All Files Functions Variables
trans3.f
Go to the documentation of this file.
1  SUBROUTINE trans3(KDUM)
2 *
3 *
4 * Transformation of physical & KS variables.
5 * ------------------------------------------
6 *
7  IMPLICIT REAL*8 (a-h,m,o-z)
8  common/azreg/ time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
9  & xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
10  & nstep3,name3(3),kz15,kz27
11  REAL*8 q1(9),p1(9),q2(9),p2(9)
12 *
13 *
14 * Decide path (initialization, termination, KS -> phys, switching).
15  IF (kdum.EQ.4) go to 4
16  IF (kdum.GT.1) go to 20
17 *
18 * Obtain total kinetic & potential energy.
19  zke = 0.0d0
20  pot = 0.0d0
21  DO 3 i = 1,3
22  zke = zke + m(i)*(xdot3(1,i)**2 + xdot3(2,i)**2 +
23  & xdot3(3,i)**2)
24  DO 2 j = i+1,3
25  pot = pot - m(i)*m(j)/sqrt((x3(1,i) - x3(1,j))**2 +
26  & (x3(2,i) - x3(2,j))**2 +
27  & (x3(3,i) - x3(3,j))**2)
28  2 CONTINUE
29  3 CONTINUE
30 *
31 * Save twice the initial energy for routine DERQP3.
32  energy = zke + 2.0d0*pot
33  IF (kdum.EQ.0) go to 40
34 *
35 * Introduce physical coordinates & momenta.
36  4 DO 6 i = 1,3
37  DO 5 k = 1,3
38  i1 = 3*i + k - 3
39  q1(i1) = x3(k,i)
40  p1(i1) = m(i)*xdot3(k,i)
41  5 CONTINUE
42  6 CONTINUE
43 *
44 * Set mass factors for routine DERQP3.
45  c11 = 0.25d0/m(1) + 0.25d0/m(3)
46  c12 = 0.25d0/m(2) + 0.25d0/m(3)
47  c19 = 2.0d0*m(2)*m(3)
48  c20 = 2.0d0*m(1)*m(3)
49  c24 = 0.25d0/m(3)
50  c25 = 2.0d0*m(1)*m(2)
51 *
52 * Define relative coordinates & absolute momenta (eqn (45)).
53  DO 7 k = 1,3
54  q2(k) = q1(k) - q1(k+6)
55  q2(k+3) = q1(k+3) - q1(k+6)
56  p2(k+3) = p1(k+3)
57  7 CONTINUE
58 *
59 * Expand the variables by relabelling (eqn (46)).
60  DO 8 k = 1,3
61  q1(k) = q2(k)
62  q1(k+4) = q2(k+3)
63  p1(k+4) = p2(k+3)
64  8 CONTINUE
65 *
66 * Initialize the redundant variables (eqn (47)).
67  q1(4) = 0.0d0
68  q1(8) = 0.0d0
69  p1(4) = 0.0d0
70  p1(8) = 0.0d0
71 *
72 * Introduce regularized variables for each KS pair.
73  DO 10 kcomp = 1,2
74  k = 4*(kcomp - 1)
75 *
76 * Form scalar distance from relative separation vector.
77  rk = sqrt(q1(k+1)**2 + q1(k+2)**2 + q1(k+3)**2)
78 *
79 * Perform KS coordinate transformation (eqn (48) or (49)).
80  IF (q1(k+1).GE.0.0) THEN
81  q(k+1) = sqrt(0.5d0*(rk + q1(k+1)))
82  q(k+2) = 0.5d0*q1(k+2)/q(k+1)
83  q(k+3) = 0.5d0*q1(k+3)/q(k+1)
84  q(k+4) = 0.0d0
85  ELSE
86  q(k+2) = sqrt(0.5d0*(rk - q1(k+1)))
87  q(k+1) = 0.5d0*q1(k+2)/q(k+2)
88  q(k+4) = 0.5d0*q1(k+3)/q(k+2)
89  q(k+3) = 0.0d0
90  END IF
91 *
92 * Set regularized momenta (eqn (50)).
93  p(k+1) = 2.0d0*(+q(k+1)*p1(k+1) + q(k+2)*p1(k+2) +
94  & q(k+3)*p1(k+3))
95  p(k+2) = 2.0d0*(-q(k+2)*p1(k+1) + q(k+1)*p1(k+2) +
96  & q(k+4)*p1(k+3))
97  p(k+3) = 2.0d0*(-q(k+3)*p1(k+1) - q(k+4)*p1(k+2) +
98  & q(k+1)*p1(k+3))
99  p(k+4) = 2.0d0*(+q(k+4)*p1(k+1) - q(k+3)*p1(k+2) +
100  & q(k+2)*p1(k+3))
101  10 CONTINUE
102 *
103  go to 40
104 *
105 * Transform each KS pair from regularized to physical variables.
106  20 DO 22 kcomp = 1,2
107  k = 4*(kcomp - 1)
108 *
109 * Obtain relative coordinates (eqn (52)).
110  q1(k+1) = q(k+1)**2 - q(k+2)**2 - q(k+3)**2 + q(k+4)**2
111  q1(k+2) = 2.0d0*(q(k+1)*q(k+2) - q(k+3)*q(k+4))
112  q1(k+3) = 2.0d0*(q(k+1)*q(k+3) + q(k+2)*q(k+4))
113 *
114 * Form product of half Levi-Civita matrix & regularized momenta.
115  p1(k+1) = q(k+1)*p(k+1) - q(k+2)*p(k+2) - q(k+3)*p(k+3) +
116  & q(k+4)*p(k+4)
117  p1(k+2) = q(k+2)*p(k+1) + q(k+1)*p(k+2) - q(k+4)*p(k+3) -
118  & q(k+3)*p(k+4)
119  p1(k+3) = q(k+3)*p(k+1) + q(k+4)*p(k+2) + q(k+1)*p(k+3) +
120  & q(k+2)*p(k+4)
121 *
122 * Evaluate scalar distance.
123  rk = q(k+1)**2 + q(k+2)**2 + q(k+3)**2 + q(k+4)**2
124 *
125 * Set absolute momenta (eqn (53)).
126  p1(k+1) = 0.5d0*p1(k+1)/rk
127  p1(k+2) = 0.5d0*p1(k+2)/rk
128  p1(k+3) = 0.5d0*p1(k+3)/rk
129  22 CONTINUE
130 *
131 * Re-define variables in six dimensions (eqn (54)).
132  DO 25 k = 1,3
133  q1(k+3) = q1(k+4)
134  p1(k+3) = p1(k+4)
135  25 CONTINUE
136 *
137 * Obtain physical coordinates & momenta in c.m. frame.
138  DO 26 k = 1,3
139  q2(k+6) = -(m(1)*q1(k) + m(2)*q1(k+3))/(m(1) + m(2) + m(3))
140 * Physical coordinates of reference body M(3) (first eqn (55)).
141  q2(k) = q1(k) + q2(k+6)
142  q2(k+3) = q1(k+3) + q2(k+6)
143 * Physical coordinates of body M(1) & M(2) (second eqn (55)).
144  p2(k) = p1(k)
145  p2(k+3) = p1(k+3)
146 * Physical momenta of body M(1) & M(2) (third eqn (55)).
147  p2(k+6) = -(p2(k) + p2(k+3))
148 * Absolute momentum of reference body M(3) (last eqn (55)).
149  26 CONTINUE
150 *
151 * Specify coordinates & velocities in c.m. frame.
152  DO 30 i = 1,3
153  DO 28 k = 1,3
154  i1 = 3*i + k - 3
155  x3(k,i) = q2(i1)
156  xdot3(k,i) = p2(i1)/m(i)
157  28 CONTINUE
158  30 CONTINUE
159 *
160  IF (kdum.EQ.3) go to 40
161 *
162 * Evaluate total energy & relative energy error.
163  s1 = p2(1)**2 + p2(2)**2 + p2(3)**2
164  s2 = p2(4)**2 + p2(5)**2 + p2(6)**2
165  s3 = p2(7)**2 + p2(8)**2 + p2(9)**2
166  zke = 0.5d0*(s1/m(1) + s2/m(2) + s3/m(3))
167  s1 = m(1)*m(3)/sqrt((q2(7) - q2(1))**2 + (q2(8) - q2(2))**2 +
168  & (q2(9) - q2(3))**2)
169  s2 = m(2)*m(3)/sqrt((q2(7) - q2(4))**2 + (q2(8) - q2(5))**2 +
170  & (q2(9) - q2(6))**2)
171  s3 = m(1)*m(2)/sqrt((q2(4) - q2(1))**2 + (q2(5) - q2(2))**2 +
172  & (q2(6) - q2(3))**2)
173  ht = zke - s1 - s2 - s3
174 * Current total energy computed from physical variables.
175  cm(10) = (ht - 0.5d0*energy)/(0.5d0*energy)
176 * Relative energy error with respect to initial value.
177 *
178  40 RETURN
179 *
180  END