Nbody6
 All Files Functions Variables
qpmod.f
Go to the documentation of this file.
1  SUBROUTINE qpmod(IM,ITERM)
2 *
3 *
4 * Modification of chain variables for tidal dissipation.
5 * ------------------------------------------------------
6 *
7  include 'commonc.h'
8  include 'common2.h'
9  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
10  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
11  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
12  common/ksave/ k1,k2
13  REAL*8 g0(3),radius(2),de(2)
14  INTEGER is(2)
15  DATA eccm,eccm2 /0.002,0.00000399/
16 *
17 *
18 * Skip for small dissipation (ISYNC > 0 delays further calls).
19  IF (kz27.EQ.1) THEN
20  r1 = max(SIZE(k1),SIZE(k2))
21  IF (abs(qperi - 4.0*r1).LT.0.01*qperi) THEN
22  isync = 1
23  iterm = 0
24  go to 90
25  END IF
26  END IF
27 *
28 * Copy radius & stellar type for interacting bodies #K1 & K2.
29  radius(1) = SIZE(k1)
30  radius(2) = SIZE(k2)
31  is(1) = istar(k1)
32  is(2) = istar(k2)
33  kstari = 0
34 *
35 * Evaluate binding energy & semi-major axis from non-singular terms.
36  CALL erel(im,eb,semi)
37 *
38 * Define mass & reduced mass for the dominant bodies.
39  mb = m(k1) + m(k2)
40  mu = m(k1)*m(k2)/mb
41 *
42 * Set energy per unit mass, eccentricity & pericentre (assume R' = 0).
43  h = eb/mu
44  ecc = 1.0 - 1.0/(rinv(im)*semi)
45  peri = 1.0/rinv(im)
46 * PERI = SEMI*(1.0D0 - ECC)
47  am0 = semi*(1.0d0 - ecc**2)
48  ndiss1 = ndiss1 + 1
49 *
50 * Choose between chaos treatment and PT or GR formulation.
51  IF (kz27.EQ.2) THEN
52  idis = 0
53  CALL chaos2(k1,k2,ecc,h,is,mb,mu,radius,semi1,ecc1,dh,idis,
54  & kstari)
55 * Exit on reaching circular orbit or spiralling stage (also collision).
56  IF (kstari.EQ.-1) go to 2
57  go to 90
58  END IF
59 *
60 * Consider sequential circularization or GR case.
61  IF (kz27.EQ.1) THEN
62  am0 = semi*(1.0d0 - ecc**2)
63  ecc2 = eccm2
64  ecc1 = sqrt(ecc2)
65  acirc = am0/(1.0 - ecc2)
66 * Accept circularized orbit directly if ACIRC < 4*R.
67  IF (acirc.LT.4.0*r1) THEN
68  semi1 = acirc
69  ELSE
70 * Obtain E1 by (1 + E1) = AM0/(4*R1) and A1 by A1*(1 - E1) = 4*R1.
71  ecc1 = 0.25*am0/r1 - 1.0
72  ecc1 = max(ecc1,eccm)
73 * Set semi-major axis from angular momentum conservation.
74  semi1 = am0/(1.0 - ecc1**2)
75  END IF
76 * Form the corresponding energy change.
77  dh = 0.5*mb*(1.0/semi - 1.0/semi1)
78  de(1) = -mu*dh
79  de(2) = 0.0
80  ELSE
81 * Obtain the tidal energy change for GR.
82  CALL tides3(qperi,m(k1),m(k2),vstar1,h,ecc,de)
83 * Include safety check on energy loss to prevent new SEMI < R.
84  dh = -(de(1) + de(2))/mu
85  IF (h + dh.LT.-0.5*mb*rinv(im)) THEN
86  dh = -0.5*mb*rinv(im) - h
87  END IF
88  semi1 = -0.5*mb/(h + dh)
89  ecc1 = 1.0 - peri/semi1
90  ecc1 = max(ecc1,0.0d0)
91 * Note: Do not impose minimum ECCM unless with 2nd alternative C2.
92 * ECC1 = MAX(ECC1,ECCM)
93  END IF
94 *
95 * Determine new eccentricity from angular momentum conservation.
96 * ECC2 = ECC**2 + 2.0D0*AM0*DH/MB
97 * ECC2 = MAX(ECC2,ECCM2)
98 * ECC1 = SQRT(ECC2)
99 *
100 * Adopt instantaneous circularization instead of standard PT.
101 * ECC2 = ECCM2
102 * ECC1 = SQRT(ECC2)
103 * SEMI1 = AM0/(1.0 - ECC2)
104 * DH = 0.5*MB*(1.0/SEMI - 1.0/SEMI1)
105 * DE(1) = -MU*DH
106 * DE(2) = 0.0
107 
108 * Skip on final hyperbolic energy.
109  IF (h + dh.GT.0.0) go to 90
110 *
111 * Set new pericentre and update binding energy.
112  2 peri1 = semi1*(1.0d0 - ecc1)
113  hi = h
114  h = h + dh
115 *
116 * Form KS coordinate scaling factor from pericentre ratio.
117  c1 = sqrt(peri1/peri)
118 *
119 * Specify KS velocity scaling for conserved angular momentum.
120  IF (kz27.EQ.1.OR.kstari.EQ.-2) THEN
121  c2 = 1.0/c1**2
122  ELSE
123  c2 = sqrt((mb + h*peri1)/(mb + hi*peri))
124 * Note: since PERI1 = PERI this is same as used for KZ27 = 2.
125  END IF
126 *
127 * See whether circular orbit condition applies.
128  IF (ecc1.LE.eccm.AND.kz27.EQ.1) THEN
129  am = semi1*(1.0d0 - ecc1**2)
130  c2 = (am/am0)/c1**2
131  END IF
132 *
133 * Include alternative formulation of Mikkola (which is equivalent).
134  IF (kz27.LT.-1) THEN
135 * Note this derivation may contain assumption of J = const (10/6/99).
136  v02 = mb*(2.0/qperi - 1.0/semi)
137  a1 = mb/(qperi*v02)
138  a2 = mb/(v02*semi1)
139 * Adopt reduced energy loss in case of imaginary solution.
140  IF (a1**2.GT.a2) THEN
141  a3 = a1 + sqrt(a1**2 - a2)
142  ELSE
143  a3 = a1
144  WRITE (6,5) a1, a2, sqrt(1.0/a3), a3, c1
145  5 FORMAT (' WARNING! QPMOD A1 A2 SQRT(1/A3) A3 C1 ',
146  & 1p,2e10.2,0p,3f12.6)
147  END IF
148  c2 = a3
149  END IF
150 *
151 * Derive scale factor from velocity ratio for chaotic motion.
152  IF (kz27.EQ.2.AND.kstari.EQ.-1) THEN
153  c2 = sqrt((h + mb/peri1)/(hi + mb/peri))
154 * Note the use of physical velocity here vs UDOT in routine KSTIDE.
155  END IF
156 *
157 * Modify the KS coordinates for dominant bodies and update distance.
158  ks = 4*(im - 1)
159  rm = 0.0d0
160  DO 10 k = 1,4
161  q(ks+k) = c1*q(ks+k)
162  rm = rm + q(ks+k)**2
163  10 CONTINUE
164  rinv(im) = 1.0/rm
165 *
166 * Change the corresponding physical momenta (IM & IM + 1 in chain).
167  j1 = 3*(im - 1)
168  j2 = 3*im
169  DO 20 k = 1,3
170 * P11 = A3*(M(K2)*PI(J1+K) - M(K1)*PI(J2+K))/MB
171  vb = (pi(j1+k) + pi(j2+k))/mb
172 * P1K = -MU*(1.0 - C2)*(PI(J1+K)/M(K1) - PI(J2+K)/M(K2))
173  p1k = c2*(m(k2)*pi(j1+k) - m(k1)*pi(j2+k))/mb
174  pi(j1+k) = m(k1)*vb + p1k
175  pi(j2+k) = m(k2)*vb - p1k
176 * PI(J1+K) = PI(J1+K) + P1K
177 * PI(J2+K) = PI(J2+K) - P1K
178  20 CONTINUE
179 *
180 * Form physical chain momenta (first & last, then intermediate values).
181  l = 3*(n - 2)
182  DO 30 k = 1,3
183  wc(k) = -pi(k)
184  wc(l+k) = pi(l+k+3)
185  30 CONTINUE
186 *
187  DO 40 i = 2,n-2
188  l = 3*(i - 1)
189  DO 35 k = 1,3
190  wc(l+k) = wc(l+k-3) - pi(l+k)
191  35 CONTINUE
192  40 CONTINUE
193 *
194 * Re-determine all regularized momenta by KS transformation.
195  DO 50 l = 1,n-1
196  l1 = 3*(l - 1) + 1
197  l2 = l1 + 1
198  l3 = l2 + 1
199  ks1 = 4*(l - 1) + 1
200  ks2 = ks1 + 1
201  ks3 = ks2 + 1
202  ks4 = ks3 + 1
203  p(ks1) = 2.d0*(+q(ks1)*wc(l1) + q(ks2)*wc(l2) + q(ks3)*wc(l3))
204  p(ks2) = 2.d0*(-q(ks2)*wc(l1) + q(ks1)*wc(l2) + q(ks4)*wc(l3))
205  p(ks3) = 2.d0*(-q(ks3)*wc(l1) - q(ks4)*wc(l2) + q(ks1)*wc(l3))
206  p(ks4) = 2.d0*(+q(ks4)*wc(l1) - q(ks3)*wc(l2) + q(ks2)*wc(l3))
207  50 CONTINUE
208 *
209 * Evaluate potential energy due to non-chained distances (KS terms OK).
210  it = 0
211  52 pot1 = 0.0
212  ij = n - 1
213  DO 60 i = 1,n-2
214  DO 55 j = i+2,n
215  ij = ij + 1
216  pot1 = pot1 + mij(i,j)*rinv(ij)
217  55 CONTINUE
218  60 CONTINUE
219 *
220 * Obtain the potential energy of the new system (KS terms unchanged).
221  IF (it.EQ.0) THEN
222 * Save modified configuration in QK & PK (for EREL & TRANSK).
223  DO 70 i = 1,n-1
224  ks = 4*(i - 1)
225  DO 65 j = 1,4
226  qk(ks+j) = q(ks+j)
227  pk(ks+j) = p(ks+j)
228  65 CONTINUE
229  70 CONTINUE
230 *
231 * Set new values of inverse distances (only non-chained values needed).
232  CALL transk
233  pot0 = pot1
234  it = it + 1
235  go to 52
236  END IF
237 *
238 * Update collision energy & internal energy (MU*DH is not sufficient).
239  ecoll1 = ecoll1 - mu*dh + (pot1 - pot0)
240  energy = energy + mu*dh - (pot1 - pot0)
241 *
242  IF (kz30.GT.1) THEN
243 * Obtain consistent value of the total energy (for checking purposes).
244  CALL transx
245  CALL const(x,v,m,n,en1,g0,alag)
246  dpot = pot1 - pot0
247  WRITE (6,75) en1, en1-energy, mu*dh, semi, semi1, dpot
248  75 FORMAT (/,' QPMOD: EN1 ERR DE A A1 DP ',
249  & f10.6,1p,5e10.2)
250  END IF
251 *
252 * Activate indicator for synchronous orbit.
253  IF (ecc.GT.eccm.AND.ecc1.LE.eccm) THEN
254  isync = 1
255  END IF
256 *
257 * Print diagnostic if eccentricity > 0.99 or #30 > 1.
258  IF (ecc.GT.0.99.OR.kz30.GT.1) THEN
259  WRITE (6,80) namec(k1), namec(k2), semi1, ecc, ecc1, h, qperi
260  80 FORMAT (' NEW QPMOD NAM AF E0 EF H QP ',
261  & 2i6,1p,e10.2,0p,2f8.4,f9.1,1p,2e10.2)
262  END IF
263 *
264 * Perform stability test (ITERM = -1: termination; = -2: reduction).
265  CALL stablc(im,iterm,semi1)
266 *
267  90 RETURN
268 *
269  END