merge2.f

Go to the documentation of this file.

      SUBROUTINE merge2
*
*
*       Merging of double hierarchy.
*       ----------------------------
*
      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)
      CHARACTER*11  which1
*
*
*       COMMON variables for merge procedure
*       ************************************
*
*       ------------------------------------------------------------------
*       CM      Component masses of merged binary (2nd binary in CM(3&4)).
*       HM      Binding energy per unit mass of merged binary.
*       KSTARM  Roche stage indicator (KSTAR for 1st c.m.; ghost is OK).
*       NAMEG   Name of the associated ghost component.
*       NAMEM   Name of new c.m. (< 0) for identification of merger index.
*       NMERG   Total number of mergers.
*       NMERGE  Current number of merged binaries (maximum is MMAX).
*       UM      Regularized coordinates of merged binary.
*       UMDOT   Regularized velocity of merged binary.
*       VREL    Relative velocity of merged binary components.
*       XREL    Relative coordinates of merged binary components.
*       ------------------------------------------------------------------
*
*
*       Ensure that the hierarchy is retained as primary KS pair.
      IF (name(n+kspair).GT.0) THEN
          k = kspair
          kspair = jcomp - n
          jcomp = n + k
      END IF
*
*       Select deepest level as primary in case of two hierarchies.
      IF (name(n+kspair).LT.0.AND.name(jcomp).LT.0) THEN
*       Note possibility of a quartet (or even quintet) and triple.
          IF (name(jcomp).LT.name(n+kspair)) THEN
              k = kspair
              kspair = jcomp - n
              jcomp = n + k
          END IF
      END IF
*
*       Set pair index & c.m. of inner binary and save outer component.
      ipair = kspair
      i = n + ipair
      jcomp1 = jcomp
      icomp1 = i
      i1 = 2*ipair - 1
      i2 = i1 + 1
      semi = -0.5*body(i)/h(ipair)
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ecc = sqrt(ecc2)
*
*       Produce diagnostics for [[B,S],B] quintuplet or sextuplet system.
      IF (name(jcomp).GT.nzero) THEN
          ri2 = 0.0
          DO 2 k = 1,3
              ri2 = ri2 + (x(k,i) - rdens(k))**2
    2     CONTINUE
          ri = sqrt(ri2)
          eb = 0.5*body(i1)*body(i2)/semi
          eb = eb*float(n-npairs)/zkin
          eb = min(eb,999.9d0)
          zm = (body(i) + body(jcomp))*smu
          CALL findj(i,jg,im)
*       Check possible extension of next look-up time (skip GR case).
          IF (kz(19).GE.3.AND.kz(27).LT.3) THEN
              jlist(1) = i1
              jlist(2) = i2
              jlist(3) = 2*(jcomp - n) - 1
              jlist(4) = jlist(3) + 1
              jlist(5) = jg
              CALL newtev(5,ix)
          END IF
          which1 = ' QUINTUPLET'
          IF (jg.GT.n) which1 = ' SEXTUPLET '
          WRITE (6,5)  which1, time+toff, zm, name(i1), name(i2),
     &                 name(jcomp), ri, ecc, eb, semi,pcrit,gamma(ipair)
    5     FORMAT (/,' NEW',a11,'   T MT NM1 NM2 NM3 RI E0 EB0 A0 PC G0',
     &                             f10.2,f6.2,3i6,2f6.2,f6.1,1p,3e10.2)
      END IF
*
*       Check for [[B,S],S] quartet or [[B,B],S] quintuplet.
      IF (name(jcomp).GT.0.AND.name(jcomp).LE.nzero) THEN
          CALL findj(i,jg,im)
          IF (kz(19).GE.3.AND.kz(27).LT.3) THEN
              jlist(1) = i1
              jlist(2) = i2
              jlist(3) = jcomp
              jlist(4) = jg
              CALL newtev(4,ix)
          END IF
          zm = (body(i) + body(jcomp))*smu
          IF (jg.LE.n.AND.jcomp.LE.n) THEN
              WRITE (6,7)  time+toff, zm, name(i1), name(i2),
     &                     name(jcomp), ecc, semi, pcrit, gamma(ipair)
    7         FORMAT (/,' NEW QUARTET    T MT NM1 NMG NM3 E0 A0 PC G0 ',
     &                                   f9.2,f6.2,3i6,f6.2,1p,3e10.2)
          ELSE
              WRITE (6,8)  time+toff, zm, name(i1), name(i2),
     &                     name(jcomp), ecc, semi, pcrit, gamma(ipair)
    8         FORMAT (/,' NEW QUINTUP2    T MT NM1 NMG NM3 E0 A0 PC G0',
     &                                    f10.2,f6.2,3i6,f6.2,1p,3e10.2)
          END IF
      END IF
*
*       Include diagnostics for double triple as [[B,S],[B,S]].
      IF (name(jcomp).LT.0) THEN
          jpair = jcomp - n
          j1 = 2*jpair - 1
          CALL findj(i1,ji,im)
          CALL findj(j1,jj,jm)
          IF (kz(19).GE.3.AND.kz(27).LT.3) THEN
              jlist(1) = i1
              jlist(2) = i2
              jlist(3) = j1
              jlist(4) = j1 + 1
              jlist(5) = ji
              jlist(6) = jj
              CALL newtev(6,ix)
          END IF
          ai = -0.5*(cm(1,im) + cm(2,im))/hm(im)
          aj = -0.5*(cm(1,jm) + cm(2,jm))/hm(jm)
          zm = (body(i) + body(jcomp))*smu
          gx = max(gamma(ipair),gamma(jpair))
          WRITE (6,10)  time+toff, zm, name(i1), name(i2), name(j1),
     &                  name(2*jpair), ecc, ai, aj, r(ipair), r(jpair),
     &                  pcrit, gx
   10     FORMAT (/,' NEW HITRIP    T MT NM E0 AI AJ RI RJ PC GX ',
     &                              f9.2,f6.2,4i6,f6.2,1p,6e10.2)
      END IF
*
*       Ensure correct coordinates & velocities.
      CALL resolv(ipair,3)
*
*       Check optional diagnostics for hierarchy.
      IF ((kz(18).EQ.1.OR.kz(18).EQ.3).AND.kstar(i).LE.20) THEN
          CALL hiarch(ipair)
      END IF
*
*       Increase merger counter and set current merger index.
      nmerg = nmerg + 1
      nmerge = nmerge + 1
      imerge = nmerge
*
*       Save component masses and evaluate reduced mass of inner binary.
      cm(1,imerge) = body(2*ipair-1)
      cm(2,imerge) = body(2*ipair)
      cm(3,imerge) = 0.0d0
      zmu = body(2*ipair-1)*body(2*ipair)/body(n+ipair)
*
*       Set current energy of any second perturbed binary (from XVPRED).
      IF (jcomp.GT.n) THEN 
          CALL xvpred(jcomp,-1)
          jpair = jcomp - n
          IF (list(1,2*jpair-1).GT.0) h(jpair) = ht
*       Ensure that outer binary components are resolved.
          CALL ksres(jpair,j1,j2,0.0d0)
*       Enforce non-zero membership for potential correction in NBPOT.
          list(1,2*jpair-1) = 1
      ELSE
          CALL xvpred(jcomp,-1)
      END IF
*
*       Save the neighbours for correction procedure.
      nnb = list(1,i)
      DO 20 l = 1,nnb
          j = list(l+1,i)
          jpert(l) = j
   20 CONTINUE
*
*       Retain basic KS variables for explicit restart at merge termination.
      hm(imerge) = h(ipair)
      DO 25 k = 1,4
          um(k,imerge) = u(k,ipair)
          umdot(k,imerge) = udot(k,ipair)
   25 CONTINUE
*
*       Save stellar type.
      kstarm(imerge) = kstar(i)
*
*       Predict outer component to highest order if not on the block.
      IF (time - t0(jcomp1).GT.0.0d0) THEN
          CALL xvpred(jcomp1,-1)
      END IF
*
*       Set temporary KS components for hierarchical binary.
      icomp = 2*ipair - 1
      jcomp = icomp + 1
*
*       Obtain potential energy with respect to inner components.
      jlist(1) = icomp
      jlist(2) = jcomp
      CALL nbpot(2,nnb,pot1)
*
*       Save relative configuration.
      DO 30 k = 1,3
          xrel(k,imerge) = x(k,icomp) - x(k,jcomp)
          vrel(k,imerge) = x0dot(k,icomp) - x0dot(k,jcomp)
*       Initialize primary velocity of JCOMP1 (needed in KSIN2).
          x0dot(k,jcomp1) = xdot(k,jcomp1)
   30 CONTINUE
*
*       Include interaction of inner c.m. & neighbours before making new KS.
      icomp = icomp1
      jlist(1) = icomp
      CALL nbpot(1,nnb,pot2)
*
*       Copy mass of intruder to second KS component.
      body(jcomp) = body(jcomp1)
*
*       Replace name of #JCOMP with c.m. name (temporary for diagnostics).
      name2 = name(jcomp)
      name(jcomp) = name(jcomp1)
*
*       Form new c.m. and initialize KS variables (JPERT safe first call!).
      jcomp = jcomp1
      CALL ksin2(1)
*
*       See whether modifications due to second binary are needed.
      pot3 = 0.0d0
      pot4 = 0.0d0
      IF (jcomp1.LE.n) go to 50
*
*       Initialize unperturbed ghost binary of outer component.
      t0(2*jpair-1) = 1.0e+06
      list(1,2*jpair-1) = 0
*
*       Apply tidal correction for outer binary perturbers.
      jlist(1) = 2*jpair - 1
      jlist(2) = 2*jpair
      CALL nbpot(2,nnb,pot3)
      jlist(1) = jcomp1
      CALL nbpot(1,nnb,pot4)
*
*       Update the merger energy to maintain conservation.
      eb1 = body(2*jpair-1)*body(2*jpair)*h(jpair)/body(jcomp1)
      emerge = emerge + eb1
*
*       Save component masses and initialize ghost components.
      cm(3,imerge) = body(2*jpair-1)
      cm(4,imerge) = body(2*jpair)
      body(2*jpair-1) = 0.0d0
      body(2*jpair) = 0.0d0
      list(1,jcomp) = 0
*
*       Remove ghost from all neighbour lists.
   50 jlist(1) = jcomp1
      icm = n + kspair
      CALL nbrem(icm,1,nnb)
*
*       Specify JCOMP1 as ghost of zero mass.
      body(jcomp1) = 0.0d0
*
*       Initialize integration variables to prevent spurious predictions.
      DO 60 k = 1,3
          x0dot(k,jcomp1) = 0.0d0
          xdot(k,jcomp1) = 0.0d0
          f(k,jcomp1) = 0.0d0
          fdot(k,jcomp1) = 0.0d0
          d2(k,jcomp1) = 0.0d0
          d3(k,jcomp1) = 0.0d0
          d2r(k,jcomp1) = 0.0d0
          d3r(k,jcomp1) = 0.0d0
   60 CONTINUE
*
*       Set large value of T0 which avoids integration of ghost particle.
      t0(jcomp1) = 1.0e+06
*       Set large X0 & X to avoid perturber selection (no escape removal).
      x0(1,jcomp1) = 1.0e+06
      x(1,jcomp1) = 1.0e+06
*
*       Initialize c.m. & KS polynomials.
      icomp = 2*ipair - 1
      jcomp = icomp + 1
      CALL ksin2(2)
*
*       Define large negative c.m. name for identification & termination.
      name(icm) = name(icm) - 3*nzero
*
*       Set c.m. & ghost names for merger identification (include escape).
      namem(imerge) = name(icm)
      nameg(imerge) = name(jcomp1)
      name(jcomp) = name2
*
*       Copy stability limit for termination test A(1 - E) < R0 in KSINT.
      r0(ipair) = pcrit 
*
*       Update merger energy to maintain conservation.
      dphi = (pot2 - pot1) + (pot4 - pot3)
      emerge = emerge + zmu*hm(imerge) + dphi
*
*       Set IPHASE < 0 for new time-step list in routine INTGRT.
      iphase = -1
*
      RETURN
*
      END