https://github.com/florentrenaud/nbody6tt
Tip revision: 8aea50c213fd132d500c415511ae1e27eeabab80 authored by florent on 14 February 2015, 16:38:53 UTC
corrected mb typo in ttgalaxy
corrected mb typo in ttgalaxy
Tip revision: 8aea50c
tstab.f
SUBROUTINE TSTAB(I,ECC1,SEMI1,PMIN1,PCRIT1,YFAC,JCL,ITERM)
*
*
* Hierarchical stability time.
* ----------------------------
*
INCLUDE 'common6.h'
COMMON/BINARY/ CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
& HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
& NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
*
*
* Distinguish between absolute and approximate stability.
ITERM = 0
IF (PMIN1.GT.YFAC*PCRIT1.OR.KZ(19).EQ.0.OR.KZ(27).EQ.0) THEN
TMDIS(NMERGE+1) = 1.0D+10
ELSE IF (NAME(I).LT.0.OR.NAME(JCL).LT.0) THEN
* Skip double hierarchy (RESET2 uses standard stability criterion).
ITERM = 1
ELSE IF (PMIN1.GT.0.8*YFAC*PCRIT1.OR.
& (ECC1.GT.0.99.AND.ECC1.LT.1.0)) THEN
* Estimate time-scale for long-lived hierarchy.
NK = 1 + 10.0*ECC1/(1.0 - ECC1)
* Define termination time in terms of outer periods.
TK = TWOPI*SEMI1*SQRT(ABS(SEMI1)/(BODY(I) + BODY(JCL)))
* Set disruption time in new merger location for routine KSINT.
TMDIS(NMERGE+1) = TIME + ABS(NK*TK)
* WRITE (6,3) YFAC, PMIN1, PCRIT1, 0.8*YFAC*PCRIT1, NK*TK
* 3 FORMAT (' TSTAB: YFAC PMIN1 PCR PTEST N*TK ',F6.2,1P,4E9.2)
* Note that large NK implies termination by other means.
IF (ECC1.LT.1.0.AND.NK.LE.1) ITERM = 1
* Modify PCRIT1 so it becomes < SEMI1*(1 - ECC1)*(1 - 2*PERT).
IF (PCRIT1.LT.PMIN1) PCRIT1 = 0.999*PMIN1/YFAC
ELSE IF (ECC1.LT.1.0) THEN
* Specify termination in all other cases (but allow hyperbolic orbit).
ITERM = 1
END IF
*
RETURN
*
END