https://github.com/florentrenaud/nbody6tt
Raw File
Tip revision: 8aea50c213fd132d500c415511ae1e27eeabab80 authored by florent on 14 February 2015, 16:38:53 UTC
corrected mb typo in ttgalaxy
Tip revision: 8aea50c
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
back to top