decide.f

Go to the documentation of this file.

      SUBROUTINE decide(IPAIR,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 #JCOMP (KS resolved in INDUCE).
      CALL edot(i1,i2,jcomp,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(IPAIR,JCOMP,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(JCOMP), 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