Nbody6
 All Files Functions Variables
qpmod4.f
Go to the documentation of this file.
1  SUBROUTINE qpmod4(IM,ITERM)
2 *
3 *
4 * Modification of chain variables for tidal dissipation.
5 * ------------------------------------------------------
6 *
7  IMPLICIT REAL*8 (a-h,m,o-z)
8  common/creg/ m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
9  & xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
10  & ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
11  common/iconf/ i1,i2,i3,i4
12  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
13  common/ccoll/ qk(12),pk(12),icall,icoll,ndiss4
14  common/savep/ pi(12)
15  common/ksave/ k1,k2
16  REAL*8 de(2),radius(2)
17  INTEGER is(2)
18 *
19 *
20 * Obtain tidal energy loss (bodies #K1 & K2 with separation QPERI).
21  radius(1) = SIZE(k1)
22  radius(2) = SIZE(k2)
23  is(1) = ip(k1)
24  is(2) = ip(k2)
25 *
26  CALL tides(qperi,m(k1),m(k2),radius(1),radius(2),is,de)
27 *
28 * Evaluate binding energy & semi-major axis from non-singular terms.
29  CALL erel4(im,eb,semi)
30 *
31 * Define mass & reduced mass for the dominant bodies.
32  mb = m(k1) + m(k2)
33  mu = m(k1)*m(k2)/mb
34 *
35 * Set energy per unit mass, eccentricity & pericentre (assume R' = 0).
36  h = eb/mu
37  ecc = 1.0 - r(im)/semi
38  peri = semi*(1.0d0 - ecc)
39 *
40 * Determine new eccentricity from angular momentum conservation.
41  dh = -(de(1) + de(2))/mu
42  am0 = semi*(1.0d0 - ecc**2)
43  ecc2 = ecc**2 + 2.0d0*am0*dh/mb
44  IF (ecc2.GT.0.0d0) THEN
45  ecc1 = sqrt(ecc2)
46  ELSE
47  ecc1 = 0.001
48  END IF
49 *
50 * Update binding energy and set new semi-major axis & pericentre.
51  h1 = h + dh
52  semi1 = -0.5d0*mb/h1
53  peri1 = semi1*(1.0d0 - ecc1)
54 *
55 * Form KS coordinate scaling factor from pericentre ratio.
56  c1 = sqrt(peri1/peri)
57 *
58 * Specify velocity scaling from angular momentum conservation.
59  c2 = 1.0/c1
60 *
61 * See whether circular orbit condition applies.
62  IF (ecc1.LE.0.001) THEN
63  am = semi1*(1.0d0 - ecc1**2)
64  c2 = sqrt(am/am0)/c1
65  END IF
66 *
67 * Modify the KS coordinates for dominant bodies and update distance.
68  ks = 4*(im - 1)
69  r(im) = 0.0d0
70  DO 20 k = 1,4
71  q(ks+k) = c1*q(ks+k)
72  r(im) = r(im) + q(ks+k)**2
73  20 CONTINUE
74 *
75 * Set modified physical momenta for the two critical bodies (K1 & K2).
76  j1 = 3*(k1 - 1)
77  j2 = 3*(k2 - 1)
78 * Note PI is associated with the bodies (hence J1 # 3*(IM - 1).
79  DO 30 k = 1,3
80  p1k = -mu*(1.0 - c2**2)*(pi(j1+k)/m(k1) - pi(j2+k)/m(k2))
81  pi(j1+k) = pi(j1+k) + p1k
82  pi(j2+k) = pi(j2+k) - p1k
83  30 CONTINUE
84 *
85 * Form physical chain momenta.
86  ip1 = 3*(i1 - 1)
87  ip2 = 3*(i2 - 1)
88  ip4 = 3*(i4 - 1)
89  DO 40 k = 1,3
90  w(k ) = -pi(ip1+k)
91  w(k+3) = -pi(ip1+k) - pi(ip2+k)
92  w(k+6) = +pi(ip4+k)
93  40 CONTINUE
94 *
95 * Re-determine all regularized momenta by KS transformation.
96  DO 50 l = 1,3
97  l1 = 3*(l - 1) + 1
98  l2 = l1 + 1
99  l3 = l2 + 1
100  lq1 = 4*(l - 1) + 1
101  lq2 = lq1 + 1
102  lq3 = lq2 + 1
103  lq4 = lq3 + 1
104  p(lq1) = 2.d0*(+q(lq1)*w(l1) + q(lq2)*w(l2) + q(lq3)*w(l3))
105  p(lq2) = 2.d0*(-q(lq2)*w(l1) + q(lq1)*w(l2) + q(lq4)*w(l3))
106  p(lq3) = 2.d0*(-q(lq3)*w(l1) - q(lq4)*w(l2) + q(lq1)*w(l3))
107  p(lq4) = 2.d0*(+q(lq4)*w(l1) - q(lq3)*w(l2) + q(lq2)*w(l3))
108  50 CONTINUE
109 *
110 * Obtain consistent value of the total energy (instead of correcting).
111  CALL endreg
112  e0 = energy
113  CALL newsys(x,xd,m,4,energy,gam)
114 *
115 * Update diagnostic variables (Note: DE(1) + DE(2) is not sufficient).
116 * ECOLL3 = ECOLL3 + (DE(1) + DE(2))
117  ecoll3 = ecoll3 + (e0 - energy)
118  ndiss4 = ndiss4 + 1
119 *
120 * Perform stability test (ITERM < 0 denotes termination).
121  CALL stabl4(iterm)
122 *
123 * Print diagnostic if eccentricity > 0.99.
124  IF (ecc.GT.0.99) THEN
125  WRITE (6,60) name4(k1), name4(k2), semi1, ecc, ecc1, h, qperi
126  60 FORMAT (' NEW QPMOD4 NAM AF E0 EF H QP ',
127  & 2i5,1p,e10.2,0p,2f8.3,f9.1,1p,e10.2)
128  END IF
129 *
130  RETURN
131 *
132  END