escape.f

Go to the documentation of this file.

      SUBROUTINE escape
*
*
*       Escaper removal.
*       ----------------
*
      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),iflag(mmax)
      COMMON /echain/ ech
      CHARACTER*11  which1
      LOGICAL  first
      SAVE  first
      DATA  first /.true./
*
*
*       Adopt twice the tidal radius as escape condition.
      resc2 = 4.0*rtide**2
      rtide2 = rtide**2
      ncorr = 0
      ncrit1 = 0
      i3hi = 0
      DO 1 k = 1,3
          cmr(k) = 0.0d0
          cmrdot(k) = 0.0d0
    1 CONTINUE
*
*       Set the distance (squared) with respect to the density centre.
      i = ifirst
    5 ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 + 
     &                               (x(3,i) - rdens(3))**2
      IF (ri2.LT.rtide2) ncrit1 = ncrit1 + 1
      vi2 = xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
*
*       See whether escape is indicated (retain merger ghost particles).
      IF (ri2.GT.resc2.AND.ri2.LT.1.0d+10) THEN
          rjmin2 = 1.0d+10
*       Find distance to the nearest neighbour and calculate potential.
          poti = 0.0d0
          DO 8 j = ifirst,ntot
              IF (j.EQ.i) go to 8
              rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
     &                                      (x(3,i) - x(3,j))**2
              poti = poti + body(j)/sqrt(rij2)
              IF (rij2.LT.rjmin2) THEN
                  rjmin2 = rij2
                  jmin = j
              END IF
    8     CONTINUE
*       Check escape criterion for external fields or isolated system.
          ei = 0.5*vi2 - poti
          IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
              ei = ei - mp/sqrt(ri2 + ap2)
              IF (body(i).EQ.0.0d0) go to 30
          END IF
          IF ((kz(14).GT.0.AND.kz(14).NE.4).OR.ei.GT.0.0) go to 30
      END IF
*
   10 i = i + 1
   12 IF (n.EQ.2) go to 13
*
      IF (i.LE.ntot) go to 5
      IF (ncorr.EQ.0) go to 25
*
*       Form centre of mass terms.
   13 DO 15 i = 1,n
          DO 14 k = 1,3
              cmr(k) = cmr(k) + body(i)*x(k,i)/zmass
              cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)/zmass
   14     CONTINUE
   15 CONTINUE
      cmr(4) = sqrt(cmr(1)**2 + cmr(2)**2 + cmr(3)**2)
*
      stepi = min(stepi,1.0d0)
      jlast = min(ncorr,nmax)
      WRITE (6,18)  n, nsesc, nbesc, zmass, be(3), cmr(4), resc, stepi,
     &              rsi, zmass/float(n), ncrit1, (jlist(j),j=1,jlast)
   18 FORMAT (/,' ESCAPE    N =',2i6,i4,f8.4,f12.6,2f7.2,f7.3,f7.2,f9.5,
     &                                        i6,2x,6i6,/,5(6x,20i6,/))
*
*     IF (KZ(23).EQ.2.AND.KZ(14).EQ.1) THEN
*         JLAST = MIN(2*NCORR,LMAX)
*         WRITE (6,20)  (ILIST(J),J=1,JLAST)
*  20     FORMAT (/,' ESCAPE ANGLES ',11(2I4,2X),9(/,15X,11(2I4,2X)))
*     END IF
*
*       Check updating of global index for chain c.m.
      IF (nch.GT.0) THEN
          CALL chfind
      END IF
*
*       Set phase indicator < 0 for new time-step list in INTGRT.
      iphase = -1
*
   25 RETURN
*
   30 a2 = (x(1,jmin) - rdens(1))**2 + (x(2,jmin) - rdens(2))**2 +
     &                                 (x(3,jmin) - rdens(3))**2
*       See whether nearest body satisfies escape condition or RIJ > 10*RMIN.
      IF (a2.GT.resc2.AND.kz(14).GT.0) go to 40
      IF (rjmin2.GT.100.0*rmin2) go to 40
      a3 = xdot(1,jmin)**2 + xdot(2,jmin)**2 + xdot(3,jmin)**2
      a4 = (xdot(1,i) - xdot(1,jmin))**2 +
     &     (xdot(2,i) - xdot(2,jmin))**2 +
     &     (xdot(3,i) - xdot(3,jmin))**2
      a5 = (body(i) + body(jmin))/sqrt(rjmin2)
*       Check velocity of binary component in case of bound pair.
      a6 = 2.0*a5/a4
      IF (a6.GT.1.0.AND.a3.GT.2.0*zmass/sqrt(a2)) go to 40
*       Retain escaper if dynamical effect on neighbour is significant.
      IF (a6.GT.0.01) go to 10
*
*       Form optional output diagnostics.
   40 x1 = x(1,i) - rdens(1)
      y1 = x(2,i) - rdens(2)
      z1 = x(3,i) - rdens(3)
      a3 = abs(y1/x1)
      a4 = abs(z1)/sqrt(x1**2 + y1**2)
      ilist(2*ncorr+1) = datan(a3)*180.0/3.14159
      ilist(2*ncorr+2) = datan(a4)*180.0/3.14159
*       Escape angles with respect to the X-axis and XY-plane.
      resc = sqrt(ri2)
      stepi = step(i)
      rsi = rs(i)
*
*       Accumulate escaper names and save current name in case I > N.
      ncorr = ncorr + 1
      jlist(ncorr) = name(i)
      namei = name(i)
      kstari = kstar(i)
*       Include termination for escaping chain c.m. system (mainly NBODY7).
      IF (name(i).EQ.0) THEN
          nsub = max(nsub - 1,0)
          nch = 0
          WRITE (6,42)  ech, ei
   42     FORMAT (' CHAIN ESCAPE!!!!!    ECH EI  ',2f10.6)
          ecoll = ecoll + ech
          ech = 0.0
      END IF
      IF (name(i).LT.0) kstari = kstari + 30
      IF (body(i).LE.0.0d0) go to 50
*
*       Check initialization of tidal tail integration for 3D model.
      IF (kz(23).GE.3.AND.kz(14).EQ.3) THEN
          CALL tail0(i)
      END IF
*
*       Obtain binding energy of body #I and update optional tidal radius.
      zk = 0.5d0*body(i)*vi2
      IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
          CALL xtrnlv(i,i)
          zk = zk + ht
          rtide = (zmass/tidal(1))**0.3333
      ELSE IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
          rtide = rtide0*zmass**0.3333
      END IF
      ei = zk - body(i)*poti
*
*       Correct total energy (also ZKIN & POT for consistency).
      be(3) = be(3) - ei
      zkin = zkin - zk
      pot = pot - body(i)*poti
*
*       Update total mass and save energy & number of single escapers.
      IF (i.LE.n) THEN
          zmass = zmass - body(i)
          e(4) = e(4) + ei
          npop(4) = npop(4) + 1
          nsesc = nsesc + 1
      ELSE
          ipair = i - n
      END IF
      esesc = esesc + ei
*
*       Include optional escape output on unit 11.
      IF (kz(23).GT.1) THEN
          IF (first) THEN
              OPEN (unit=11,status='NEW',form='FORMATTED',file='ESC')
              first = .false.
          END IF
          tesc = tscale*ttot
          eesc = ei/body(i)
          vb2 = 2.0*zkin/zmass
*       Distinguish between tidal field and isolated system (ECRIT at RTIDE).
          IF (kz(14).GT.0.AND.kz(14).LE.2) THEN
              ecrit = -1.5*(tidal(1)*zmass**2)**0.3333
              eesc = 2.0*(eesc - ecrit)/vb2
              eesc = min(eesc,999.0d0)
              besc = zmbar*body(i)
          ELSE
              eesc = 2.0*eesc/vb2
              besc = body(i)/bodym
          END IF
          vkm = sqrt(vi2)*vstar
          WRITE (11,45)  tesc, besc, eesc, vkm, kstari, name(i)
   45     FORMAT (f8.1,3f6.1,i4,i7)
          CALL flush(11)
      END IF
*
*       Reduce particle number & total membership and check NNBMAX.
   50 n = n - 1
      ntot = ntot - 1
      nnbmax = min((n - npairs)/2,nnbmax)
      znbmax = 0.9*nnbmax
      znbmin = max(0.2*nnbmax,1.0)
*       Set indicator for removal of c.m. or KS components.
      kcase = 1
*
*       Remove any circularized KS binary from the CHAOS/SYNCH table.
      IF (i.GT.n + 1.AND.kstar(i).GE.10) THEN
          ii = -(n + 1 + ipair)
          CALL spiral(ii)
      END IF
*
*       Update COMMON arrays to remove the escaper and correct F & FDOT.
      CALL remove(i,1)
*
*       Delete escaper from neighbour lists and reduce higher locations.
   60 DO 150 j = 1,ntot
          nnb = list(1,j)
          IF (nnb.EQ.0) go to 150
          l = 2
   70     IF (list(l,j).NE.i) go to 130
*
*       Move up the remaining list members and reduce membership by one.
          DO 80 k = l,nnb
              list(k,j) = list(k+1,j)
   80     CONTINUE
          list(1,j) = list(1,j) - 1
*       Reduce the steps to minimize error effect (do not allow DT < 0).
*         STEP(J) = MAX(0.5D0*STEP(J),TIME - T0(J))
*         STEPR(J) = MAX(0.5D0*STEPR(J),TIME - T0R(J))
          IF (list(1,j).GT.0) go to 130
*
*       Add a distant body as neighbour if list only contains escaper.
          k = ifirst - 1
  100     k = k + 1
          rjk2 = (x(1,j) - x(1,k))**2 + (x(2,j) - x(2,k))**2 +
     &                                  (x(3,j) - x(3,k))**2
          IF (rjk2.LT.resc2.AND.k.LT.n) go to 100
          list(1,j) = 1
          list(2,j) = k
          go to 150
*
*       Reduce higher particle locations by one.
  130     IF (list(l,j).GT.i) list(l,j) = list(l,j) - 1
          l = l + 1
          IF (l.LE.list(1,j) + 1) go to 70
  150 CONTINUE
*
*       Update list of old KS components (remove #I and rename > I).
      nnb = listr(1)
      DO 170 l = 2,nnb+1
          IF (listr(l).EQ.i) THEN
*       Remove both components of pair and reduce membership by two.
              j = 2*kvec(l-1)
              DO 165 k = j,nnb-1
                  listr(k) = listr(k+2)
  165         CONTINUE
              listr(1) = listr(1) - 2
          END IF
  170 CONTINUE
*
*       Reduce higher particle locations by one (separate loop for pairs).
      DO 180 l = 2,nnb+1
          IF (listr(l).GT.i) listr(l) = listr(l) - 1
  180 CONTINUE
*
*       Update list of high velocity particles (remove #I and rename > I).
      nnv = listv(1)
      DO 190 l = 2,nnv+1
          IF (listv(l).EQ.i) THEN
*       Remove escaper and reduce the membership.
              DO 185 k = l,nnv
                  listv(k) = listv(k+1)
  185         CONTINUE
              listv(1) = listv(1) - 1
          END IF
*       Reduce higher particle locations by one (or three for c.m.).
          IF (listv(l).GT.i) THEN
              listv(l) = listv(l) - 1
              IF (i.GT.n + 1) listv(l) = listv(l) - 2
          END IF
  190 CONTINUE
*
*       Check special case of second KS component removal.
      IF (kcase.GT.1) go to 205
*
*       See whether the escaper is a single particle or c.m.
      IF (i.LE.n + 1) go to 12
*
*       Skip correction if ghost is also merged binary (NAMEI = 0 below).
      IF (namei.NE.0.AND.body(2*ipair-1).GT.0.0d0) THEN
          zmb = body(2*ipair-1) + body(2*ipair)
          zm2 = body(2*ipair)
          eb = h(ipair)*body(2*ipair-1)*body(2*ipair)/zmb
          semi = -0.5*zmb/h(ipair)
          ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(semi*zmb)
          ecc = sqrt(ecc2)
          pmin  = semi*(1.0 - ecc)
          ratio = max(radius(2*ipair-1),radius(2*ipair))/pmin
          name1 = name(2*ipair-1)
          kstar1 = kstar(2*ipair-1)
          pb = days*semi*sqrt(semi/zmb)
      ELSE
*       Obtain two-body elements of ghost binary and update energies.
          eb = 0.0d0
          ratio = 0.0
          bodycm = cm(3,im) + cm(4,im)
          semi = -0.5*bodycm/h(jpair)
          zmu = cm(3,im)*cm(4,im)/bodycm
          pb = days*semi*sqrt(semi/bodycm)
          ecc2 = (1.0-r(jpair)/semi)**2 + tdot2(jpair)**2/(bodycm*semi)
          ecc = sqrt(ecc2)
          eb1 = zmu*h(jpair)
          be(3) = be(3) - eb1
          emerge = emerge - eb1
          deb = deb + eb1
      END IF
*
*       Update total energy (ECOLL with EB < -1 & #27 > 0 affects BINOUT).
      be(3) = be(3) - eb
      ebin = ebin - eb
*
*       Check optional diagnostics for hierarchical systems.
      IF (namei.LE.0.AND.(kz(18).EQ.1.OR.kz(18).EQ.3)) THEN
          iphase = 7
          CALL hiarch(ipair)
      END IF
*
*       Distinguish between actual and ghost binary (mergers come later).
      IF (body(2*ipair-1).GT.0.0d0) THEN
          name2 = name(2*ipair)
*
*       Include rare case of higher-order system (4, 5 or 6 members).
          IF (namei.LT.-2*nzero) THEN
              jm = 1
              DO 195 k = 1,nmerge
                  IF (namem(k).EQ.namei) jm = k
  195         CONTINUE
              i3hj = n
              DO 196 j = ifirst,ntot
                  IF (name(j).EQ.nameg(jm)) i3hj = j
  196         CONTINUE
*       Identify quartet, quintuplet or sextuplet.
              IF (name2.LE.nzero.AND.name(i3hj).LE.nzero) THEN
                  which1 = '  QUARTET  '
*       Include both types of quintuplet: [[B,S],B] and [[B,B],S].
              ELSE IF (name2.LE.nzero.OR.name(i3hj).LE.nzero) THEN
                  which1 = ' QUINTUPLET'
              ELSE
                  which1 = ' SEXTUPLET '
              END IF
              zm1 = cm(1,jm)*zmbar
              zm2 = cm(2,jm)*zmbar
              zm3 = (cm(3,jm) + cm(4,jm))*zmbar
              eb3 = cm(1,jm)*cm(2,jm)*hm(jm)/(cm(1,jm) + cm(2,jm))
              semi3 = -0.5*(cm(1,jm) + cm(2,jm))/hm(jm)
              pb3 = days*semi3*sqrt(semi3/(cm(1,jm) + cm(2,jm)))
              WRITE (6,199)  which1, name(2*ipair-1), name(2*ipair),
     &                       name(i3hj), zm1, zm2, zm3, eb3, semi3, pb3
  199         FORMAT (/,a11,' ESCAPE    NM =',3i6,'  M =',3f6.2,
     &                    '  EB3 =',f8.4,'  A3 =',1p,e8.1,'  P3 =',e8.1)
*       Copy actual name of outer component for KS binary output.
              name2 = name(i3hj)
          END IF
*
          vi = sqrt(0.5*vi2*zmass/abs(zkin))
          WRITE (6,200)  ipair, name(2*ipair-1), name2,
     &                   kstar(2*ipair-1), kstar(2*ipair), kstari,
     &                   list(2,2*ipair), body(2*ipair-1)*zmbar,
     &                   body(2*ipair)*zmbar, eb, ratio, vi, ecc, ei, pb
  200     FORMAT (/,' BINARY ESCAPE    KS =',i5,'  NM =',2i6,
     &                '  K* =',4i3,'  M =',2f5.1,'  EB =',f8.4,
     &                '  R*/PM =',f5.2,'  V/<V> =',f5.2,'  E =',f5.2,
     &                '  EI =',f8.5,'  P =',1p,e8.1)
      ELSE
          WRITE (6,202)  ipair, name(2*ipair-1), name(2*ipair),
     &                   kstar(2*ipair-1), kstar(2*ipair), kstar(i),
     &                   cm(3,im)*smu, cm(4,im)*smu, eb1, ecc, semi, pb
  202     FORMAT (/,' QUAD ESCAPE    KS =',i5,'  NM =',2i6,
     &                '  K* =',3i3,'  M =',2f5.1,'  EB =',f8.4,
     &                '  E =',f7.3,'  A =',1p,e8.1,'  P =',e8.1)
          eb = eb + eb1
          ei = 0.0
      END IF
*
*       Accumulate escaping binary energies and increase the counter.
      IF (list(2,2*ipair).EQ.-1) THEN
          e(5) = e(5) + eb
          e(6) = e(6) + ei
          npop(5) = npop(5) + 1
      ELSE
          e(7) = e(7) + eb
          e(8) = e(8) + ei
          npop(6) = npop(6) + 1
      END IF
      ebesc = ebesc + eb
*
      IF (namei.GT.0) THEN
          nbesc = nbesc + 1
      ELSE
          nmesc = nmesc + 1
      END IF
*
*       Reduce particle number, pair index & single particle index. 
      n = n - 1
      npairs = npairs - 1
      ifirst = 2*npairs + 1
*
*       Move up all tables of regularized pairs below IPAIR.
      IF (ipair.LE.npairs) THEN
          CALL remove(ipair,2)
      END IF
*
*       Increase index for removing KS components.
  205 kcase = kcase + 1
      IF (kcase.LE.3) THEN
*       Remove COMMON arrays of the second component before the first.
          i = 2*ipair + 2 - kcase
          ntot = ntot - 1
*       Reduce NTOT by 3 and N by 2 when KS pair escapes.
          CALL remove(i,3)
          go to 60
      END IF
*
*       Check selection of possible ghost in higher-order system.
      IF (i3hi.GT.0) THEN
          i = 0
          DO 208 j = ifirst,ntot
              IF (name(j).EQ.nameg(jm)) THEN
                  i = j
              END IF
  208     CONTINUE
          IF (i.GT.0) THEN
              i3hi = 0
              go to 50
          END IF
      END IF
*
*       Include the case of escaping merger.
      IF (namei.GE.0) go to 250
*
*       Locate current position in the merger table (standard case).
      im = 0
      DO 210 k = 1,nmerge
          IF (namem(k).EQ.namei) im = k
  210 CONTINUE
*       Skip on failed detection just in case.
      IF (im.EQ.0) go to 250
*
*       Include case of higher-order system (outer single or binary).
      deb = 0.0
      IF (namei.LT.0) THEN
*       Determine the ghost index.
          i3hi = 0
          DO 215 j = ifirst,ntot
              IF (name(j).EQ.nameg(im)) i3hi = j
  215     CONTINUE
*       Specify nominal escape distance for ghost removal.
          x(1,i3hi) = 1.0d+04
*       Define possible KS pair index for quadruple system correction.
          jpair = i3hi - n
      END IF
*
*       Consider current ghost unless deeper hierarchy is present.
      jm = im
      IF (i3hi.GT.0) THEN
*       Search for c.m. name one level below current (NAMEI < 0).
          IF (namei.LT.-2*nzero) THEN
              DO 220 k = 1,nmerge
                  IF (namem(k).EQ.namei + 2*nzero) jm = k
  220         CONTINUE
*       Use previous merger index to look for binary ghost at earlier level.
              i30 = i3hi
              DO 225 j = ifirst,ntot
                  IF (name(j).EQ.nameg(jm)) i3hi = j
  225         CONTINUE
              IF (i3hi.GT.n) THEN
                  jpair = i3hi - n
              ELSE
*       Adopt original solution on failure to identify binary.
                  i3hi = i30
              END IF
*       Set nominal escape distance of any new ghost (I3HI <= N possible).
              x(1,i3hi) = 1.0d+04
          END IF
      END IF
*
*       Copy merger energy to respective energy bins (primordial or new).
      zmu = cm(1,jm)*cm(2,jm)/(cm(1,jm) + cm(2,jm))
      IF (min(name1,nameg(jm)).LE.2*nbin0) THEN
          e(5) = e(5) + zmu*hm(jm) + deb
      ELSE
          e(7) = e(7) + zmu*hm(jm) + deb
      END IF
      emesc = emesc + zmu*hm(jm)
*
*       Reduce membership if IM is last (otherwise done in RESET).
      IF (im.EQ.nmerge) THEN
          nmerge = nmerge - 1
      END IF
*
*       Identify merged ghost particle (single body or binary c.m.).
      jcomp = -1
      DO 230 j = ifirst,ntot
          IF (body(j).EQ.0.0d0.AND.name(j).EQ.nameg(im)) jcomp = j
  230 CONTINUE
*       Include search over lower level on failed identification.
      IF (jcomp.EQ.-1.AND.i3hi.GT.0) THEN
          DO 232 j = ifirst,ntot
              IF (body(j).EQ.0.0d0.AND.name(j).EQ.nameg(jm)) jcomp = j
  232     CONTINUE
      END IF
*
*       Skip if correct ghost not identified (note I3HI # JCOMP if JM # IM).
      IF (jcomp.GT.0) THEN
          i = i3hi
*       Form two-body elements and period of inner binary.
          namei = 0
          npop(7) = npop(7) + 1
          bodycm = cm(1,jm) + cm(2,jm)
          semi0 = -0.5*bodycm/hm(jm)
          zmu = cm(1,jm)*cm(2,jm)/bodycm
          eb = zmu*hm(jm)
          pb = days*semi0*sqrt(semi0/bodycm)
          be(3) = be(3) - eb
          emerge = emerge - eb
*       Include extra corrections for mergers between binary pairs.
          IF (jm.LT.im) THEN
              zmu2 = cm(1,im)*cm(2,im)/(cm(1,im) + cm(2,im))
              eb2 = zmu2*hm(im)
              be(3) = be(3) - eb2
              emerge = emerge - eb2
              IF (jm.EQ.nmerge) nmerge = nmerge - 1
          END IF
          rj = 0.0
          td2 = 0.0
          DO 235 k = 1,4
              rj = rj + um(k,jm)**2
              td2 = td2 + 2.0*um(k,jm)*umdot(k,jm)
  235     CONTINUE
          ecc2 = (1.0 - rj/semi0)**2 + td2**2/(bodycm*semi0)
          ecc0 = sqrt(ecc2)
          pcrit = stability(cm(1,jm),cm(2,jm),body(jcomp),ecc0,ecc,
     &                                                       0.0d0)
          pcrit = pcrit*semi0
*
          WRITE (6,240)  name1, nameg(jm), kstar1, kstar(jcomp),
     &                   kstarm(jm), cm(1,jm)*zmbar, cm(2,jm)*zmbar,
     &                   eb, ecc0, pmin/pcrit, semi0, pb
  240     FORMAT (/,' HIARCH ESCAPE    NM =',2i6,'  K* =',3i3,
     &              '  M =',2f5.1,'  EB =',f8.4,'  E =',f7.3,
     &              '  PM/PC =',1p,e8.1,'  A =',e8.1,'  P =',e8.1)
*
*       Remove the ghost particle (NAME = 0 & EB = 0 for second binary).
          go to 50
      END IF
*
  250 i = ipair + n
      go to 12
*
      END