https://github.com/florentrenaud/nbody6tt
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
decide.f
SUBROUTINE DECIDE(IPAIR,JCL,SEMI,ECC,EMAX,EMIN,TC,TG,EDAV,IQ)
*
*
* Merger decision.
* ----------------
*
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)
*
*
* Initialize the merger skip indicator and define KS indices.
IQ = 0
I1 = 2*IPAIR - 1
I2 = I1 + 1
* Prevent over-writing in NMERGE+1.
IF(NMERGE.GE.MMAX-1)THEN
IQ = 1
GOTO 40
ENDIF
*
* Accept active spiral but ensure some updating to increase pericentre.
DMR = 0.D0
IF (KSTAR(N+IPAIR).EQ.-2) THEN
PMIN = SEMI*(1.0 - ECC)
ICIRC = -1
CALL TCIRC(PMIN,ECC,I1,I2,ICIRC,TC1)
* Restrict circularization time to account for stellar evolution.
TC1 = MIN(TC1,500.0D0)
* Adopt termination interval 1/2 of predicted t_{circ} (at const R*).
DT1 = 0.5*(TC1 + 0.1)/TSTAR
IF (ECC.LT.0.1.AND.MAX(KSTAR(I1),KSTAR(I2)).LE.1) THEN
DT1 = 10.0*DT1
END IF
* Rectify SPIRAL just in case.
IF(LIST(1,I1).GT.0)THEN
* IQ = 1
* GO TO 30
* Accept perturbed spiral case for now!
DMR = -1.D0
CALL CHRECT(IPAIR,DMR)
ELSE
CALL CHRECT(IPAIR,DMR)
ENDIF
GO TO 30
END IF
*
* Obtain eccentricity derivative due to #JCL (KS resolved in INDUCE).
CALL EDOT(I1,I2,JCL,SEMI,ECC,ECCDOT)
*
* Define c.m. index and set maximum stellar radius and peri.
I = N + IPAIR
RM = MAX(RADIUS(I1),RADIUS(I2))
PM = SEMI*(1.0 - ECC)/RM
PM = MIN(PM,99.9D0)
*
* WRITE (75,1) KSTAR(I1), KSTAR(I2), KSTAR(I), LIST(1,I1), NAME(I),
* & TTOT, ECC, EMIN, EMAX, TG, EDAV, PM
* 1 FORMAT (' DECIDE: K* NP NM T E EM EX TG ED PM ',
* & 4I4,I6,F10.3,F8.4,2F6.2,2F7.2,F6.1)
* CALL FLUSH(75)
*
* See whether to deform binary orbit (skip chaos case).
IF (KSTAR(I).NE.-1.AND.TG.LT.-10.0) THEN
* Change eccentricity slowly by a fraction 0.1*(1 - E) of decay time.
ECC0 = ECC
DT1 = 0.1*(1.0 - ECC)*TG/TSTAR
ECC = ECC + EDAV*DT1
ECC = MAX(ECC,EMIN)
* Impose temporary limit of 0.99*EMAX to allow for apsidal motion.
ECC = MIN(ECC,0.99*EMAX)
*
* Evaluate t_{circ} and (de/dt)_{circ} except for near-collision.
PMIN = SEMI*(1.0 - ECC)
IF (PMIN.GT.RM) THEN
ICIRC = -1
CALL TCIRC(PMIN,ECC,I1,I2,ICIRC,TC1)
* CALL EDOT(I1,I2,JCL,SEMI,ECC,ECCDOT)
CALL ECIRC(PMIN,ECC,I1,I2,ICIRC,TG,TC2,ECC2,EDT)
* WRITE (6,10) ECCDOT, R(IPAIR), SEMI, TDOT2(IPAIR)
* 10 FORMAT (' ECIRC: ECDOT R A TD2 ',1P,4E10.2)
ELSE
TC1 = TG
TC2 = TC1
EDT = 2.0*EDAV
END IF
*
* Compare eccentricity derivatives and check t_{circ} (standard case).
IF (ABS(EDT).GT.ABS(EDAV).AND.TC1.LT.50.0) THEN
* Enforce new spiral on short circularization time (IQ > 0: no merger).
IF (KSTAR(I).GE.0) IQ = 1
ECC = ECC0
IF (KSTAR(I).EQ.-2) GO TO 40
GO TO 30
END IF
*
* Rectify spiral orbit before changing eccentricity.
IF (KSTAR(I).EQ.-2) THEN
CALL CHRECT(IPAIR,DMR)
END IF
*
* Transform to exact apocentre unless done already.
IF (ABS(TDOT2(IPAIR)).GT.1.0D-12) THEN
IF (R(IPAIR).GT.SEMI.AND.TDOT2(IPAIR).LT.0.0) THEN
CALL KSAPO(IPAIR)
END IF
CALL KSPERI(IPAIR)
CALL KSAPO(IPAIR)
* Check spiral case already at peri.
ELSE IF (KSTAR(I).EQ.-2.AND.R(IPAIR).LT.SEMI.AND.
& ABS(TDOT2(IPAIR)).LT.1.0D-12) THEN
CALL KSAPO(IPAIR)
END IF
*
* Deform orbit (H = const at apo), rectify and transform to X & XDOT.
CALL DEFORM(IPAIR,ECC0,ECC)
CALL KSRECT(IPAIR)
CALL RESOLV(IPAIR,1)
*
* Rectify spiral orbit (terminate on collision).
IF (KSTAR(I).EQ.-2) THEN
CALL CHRECT(IPAIR,DMR)
IF (IPHASE.LT.0) THEN
IQ = 2
GO TO 40
END IF
END IF
*
* Re-initialize KS polynomials for perturbed motion (also in DEFORM).
IF (LIST(1,I1).GT.0) THEN
IF (R(IPAIR).LT.SEMI.AND.
& ABS(TDOT2(IPAIR)).LT.1.0D-12) THEN
CALL KSAPO(IPAIR)
END IF
T0(I1) = TIME
IMOD = KSLOW(IPAIR)
CALL KSPOLY(IPAIR,IMOD)
TDOT2(IPAIR) = -1.0D-20
END IF
*
* Produce diagnostic output for large eccentricity.
IF (ECC.GT.0.9) THEN
WRITE (75,20) NAME(I1), TTOT, ECC0, ECC, EMIN, EMAX,
& ECCDOT, EDT, TG, TC1, EDAV, PM
20 FORMAT (' DECIDE: NM T E0 E1 EM EX ED EDT TG TC EDAV ',
& 'PM ', I5,F9.2,4F7.3,1P,6E9.1)
CALL FLUSH(75)
END IF
ELSE
DT1 = 10.0/TSTAR
END IF
*
* See whether growth time exceeds merger disruption time from TSTAB.
30 TG1 = TIME + DT1
IF (IQ.EQ.0) THEN
TMDIS(NMERGE+1) = MIN(TG1,TMDIS(NMERGE+1))
* WRITE (6,55) TIME, NAME(I1), NAME(JCL), ECC, EMAX,DT1, TG, SEMI
* 55 FORMAT (' WATCH! T NM NMJ E EX DT1 TG A ',
* & F10.4,2I7,2F8.4,1P,3E10.2)
END IF
* WRITE (6,35) LIST(1,I1), ECC, EDAV*DT1, EMAX, TG1+TOFF, GAMMA(IPAIR)
* 35 FORMAT (' DECIDE: NP E DE EMAX TG1 G ',I4,3F7.3,F10.3,1P,E9.1)
*
40 RETURN
*
END