UNOFFICIAL COPY FROM 2006 07 18 -------------------------------- Refer to Sverre Aarseth's download site at http://www.ast.cam.ac.uk/~sverre/web/pages/nbody.htm for the latest version. * * C H A I N * ********* * * Chain regularization with slow-down and collisions. * --------------------------------------------------- * * Method of Mikkola & Aarseth, Cel. Mec. 57, 439 & 64, 197. * ......................................................... * * Developed by Seppo Mikkola & Sverre Aarseth, IOA, Cambridge. * ............................................................ * PROGRAM CHAIN * INCLUDE 'commonc.h' INCLUDE 'common2.h' LOGICAL CHECK,KSLOW,KCOLL,ICASE REAL*8 KSCH COMMON/SLOW1/ TK2(0:NMX),EJUMP,KSCH(NMX),KSLOW,KCOLL COMMON/CALLS/ TPR,TKPR,STEP,IPERT,ICALL,NFN,ITER COMMON/CCOLL/ QK(NMX4),PK(NMX4),RIJ(NMX,NMX),SIZE(NMX), & RCOLL,ECOLL,QPERI,ICOLL COMMON/INTFAC/ LX,LE,LP,LV,LT,NHALF2 COMMON/BSSAVE/ DSC,JC COMMON/KSAVE/ K1,K2 COMMON/SLOW3/ GCRIT COMMON/EXTRA/ XN(3,NMX),VN(3,NMX),RGRAV,RESC,EBIN,NAMEC(NMX) COMMON/PLUMMER/ ZMP,EPS2 REAL*8 G0(3),Y(NMX8),R2(NMX,NMX),X0(3,NMX),V0(3,NMX),XX(3,NMX) INTEGER IJ(NMX) * * * Main variables for chain regularization * *************************************** * * ------------------------------------------------------------------- * ENERGY Total energy. * EPS Tolerance for DIFSY1. * KSLOW Indicator for slow-down. * M Particle mass. * NAMEC Original particle name. * NEQ Total number of equations (8*N). * NFN Number of function calls. * NREG Number of regularization switches. * P Regularized momenta. * Q Regularized coordinates. * RCOLL Minimum two-body separation. * RESC Escape check distance (in units of RGRAV). * RGRAV Gravitational radius ((sum M(I)*M(J))/ABS(ENERGY)). * TCR Crossing time ((sum M(I))**2.5/ABS(2*ENERGY)**1.5). * TCRIT Termination time (in units of crossing time). * X Particle coordinates (X(1), X(2), X(3) is X, Y, Z). * V Velocity components. * ------------------------------------------------------------------- * * Input Parameters * ------------------------------------------------------------ * N Particle number. * TOL Tolerance for Bulirsch-Stoer integrator. * TCRIT Maximum time (in units of crossing time). * RESC Escape distance (scaled by gravitational radius). * ECRIT Error criterion (rms for time reversal check). * IREV Time reversal option (IREV = 0: standard). * JCOLL Close encounter (=1) and collision (>1) indicator. * ISLOW Slow-down indicator. * IMOVE Movie indicator (> 0: read NWAIT, XMAX, IWIN). * * ZMP Mass of Plummer model (= 0 for isolated system). * EPS2 Plummer softening (square is saved). * * NWAIT Loop index for CALL WAIT to increase viewing time. * XMAX Movie frame size (centre to edge). * IWIN Window size in pixels. * * M,R,V Mass, coordinates & velocities (routine CHINIT). * * SIZE Individual radii (JCOLL > 0; = 0 without collision). * * Template (N = 3 with movie, zero Plummer mass, no collision): * 3 1.0E-12 60.0 10.0 0.002 0 0 0 1 * 0.0 1.0 * 200000 6.0 800 * 3.0 1.0 4.0 0.0 0.0 0.0 0.0 * 4.0 -2.0 0.0 0.0 0.0 0.0 0.0 * 5.0 1.0 0.0 0.0 0.0 0.0 0.0 * 1.0D-02 1.0D-02 1.0D-02 8.0D-04 * ------------------------------------------------------------ * * Initialize the timer. CALL CPUTIM(TCOMP) * * Read and print initial parameters. READ (5,*) N, TOL, TCRIT, RESC, ECRIT, IREV, JCOLL, ISLOW, IMOVE WRITE (6,1) N, TOL, TCRIT, RESC, ECRIT, IREV, JCOLL, ISLOW, IMOVE 1 FORMAT (//,' CHAIN PARAMETERS: N =',I3,' TOL =',1P,E9.1, & ' TCRIT =',0P,F6.1,' RESC =',F6.1,' ECRIT =',F8.4, & ' IREV =',I3,' JCOLL =',I3,' ISLOW =',I3,' IMOVE =',I3) * * Read parameters for Plummer model. READ (5,*) ZMP, EPS2 EPS2 = EPS2**2 * IF (IMOVE.GT.0) THEN READ (5,*) NWAIT, XMAX, IWIN CALL MOVIE0(XMAX,IWIN) END IF * * Generate initial conditions in the c.m. rest frame. CALL CHINIT * * Save the tolerance and initial conditions. EPS = TOL LK = 0 DO 3 L = 1,N NAMEC(L) = L DO 2 K = 1,3 LK = LK + 1 X0(K,L) = X(LK) V0(K,L) = V(LK) 2 CONTINUE 3 CONTINUE * Check reading of radii. IF (JCOLL.GT.0) READ (5,*) (SIZE(K),K=1,N) * * Set large values for the two-body separations. DO 6 J = 1,10 DO 5 K = 1,10 RIJ(J,K) = 100.0 5 CONTINUE 6 CONTINUE * * Initialize diagnostic & COMMON variables. 8 NSTEP = 0 NSLOW = 0 NSTEPX = 0 NEXT = 0 NREG = 0 NFN = 0 KCASE = 0 IEND = -1 IRUN = 0 ITER = 0 IPERT = 0 IF (ZMP.GT.0.0D0) IPERT = 1 ICALL = 0 ICOLL = 0 ICASE = .FALSE. JC = 0 NZERO = N DSC = 1.0 ZKX = 1.0 RCOLL = 100.0 QPERI = 100.0 RESC0 = RESC ZMI1 = 100.0 ZMI2 = 100.0 ZMI3 = 100.0 * * Specify integration counter and initialize subsystem time. KSMAX = 50 CHTIME = 0.0D0 * Set independent dummy variable for DIFSY1. STIME = 0.0D0 * * Initialize chain regularization (new case or modified chain). 10 NEQ = 8*N Y(NEQ) = CHTIME * * Define indices for DIFSY1 calls to derivative routine. NC = N - 1 LX = 4*NC + 1 LE = 4*NC + 4 LP = 4*NC + 5 LV = 8*NC + 5 LT = 8*NC + 8 NHALF2 = (NEQ/2)*2 * * Evaluate constants of the motion. CALL CONST(X,V,M,N,ENER0,G0,ALAG) * * Select chain indices and transform from chain to KS variables. CALL SELECT CALL TRANSQ * * Define total energy (including any c.m. motion just in case). ENERGY = ENER0 - 0.5D0*MASS*(CMV(1)**2 + CMV(2)**2 + CMV(3)**2) ETOT0 = ENERGY * * Copy whole input array for the integrator. CALL YCOPY(Y) * * Find sum of mass products and set current gravitational radius. SUM = 0.0D0 DO 20 L = 1,N-1 DO 15 K = L+1,N SUM = SUM + MIJ(L,K) 15 CONTINUE 20 CONTINUE RGRAV = SUM/ABS(ENERGY) RSLOW = 0.5*RGRAV IF (ISLOW.EQ.0) RSLOW = 0.0 RSTAR = 0.1*RGRAV * * Define escaper check distance and escape criterion. RESC = RESC0*RGRAV * * Specify current local crossing time and nominal termination time. TCR = MASS**2.5/ABS(2.0D0*ENERGY)**1.5 IF (CHTIME.EQ.0.0D0) TMAX = TCRIT*TCR * * Determine the smallest two-body time-scale from parabolic orbit. CALL R2SORT(IJ,R2) I1 = IJ(1) I2 = IJ(2) RM = SQRT(R2(I1,I2)) VP2 = 2.0*(M(I1) + M(I2))/RM TP = RM/SQRT(VP2) TSTAR = MIN(TP,TCR) * * Initialize the slow-down vector and logical variables. DO 22 I = 1,N-1 KSCH(I) = 1.0D0 22 CONTINUE KSLOW = .false. KCOLL = .false. EJUMP = 0.0D0 GCRIT = 1.0D-05 * * Convert to regularized time units (note: ds = L dt is approximate). STEP = EPS**0.2*TSTAR*ALAG * WRITE (6,25) CHTIME, RGRAV, TCR, TMAX, ENERGY 25 FORMAT (/,' NEW CHAIN: T =',F7.2,' RG =',F7.3,' TCR =',F6.2, & ' TMAX =',F7.1,' ENERGY =',1P,E12.4) * * Evaluate initial binary energies. IF (N.EQ.NZERO) THEN CALL BINARY(IEND) END IF * * Perform regularized integration until termination or modification. KSTEPS = 0 30 KSTEPS = KSTEPS + 1 STEP0 = STEP TIME0 = CHTIME * * Advance solution by Bulirsch--Stoer integrator (note movie option). IF (IMOVE.GE.3) STEP = 0.5*STEP IF (IMOVE.EQ.2) STEP = 0.2*STEP IF (IMOVE.EQ.1) STEP = 0.1*STEP CALL DIFSY1(NEQ,EPS,STEP,STIME,Y) * * Set physical time and save COMMON variables. CHTIME = Y(NEQ) CALL YSAVE(Y) CALL TRANSX NSTEP = NSTEP + 1 * * Collect slow-down diagnostics if active. IF (ISLOW.GT.0) THEN ZK = 1.0 DO 32 I = 1,N-1 ZK = MAX(KSCH(I),ZK) 32 CONTINUE IF (ZK.GT.1.0D0) NSLOW = NSLOW + 1 ZKX = MAX(ZK,ZKX) END IF * * Save new step during standard integration for subsequent restart. IF (ICALL.EQ.0.AND.ICOLL.EQ.0) THEN SAVEIT = STEP ICASE = .TRUE. * Reset collision indicator and activate IPREV on failed iteration. ELSE IF (ICOLL.LT.0.AND..NOT.KCOLL) THEN ICOLL = 0 IPREV = 1 ICASE = .TRUE. KCOLL = .FALSE. NEXT = NSTEP + 2 STEP = SAVEIT * Restore the original configuration and copy to input array. TPR = TKPR DO 35 I = 1,N-1 KS = 4*(I - 1) DO 34 J = 1,4 Q(KS+J) = QK(KS+J) P(KS+J) = PK(KS+J) 34 CONTINUE 35 CONTINUE CALL YCOPY(Y) GO TO 30 END IF * * Check slow-down procedure for small minimum distance. IF (ISLOW.GT.0) then RM = 0.0 DO 40 I = 1,N-1 IF (RM.LT.RINV(I)) THEN RM = RINV(I) IM = I END IF 40 CONTINUE RM = 1.0/RM IF (RM.LT.RSLOW.AND.ICOLL.EQ.0) THEN CALL SLOW(Y) END IF END IF * * Include safety test on large or small steps. IF (STEP.GT.10.0*STEP0.AND.ICOLL.EQ.0) THEN STEP = 10.*STEP0 ELSE IF (STEP.EQ.0.0D0) THEN WRITE (6,*) ' Stepsize = 0!',char(7) STOP END IF * * Determine two-body distances for stability test and collision search. IF (CHTIME.GT.TIME0.AND.JC.EQ.0) THEN RM = 0.0 RC2 = 0.0 DO 42 K = 1,N-1 * Find minimum separation for collision and stability test. IF (RINV(K).GT.RM) THEN RM = RINV(K) IX = K END IF RC2 = RC2 + RINV(K)**2 42 CONTINUE * Update sum of 1/R^2 during forward integration (R^2 before 12/99). RM = 1.0/RM ZMI = RC2 ZMI1 = ZMI2 ZMI2 = ZMI3 ZMI3 = ZMI * Set search distance for closest separation. I1 = INAME(IX) I2 = INAME(IX+1) SX = SIZE(I1) + SIZE(I2) * Include an extra factor in search criterion. SX = 8.0*SX * Adopt more generous search criterion without collision test. IF (JCOLL.EQ.1) SX = MAX(RSTAR,SX) END IF * * Switch on search indicator during iteration or just after pericentre. IF (ICOLL.LT.0) ICALL = 1 IF (RM.LT.SX.AND.NSTEP.GT.NEXT) THEN IF (ZMI3.LT.ZMI2.AND.ZMI2.GT.ZMI1) THEN ICALL = 1 IF (ISLOW.GT.0.AND.KSLOW) THEN GSAVE = GCRIT GCRIT = 0.0 CALL SLOW(Y) GCRIT = GSAVE END IF END IF END IF * * Check for collision during the last step. IF (ICOLL.LE.0) GO TO 60 IPREV = ICOLL * * Restore the minimum configuration from DERQP. TPR = TKPR DO 50 I = 1,N-1 KS = 4*(I - 1) DO 45 J = 1,4 Q(KS+J) = QK(KS+J) P(KS+J) = PK(KS+J) 45 CONTINUE 50 CONTINUE * * Delay next search by two steps to avoid the same pericentre. NEXT = NSTEP + 2 ICALL = 0 * * Transform to physical variables. CALL TRANSX * WRITE (6,55) NAMEC(K1), NAMEC(K2), ENERGY-ECOLL, RCOLL 55 FORMAT (/,' COLLISION! NAMES ENERGY RCOLL ', & 2I4,F12.4,1P,2E10.2) * * Combine the two close bodies inelastically and reduce membership. CALL NEWSYS(K1,K2) RCOLL = 100.0 ICOLL = 0 GO TO 10 * * Include optional movie using PGPLOT or X11 graphics. 60 IF (IMOVE.GT.0) THEN * Copy chain variables to standard form. LK = 0 DO 65 L = 1,N DO 64 K = 1,3 LK = LK + 1 XX(K,L) = X(LK) 64 CONTINUE 65 CONTINUE CALL MOVIE(N,XX,CHTIME) ZZ = 0.0 DO 66 KK = 1,NWAIT ZZ = ZZ + SQRT(FLOAT(KK))/SQRT(FLOAT(KK+1)) 66 CONTINUE IF (ZZ.LT.0.0) STOP END IF * * Check switching condition (RINV now updated after switch). ISW = 0 CALL SWCOND(CHECK) IF (CHECK) THEN CALL SWITCH(Y) NREG = NREG + 1 ISW = 1 END IF * * Check termination criteria (T > TMAX or RSUM > RESC). IF ((CHTIME.GT.TMAX).OR.(KSTEPS.GT.KSMAX) & .OR.(RSUM.GT.RESC).OR.(ITER.GT.0)) THEN CALL YSAVE(Y) CALL TRANSX KSTEPS = 0 ELSE GO TO 30 END IF * * Test for time reversal iteration. IF (IRUN.GT.0) THEN IF (CHTIME.GT.TMAX.OR.ITER.GT.0) THEN IF (ABS(CHTIME - TMAX).LT.1.0E-12) GO TO 140 STEP = (TMAX - CHTIME)/TPR ITER = ITER + 1 IF (ITER.LE.9) GO TO 30 GO TO 140 END IF END IF * * Search for distant escaper or check dominant binary (RSUM > RESC). IF (RSUM.GT.RESC.AND.ISW.EQ.0) THEN CALL CONST(X,V,M,N,ETOT,G0,ALAG) CALL ESCAPE(KCASE) IF (KCASE.GT.0) THEN ETOT = ETOT - 0.5*MASS*(CMV(1)**2 + CMV(2)**2 + CMV(3)**2) WRITE (6,70) CHTIME, ETOT, (ETOT-ENERGY)/ENERGY 70 FORMAT (' CHECK: T E DE/E ',F9.2,F12.5,1P,E10.2) ETOT0 = ETOT END IF * Distinguish between N <= 2, new chain or continued integration. IF (N.LE.2) GO TO 80 IF (KCASE.GT.0) GO TO 10 IF (KCASE.LT.0) GO TO 75 IF (CHTIME.LT.TMAX) GO TO 30 ELSE * Check termination by energetic binary. IF (RSUM.GT.RESC.OR.RSUM.GT.100.0*RM) THEN CALL BINARY(IEND) IF (IEND.GT.0) GO TO 80 END IF END IF * * Examine three-body stability (relevant to decayed system). IF (RSUM.GT.4.0*RM.AND.N.EQ.3) THEN IF (NZERO.GT.3) THEN CALL CHSTAB(ITERM) IF (ITERM.LT.0) GO TO 75 END IF END IF * * Continue integration until termination time if no dominant binary. IF (N.GT.2.AND.CHTIME.LT.TMAX) GO TO 30 * 75 IF (IEND.EQ.0) THEN CALL BINARY(IEND) END IF * * Generate some relevant output. 80 TC = CHTIME/TCR CALL CPUTIM(TCOMP) CALL CONST(X,V,M,N,ETOT,G0,ALAG) ETOT = ETOT - 0.5D0*MASS*(CMV(1)**2 + CMV(2)**2 + CMV(3)**2) ERR = (ETOT - ENERGY)/ENERGY WRITE (6,85) TCOMP, TC, NSTEP, NREG, NFN, ERR, (NAMEC(K),K=1,N) 85 FORMAT (/,' END CHAIN ',' TCOMP =',F6.2,' TC =',F6.1, & ' #',I6,I4,' NFN =',I8,' ERR =',1P,E9.1, & ' NAME =',10I3) WRITE (6,90) CHTIME, NSLOW, ZKX, RCOLL, (1.0/RINV(K),K=1,N-1) 90 FORMAT (' TIME =',F7.1,' NSL =',I5,' KAPPA =',F7.2, & ' RCOLL =',1P,E9.1,' R =',10E9.1) * * Check time reversal indicator or terminate run. IF (IREV.GT.0) THEN IF (IRUN.EQ.0) THEN TMAX = CHTIME CHTIME = 0.0D0 NSTEP = 0 KSTEPS = 0 NREG = 0 IEND = 0 DO 95 K = 1,4*NC P(K) = -P(K) 95 CONTINUE DO 99 K = 1,3 CMV(K) = -CMV(K) 99 CONTINUE CALL YCOPY(Y) IRUN = IRUN + 1 END IF GO TO 30 ELSE IF (IMOVE.GT.0) THEN CALL MOVIE1 END IF STOP END IF * * Evaluate time reversal errors. 140 ERRX = 0.0 ERRV = 0.0 LK = 0 DO 150 L = 1,N DO 145 K = 1,3 LK = LK + 1 ERRX = ERRX + (X(LK) - X0(K,L))**2 ERRV = ERRV + (V(LK) + V0(K,L))**2 145 CONTINUE 150 CONTINUE ERRX = SQRT(ERRX/FLOAT(N)) ERRV = SQRT(ERRV/FLOAT(N)) * * Print diagnostic (successful run on unit #7, unsuccessful on #8). IF (MAX(ERRX,ERRV).LT.ECRIT) THEN WRITE (7,160) CHTIME, NSTEP, RCOLL, ERRX, ERRV, EPS 160 FORMAT (' REVERSE: ',F9.4,I6,1P,4E10.2) ELSE WRITE (8,160) CHTIME, NSTEP, RCOLL, ERRX, ERRV, EPS * * Reduce tolerance by 100 and try again (until limit of 1.0E-14). IF (EPS.GT.0.99E-12) THEN EPS = 0.01*EPS LK = 0 DO 170 L = 1,N DO 165 K = 1,3 LK = LK + 1 X(LK) = X0(K,L) V(LK) = V0(K,L) 165 CONTINUE 170 CONTINUE GO TO 8 END IF END IF * CALL CONST(X,V,M,N,ETOT,G0,ALAG) ETOT = ETOT - 0.5D0*MASS*(CMV(1)**2 + CMV(2)**2 + CMV(3)**2) WRITE (6,70) CHTIME, ETOT0, (ETOT-ENERGY)/ENERGY * STOP * END