Nbody6
 All Files Functions Variables
rchain.f
Go to the documentation of this file.
1  SUBROUTINE rchain(IQ)
2 *
3 *
4 * Regularized integration.
5 * ------------------------
6 *
7  IMPLICIT REAL*8 (a-h,o-z)
8  REAL*8 m,mij,y(25),savex(12),savexd(12)
9  LOGICAL switch,gtype,gtype0
10  common/creg/ m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
11  & xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
12  & ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
13  common/tpr/ switch,gtype,gtype0
14  common/ind6/ ind(6)
15  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
16  common/ccoll/ qk(12),pk(12),icall,icoll,ndiss4
17 * EQUIVALENCE (Y(1),P(1))
18  SAVE
19 *
20 *
21  DO 3 i = 1,12
22  y(i) = p(i)
23  y(i+12) = q(i)
24  3 CONTINUE
25  y(25) = time4
26 *
27  IF (time4.GT.0.0d0) go to 5
28  DO 1 j = 1,12
29  savex(j) = x(j)
30  savexd(j) = xd(j)
31  1 CONTINUE
32  tsave = time4
33  tstep0 = tstep
34  istart = 0
35  next = 0
36  jc = 0
37  s = 0.0d0
38  2 zmi1 = 100.0
39  zmi2 = 100.0
40  zmi3 = 100.0
41 *
42  gtype = .false.
43  istart = istart + 1
44  runav = 0.0001*rmax4
45  epsi = eps
46  rstar = sqrt(epsr2)
47 *
48  5 DO 50 isteps = 1,10000
49  gtype0 = gtype
50  dso = ds
51  time0 = time4
52 *
53 * Advance the solution one step and order the distances.
54  CALL difsy4(25,epsi,ds,s,y)
55  DO 10 i = 1,12
56  p(i) = y(i)
57  q(i) = y(i+12)
58  10 CONTINUE
59  time4 = y(25)
60  nstep4 = nstep4 + 1
61  CALL rsort(r,switch,ind)
62 *
63 * Check time transformation type.
64  triple = (r(1) + r(2))*(r(2) + r(3))
65  IF (triple.LT.epsr2) THEN
66  gtype = .true.
67  ELSE
68  gtype = .false.
69  END IF
70 *
71 * Update smallest moment of inertia during forward integration.
72  IF (time4.GT.time0.AND.jc.EQ.0) THEN
73  imin = ind(1)
74  zmi = r(imin)**2
75  zmi1 = zmi2
76  zmi2 = zmi3
77  zmi3 = zmi
78  runav = 0.9*runav + 0.1*r(imin)
79 * Check termination during last few steps (until R(IM) > SEMI).
80  IF (iq.LT.0) go to 100
81  END IF
82 *
83 * Switch on search indicator during iteration or just after pericentre.
84  IF (icoll.LT.0) icall = 1
85  IF (r(imin).LT.rstar.AND.nstep4.GT.next) THEN
86  IF (zmi3.GT.zmi2.AND.zmi2.LT.zmi1) THEN
87  icall = 1
88  END IF
89  END IF
90 *
91 * Delay next search a few steps to avoid the same pericentre.
92  IF (icoll.GT.0) THEN
93  next = nstep4 + 2
94  go to 100
95  END IF
96 *
97 * Check the termination criterion.
98  im2 = ind(5)
99  IF (r(im2).GT.rmax4.OR.time4.GT.tmax) THEN
100  IF (r(imin).LT.0.01*runav) go to 50
101  iq = -1
102  go to 100
103  END IF
104 *
105 * Save old time transformation if change of type or new chain.
106  IF ((gtype.NEQV.gtype0).OR.switch) THEN
107  tpr0 = r(1)*r(2)*r(3)
108  IF (gtype0) tpr0 = tpr0/(r(1)*r(2) + (r(1) + r(2))*r(3))
109  END IF
110 *
111 * Check the switching condition.
112  IF (switch) THEN
113  CALL endreg
114  CALL newreg
115  nreg = nreg + 1
116  END IF
117 *
118 * Modify new step after change of type or new chain.
119  IF ((gtype.NEQV.gtype0).OR.switch) THEN
120  tpr = r(1)*r(2)*r(3)
121  IF (gtype) tpr = tpr/(r(1)*r(2) + (r(1) + r(2))*r(3))
122 * Scale regularized step by ratio of old & new time derivative.
123  ds = ds*tpr0/tpr
124  END IF
125 *
126 * Check rare case of zero regularized step.
127  IF (isteps.LT.3.AND.ds.EQ.0.0d0) ds = 1.0e-2*dso
128  IF (ds.EQ.0.0d0) go to 120
129  50 CONTINUE
130 *
131  100 RETURN
132 *
133 * Make one restart with reduced tolerance.
134  120 IF (istart.GT.1) RETURN
135  WRITE (6,125) nstep4, time4, r(ind(6)), rmax4
136  125 FORMAT (5x,' RCHAIN RESTART: # T R6 RMAX4 ',i5,3f12.6)
137  time4 = tsave
138  DO 130 j = 1,12
139  x(j) = savex(j)
140  xd(j) = savexd(j)
141  130 CONTINUE
142  switch = .false.
143  CALL newreg
144  epsi = 0.1*eps
145  ds = 0.1*tstep0/(r(1)*r(2)*r(3))
146  WRITE (6,135) icall,jc,icoll,ds
147  135 FORMAT (' RESTART: ICALL JC ICOLL FACM DS ',3i4,1p,e9.1)
148  icall = 0
149  icoll = 0
150  jc = 0
151  go to 2
152 *
153  END