https://github.com/nbodyx/Nbody6
Tip revision: dff3b5c68673d4d4c4c9f54c8ba110b8416098f2 authored by nitadori on 31 May 2020, 13:04:00 UTC
Avoided touching uninitialized SEMIX
Avoided touching uninitialized SEMIX
Tip revision: dff3b5c
rchain.f
SUBROUTINE RCHAIN(IQ)
*
*
* Regularized integration.
* ------------------------
*
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 M,MIJ,Y(25),SAVEX(12),SAVEXD(12)
LOGICAL SWITCH,GTYPE,GTYPE0
COMMON/CREG/ M(4),X(12),XD(12),P(12),Q(12),TIME4,ENERGY,EPSR2,
& XR(9),W(9),R(6),TA(6),MIJ(6),CM(10),RMAX4,TMAX,
& DS,TSTEP,EPS,NSTEP4,NAME4(4),KZ15,KZ27,NREG,NFN
COMMON/TPR/ SWITCH,GTYPE,GTYPE0
COMMON/IND6/ IND(6)
COMMON/CLOSE/ RIJ(4,4),RCOLL,QPERI,SIZE(4),ECOLL3,IP(4)
COMMON/CCOLL/ QK(12),PK(12),ICALL,ICOLL,NDISS4
* EQUIVALENCE (Y(1),P(1))
SAVE
*
*
DO 3 I = 1,12
Y(I) = P(I)
Y(I+12) = Q(I)
3 CONTINUE
Y(25) = TIME4
*
IF (TIME4.GT.0.0D0) GO TO 5
DO 1 J = 1,12
SAVEX(J) = X(J)
SAVEXD(J) = XD(J)
1 CONTINUE
TSAVE = TIME4
TSTEP0 = TSTEP
ISTART = 0
NEXT = 0
JC = 0
S = 0.0D0
2 ZMI1 = 100.0
ZMI2 = 100.0
ZMI3 = 100.0
*
GTYPE = .FALSE.
ISTART = ISTART + 1
RUNAV = 0.0001*RMAX4
EPSI = EPS
RSTAR = SQRT(EPSR2)
*
5 DO 50 ISTEPS = 1,10000
GTYPE0 = GTYPE
DSO = DS
TIME0 = TIME4
*
* Advance the solution one step and order the distances.
CALL DIFSY4(25,EPSI,DS,S,Y)
DO 10 I = 1,12
P(I) = Y(I)
Q(I) = Y(I+12)
10 CONTINUE
TIME4 = Y(25)
NSTEP4 = NSTEP4 + 1
CALL RSORT(R,SWITCH,IND)
*
* Check time transformation type.
TRIPLE = (R(1) + R(2))*(R(2) + R(3))
IF (TRIPLE.LT.EPSR2) THEN
GTYPE = .TRUE.
ELSE
GTYPE = .FALSE.
END IF
*
* Update smallest moment of inertia during forward integration.
IF (TIME4.GT.TIME0.AND.JC.EQ.0) THEN
IMIN = IND(1)
ZMI = R(IMIN)**2
ZMI1 = ZMI2
ZMI2 = ZMI3
ZMI3 = ZMI
RUNAV = 0.9*RUNAV + 0.1*R(IMIN)
* Check termination during last few steps (until R(IM) > SEMI).
IF (IQ.LT.0) GO TO 100
END IF
*
* Switch on search indicator during iteration or just after pericentre.
IF (ICOLL.LT.0) ICALL = 1
IF (R(IMIN).LT.RSTAR.AND.NSTEP4.GT.NEXT) THEN
IF (ZMI3.GT.ZMI2.AND.ZMI2.LT.ZMI1) THEN
ICALL = 1
END IF
END IF
*
* Delay next search a few steps to avoid the same pericentre.
IF (ICOLL.GT.0) THEN
NEXT = NSTEP4 + 2
GO TO 100
END IF
*
* Check the termination criterion.
IM2 = IND(5)
IF (R(IM2).GT.RMAX4.OR.TIME4.GT.TMAX) THEN
IF (R(IMIN).LT.0.01*RUNAV) GO TO 50
IQ = -1
GO TO 100
END IF
*
* Save old time transformation if change of type or new chain.
IF ((GTYPE.NEQV.GTYPE0).OR.SWITCH) THEN
TPR0 = R(1)*R(2)*R(3)
IF (GTYPE0) TPR0 = TPR0/(R(1)*R(2) + (R(1) + R(2))*R(3))
END IF
*
* Check the switching condition.
IF (SWITCH) THEN
CALL ENDREG
CALL NEWREG
NREG = NREG + 1
END IF
*
* Modify new step after change of type or new chain.
IF ((GTYPE.NEQV.GTYPE0).OR.SWITCH) THEN
TPR = R(1)*R(2)*R(3)
IF (GTYPE) TPR = TPR/(R(1)*R(2) + (R(1) + R(2))*R(3))
* Scale regularized step by ratio of old & new time derivative.
DS = DS*TPR0/TPR
END IF
*
* Check rare case of zero regularized step.
IF (ISTEPS.LT.3.AND.DS.EQ.0.0D0) DS = 1.0E-2*DSO
IF (DS.EQ.0.0D0) GO TO 120
50 CONTINUE
*
100 RETURN
*
* Make one restart with reduced tolerance.
120 IF (ISTART.GT.1) RETURN
WRITE (6,125) NSTEP4, TIME4, R(IND(6)), RMAX4
125 FORMAT (5X,' RCHAIN RESTART: # T R6 RMAX4 ',I5,3F12.6)
TIME4 = TSAVE
DO 130 J = 1,12
X(J) = SAVEX(J)
XD(J) = SAVEXD(J)
130 CONTINUE
SWITCH = .FALSE.
CALL NEWREG
EPSI = 0.1*EPS
DS = 0.1*TSTEP0/(R(1)*R(2)*R(3))
WRITE (6,135) ICALL,JC,ICOLL,DS
135 FORMAT (' RESTART: ICALL JC ICOLL FACM DS ',3I4,1P,E9.1)
ICALL = 0
ICOLL = 0
JC = 0
GO TO 2
*
END