Nbody6
invert.f
Go to the documentation of this file.
1  SUBROUTINE invert(DT,DTAU)
2 *
3 *
4 * Inversion of physical time interval.
5 * ------------------------------------
6 *
7  include 'commonc.h'
8  include 'common2.h'
9  LOGICAL kslow,kcoll
10  REAL*8 vi(nmx3),vc(nmx3),ksch,c(5)
11  common/slow1/ tk2(0:nmx),ejump,ksch(nmx),kslow,kcoll
12  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
13  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
14  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
15  common/calls/ tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
16  common/clump/ bodys(nmx,5),t0s(5),ts(5),steps(5),rmaxs(5),
17  & names(nmx,5),isys(5)
18  common/ksave/ k1,k2
19  SAVE
20  DATA iwarn/0/
21 *
22 *
23 * Transform to current chain coordinates and set RINV.
24  DO i=1,n-1
25  l1=3*(i-1)+1
26  ks1=4*(i-1)+1
27  CALL ksphys(q(ks1),p(ks1),xc(l1),wc(l1))
28  rk = q(ks1)**2 + q(ks1+1)**2 + q(ks1+2)**2 + q(ks1+3)**2
29  rinv(i) = 1.0/rk
30  END DO
31 *
32 * Obtain physical velocities (first absolute, then relative).
33  l = 3*(n-2)
34  DO k = 1,3
35  vi(k) = -wc(k)/mc(1)
36  vi(l+k+3) = wc(l+k)/mc(n)
37  END DO
38  DO i = 2,n-1
39  l = 3*(i-1)
40  DO k = 1,3
41  vi(l+k) = (wc(l+k-3) - wc(l+k))/mc(i)
42  END DO
43  END DO
44  do j = 1,3*(n-1)
45  vc(j) = vi(j+3) - vi(j)
46  end do
47 *
48 * Determine the smallest two-body distance and chain index.
49  rm = 0.0
50  do i = 1,n-1
51  if (rinv(i).gt.rm) then
52  rm = rinv(i)
53  i1 = i
54  end if
55  end do
56  rm = 1.0/rm
57 *
58 * Form semi-major axes for the two closest distances.
59  dm = 0.0
60  do i = 1,n-1
61  l = 3*(i - 1)
62  if (i.eq.i1) then
63  w2 = vc(l+1)**2 + vc(l+2)**2 + vc(l+3)**2
64  mb = mc(i) + mc(i+1)
65  amax = 2.0*rinv(i) - w2/mb
66  a1 = 1.0/amax
67  v21 = w2
68 * Obtain the largest relative perturbation from nearest neighbour.
69  if (i.eq.1) then
70  g1 = 2.0*mc(3)/mb*(rinv(2)/rinv(1))**3
71  else
72  g1 = 2.0*mc(i-1)/mb*(rinv(i-1)/rinv(i))**3
73  if (i.lt.n-1) then
74  g1b = 2.0*mc(i+2)/mb*(rinv(i+1)/rinv(i))**3
75  g1 = max(g1b,g1)
76  end if
77  end if
78 * Treat the second smallest distance similarly.
79  else if (rinv(i).gt.dm) then
80  w2 = vc(l+1)**2 + vc(l+2)**2 + vc(l+3)**2
81  mb = mc(i) + mc(i+1)
82  amax = 2.0*rinv(i) - w2/mb
83  a2 = 1.0/amax
84  i2 = i
85  dm = rinv(i)
86  if (i.eq.1) then
87  g2 = 2.0*mc(3)/mb*(rinv(2)/rinv(1))**3
88  else
89  g2 = 2.0*mc(i-1)/mb*(rinv(i-1)/rinv(i))**3
90  if (i.lt.n-1) then
91  g2b = 2.0*mc(i+2)/mb*(rinv(i+1)/rinv(i))**3
92  g2 = max(g2b,g2)
93  end if
94  end if
95  end if
96  end do
97  dm = 1.0/dm
98 *
99 * Obtain reliable semi-major axis for small pericentre or large EB.
100  eb = -0.5*mc(i1)*mc(i1+1)/a1
101  eb1 = eb/energy
102  IF (rm.lt.0.2*a1.or.eb1.gt.0.9) then
103 * Copy current configuration for EREL & TRANSK.
104  DO 2 i = 1,n-1
105  ks = 4*(i - 1)
106  DO 1 j = 1,4
107  qk(ks+j) = q(ks+j)
108  pk(ks+j) = p(ks+j)
109  1 CONTINUE
110  2 CONTINUE
111 *
112 * Evaluate semi-major axis of K1 & K2 from non-singular variables.
113  k1 = iname(i1)
114  k2 = iname(i1+1)
115  CALL erel(i1,eb,semi)
116 *
117  err = (a1 - semi)/semi
118  IF (abs(err).GT.1.0e-05.AND.iwarn.LT.10) THEN
119  iwarn = iwarn + 1
120  zk = 1.0
121  DO i = 1,n-1
122  zk = max(ksch(i),zk)
123  END DO
124  WRITE (6,4) nstep1, zk, g1, rm, semi, err
125  4 FORMAT (' WARNING! INVERT # K g1 RM SEMI DA/A ',
126  & i6,f7.2,1p,4e10.2)
127  END IF
128 * Replace direct value by improved determination.
129  a1 = semi
130  END IF
131 *
132 * Restrict the interval to a fractional period or crossing time.
133  IF (a1.GT.0.0) THEN
134  tks = 6.28*ksch(i1)*a1*sqrt(a1/(mc(i1) + mc(i1+1)))
135 * Adopt upper limit of TKS/2 but allow genuine small DT.
136  dt = min(0.5*tks,dt)
137  ELSE
138 * Use crossing time for hyperbolic velocity (TPR compensates small R).
139  dt = min(rm/sqrt(v21),dt)
140  END IF
141 *
142 * Check second closest pair if bound (small period; a1 < 0 possible).
143  if (a2.gt.0.0) then
144  tks = 6.28*ksch(i2)*a2*sqrt(a2/(mc(i2) + mc(i2+1)))
145  dt = min(0.5*tks,dt)
146  end if
147 *
148 * Initialize time variables.
149  dt0 = dt
150  sum1 = 1.0/tpr
151  sum2 = 0.0
152 *
153 * Skip iteration for significant perturbations.
154  IF (g1.GT.0.01) THEN
155  dtau = sum1*dt
156  go to 30
157  END IF
158 *
159 * Begin the loop (second time depends on conditions).
160  5 dt = dt0
161  mb = mc(i1) + mc(i1+1)
162 * Form useful quantities (including EB and eccentricity if needed).
163  beta = mb/a1
164  li = 3*(i1 - 1)
165  eta = xc(li+1)*vc(li+1) + xc(li+2)*vc(li+2) + xc(li+3)*vc(li+3)
166  zeta = mb - beta*rm
167 * EB = -0.5*MC(I1)*MC(I1+1)/a1
168 * ECC = SQRT((1.0 - RM/a1)**2 + ETA**2/(MB*a1))
169 *
170 * Sum the constant part of the integral.
171  sum1 = sum1 - 2.0*mc(i1)*mc(i1+1)/(ksch(i1)*rm)
172 * Reduce physical time by slow-down factor for the two-body integral.
173  dt = dt/ksch(i1)
174 *
175 * Obtain improved initial guess by factor of 2 bisecting (Chris Tout).
176  it = 0
177 * Choose a conservative starting value in case of small RM.
178  y = dt*min(1.0/rm,0.5/abs(a1))
179  8 z = beta*y**2
180  CALL cfuncs(z,c)
181  y0 = dt - ((zeta*c(3)*y + eta*c(2))*y + rm)*y
182  it = it + 1
183  IF (y0.GT.0.0) THEN
184  y = 2.0*y
185  ELSE
186  y = 0.5*y
187  END IF
188  IF (it.LT.4) go to 8
189 *
190 * Perform Newton-Raphson iteration using first and second derivatives.
191  10 z = beta*y**2
192 *
193 * Obtain new c-functions for argument Z = -2h*Y^2.
194  CALL cfuncs(z,c)
195 *
196 * Define the function Y0 and its derivatives (cf. Book eqns. 12.9).
197  y0 = dt - ((zeta*c(3)*y + eta*c(2))*y + rm)*y
198  ypr = -((zeta*c(2)*y + eta*c(1))*y + rm)
199  y2pr = -eta*(1.0 - z*c(2)) - zeta*y*c(1)
200 * Adopt safety measure to avoid negative argument (Seppo Mikkola).
201  d = y0*y2pr/ypr**2
202  y = y - y0/ypr/sqrt(0.01 + abs(0.99d0 - d))
203  dt1 = ((zeta*c(3)*y + eta*c(2))*y + rm)*y
204  it = it + 1
205  IF (abs(dt - dt1).GT.1.0d-06*abs(dt).AND.it.LT.25) go to 10
206 *
207 * Accumulate the integral of dt/R.
208  sum2 = sum2 + 2.0*mc(i1)*mc(i1+1)*y
209 *
210  IF (it.GT.15.AND.iwarn.LT.50) THEN
211  iwarn = iwarn + 1
212  ecc = sqrt((1.0 - rm/a1)**2 + eta**2/(mb*a1))
213  WRITE (6,15) it, ksch(i1), ecc, rm, a1, dt-dt1, dtau
214  15 FORMAT (' WARNING! INVERT IT K E R A ERR DTAU ',
215  & i4,f7.2,f7.3,1p,5e9.1)
216  END IF
217 *
218 * Consider the second close interaction for small perturbation.
219  IF (i2.ne.i1.and.g2.lt.0.01) THEN
220  i1 = i2
221  a1 = a2
222  rm = dm
223  go to 5
224  END IF
225 *
226 * Combine the two terms and ensure positive value.
227  dtau = sum1*dt0 + sum2
228  IF (dtau.LT.0.0) THEN
229  dtau = abs(sum1)*dt0
230  IF (iwarn.LT.100) THEN
231  iwarn = iwarn + 1
232  WRITE (6,16) ksch(i1), rm/a1, sum1*dt0, sum2, dtau
233  16 FORMAT (' INVERT NEGATIVE! KSCH R/A S1*DT0 S2 DTAU ',
234  & f7.2,1p,4e10.2)
235  END IF
236  END IF
237 *
238 * WRITE (6,20) IT, a1, DT, DT-DT1, DTAU, STEP
239 * 20 FORMAT (' INVERT: IT A DT ERR DTAU STEP ',I4,1P,5E10.2)
240 *
241  30 RETURN
242 *
243  END