Nbody6
qpmod3.f
Go to the documentation of this file.
1  SUBROUTINE qpmod3(IM,ITERM)
2 *
3 *
4 * Modification of AZ variables for tidal dissipation.
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  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
12  common/azcoll/ rk(3),qk(8),pk(8),icall,icoll,ndiss3
14  INTEGER is(2)
15 *
16 *
17 * Obtain tidal energy loss (bodies #IM & 3 with separation QPERI).
20  is(1) = ip(im)
21  is(2) = ip(3)
22 *
24 *
25 * Determine physical momenta of M(IM) & M(3-IM).
26  DO 10 kcomp = 1,2
27  k = 4*(kcomp - 1)
28 *
29 * Form product of half Levi-Civita matrix & regularized momentum.
30  p1(k+1) = q(k+1)*p(k+1) - q(k+2)*p(k+2) - q(k+3)*p(k+3) +
31  & q(k+4)*p(k+4)
32  p1(k+2) = q(k+2)*p(k+1) + q(k+1)*p(k+2) - q(k+4)*p(k+3) -
33  & q(k+3)*p(k+4)
34  p1(k+3) = q(k+3)*p(k+1) + q(k+4)*p(k+2) + q(k+1)*p(k+3) +
35  & q(k+2)*p(k+4)
36  p1(k+4) = 0.0d0
37 *
38 * Divide by the distance to obtain physical momentum (eqn (53)).
39  ri = q(k+1)**2 + q(k+2)**2 + q(k+3)**2 + q(k+4)**2
40  p1(k+1) = 0.5d0*p1(k+1)/ri
41  p1(k+2) = 0.5d0*p1(k+2)/ri
42  p1(k+3) = 0.5d0*p1(k+3)/ri
43  10 CONTINUE
44 *
45 * Define KS index for closest bodies and second interaction.
46  k1 = 4*(im - 1)
47  k2 = 4*(2 - im)
48 *
49 * Evaluate kinetic energy term of least dominant interaction.
50  p2 = 0.0d0
51  DO 15 k = 1,3
52  p2 = p2 + p1(k2+k)**2
53 * Form relative momentum of dominant motion from c.m. condition.
54  p1(k1+k) = 2.0d0*p1(k1+k) + p1(k2+k)
55  15 CONTINUE
56 *
57 * Set consistent third distance for energy relation.
58  a21 = q(1)*q(1) - q(2)*q(2) - q(3)*q(3) + q(4)*q(4)
59  & - q(5)*q(5) + q(6)*q(6) + q(7)*q(7) - q(8)*q(8)
60  a22 = q(1)*q(2) - q(3)*q(4) - q(5)*q(6) + q(7)*q(8)
61  a23 = q(1)*q(3) + q(2)*q(4) - q(5)*q(7) - q(6)*q(8)
62  a22 = a22 + a22
63  a23 = a23 + a23
64  r3 = sqrt(a21*a21 + a22*a22 + a23*a23)
65 *
66 * Define radius, mass & reduced mass for the dominant bodies.
67  rm = min(r1,r2)
68  mb = m(im) + m(3)
69  mu = m(im)*m(3)/mb
70 *
71 * Obtain binding energy from total energy and perturbing function.
72  mu2 = m(3-im)*mb/(m(3-im) + mb)
73  vp = 0.5d0*p2/mu2 - m(1)*m(2)/r3 - m(3-im)*m(3)/max(r1,r2)
74  h = (0.5d0*energy - vp)/mu
75 *
76 * Set semi-major axis, eccentricity & pericentre (assume R' = 0).
77  semi = -0.5d0*mb/h
78  ecc = 1.0 - rm/semi
79  peri = semi*(1.0d0 - ecc)
80 *
81 * Determine new eccentricity from angular momentum conservation.
82  dh = -(de(1) + de(2))/mu
83  am0 = semi*(1.0d0 - ecc**2)
84  ecc2 = ecc**2 + 2.0d0*am0*dh/mb
85  IF (ecc2.GT.0.0d0) THEN
86  ecc1 = sqrt(ecc2)
87  ELSE
88  ecc1 = 0.0
89  END IF
90 *
91 * Update binding energy and set new semi-major axis & pericentre.
92  h1 = h + dh
93  semi1 = -0.5d0*mb/h1
94  peri1 = semi1*(1.0d0 - ecc1)
95 *
96 * Form KS coordinate scaling factor from pericentre ratio.
97  c1 = sqrt(peri1/peri)
98 * Specify KS velocity scaling from angular momentum conservation.
99  c2 = 1.0/c1
100 *
101 * Modify the KS coordinates & physical momentum of dominant motion.
102  ri = 0.0d0
103  DO 20 k = 1,4
104  q(k1+k) = c1*q(k1+k)
105  ri = ri + q(k1+k)**2
106  p1(k1+k) = c2**2*p1(k1+k)
107  p1(k1+k) = 0.5d0*(p1(k1+k) - p1(k2+k))
108  20 CONTINUE
109 *
110 * Transform to new regularized momentum (eqn (50)).
111  k = k1
112  p(k+1) = 2.0d0*(+q(k+1)*p1(k+1) + q(k+2)*p1(k+2) +
113  & q(k+3)*p1(k+3))
114  p(k+2) = 2.0d0*(-q(k+2)*p1(k+1) + q(k+1)*p1(k+2) +
115  & q(k+4)*p1(k+3))
116  p(k+3) = 2.0d0*(-q(k+3)*p1(k+1) - q(k+4)*p1(k+2) +
117  & q(k+1)*p1(k+3))
118  p(k+4) = 2.0d0*(+q(k+4)*p1(k+1) - q(k+3)*p1(k+2) +
119  & q(k+2)*p1(k+3))
120 *
121 * Update the smallest distance (new pericentre has increased).
122  IF (im.EQ.1) THEN
123  r1 = ri
124  ELSE
125  r2 = ri
126  END IF
127 *
128 * Transform to physical variables and obtain consistent energy.
129  CALL trans3(3)
130  CALL trans3(0)
131 *
132 * Update diagnostic variables (note DE > 0).
133  ecoll3 = ecoll3 + (de(1) + de(2))
134  ndiss3 = ndiss3 + 1
135 *
136 * Perform stability test (ITERM < 0 denotes termination).
137  CALL stabl3(iterm)
138 *
139 * Print diagnostic if eccentricity > 0.99.
140  IF (ecc.GT.0.99) THEN
141  WRITE (6,30) name3(im), name3(3), semi1, ecc, ecc1, h, qperi
142  30 FORMAT (' NEW QPMOD3 NAM AF E0 EF H QP ',
143  & 2i5,1p,e10.2,0p,2f8.3,f9.1,1p,e10.2)
144  CALL flush(6)
145  END IF
146 *
147  RETURN
148 *
149  END