merge.f

Go to the documentation of this file.

      SUBROUTINE merge
*
*
*       Merging of hierarchical binary.
*       -------------------------------
*
      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)
      INTEGER  jsave(lmax)
*
*
*       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 stable mergers (exclude flybys).
*       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.
*       ------------------------------------------------------------------
*
*
*       Check merger limit.
      IF (nmerge.GE.mmax) THEN
          WRITE (6,5)  nmerge
    5     FORMAT (/,5x,'WARNING!   MERGER LIMIT EXCEEDED   NMERGE =',i4)
          go to 100
      END IF
*
*       Use double hierarchy procedure if one binary is already a merger.
      IF (name(n+kspair).LT.0.OR.name(jcomp).LT.0) THEN
          CALL merge2
          go to 100
      END IF
*
*       Check switching of inner and outer binary for improved restart.
      IF (jcomp.GT.n) THEN
          semi2 = -0.5*body(jcomp)/h(jcomp-n)
          IF (semi2.GT.r(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
*
*       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 energy of the 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)
*     EB = ZMU*H(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
*       Reduce index of outer component if JPAIR moved up at termination.
          IF (ipair.LT.jpair) THEN
              jpair = jpair - 1
              jcomp1 = jcomp1 - 1
          END IF
      END IF
*
*       Save the neighbours for correction procedure & rename if moved up.
      nnb = list(1,i)
      DO 15 l = 1,nnb
          j = list(l+1,i)
          IF (j.GT.i) j = j - 1
          IF (j.LE.2*npairs.AND.j.GT.2*ipair) j = j - 2
          jsave(l) = j
   15 CONTINUE
*
*       Reduce steps of perturbers (new force uses c.m. approximation).
*     NP1 = LIST(1,2*IPAIR-1) + 1
*     DO 20 L = 2,NP1
*         J = LIST(L,2*IPAIR-1)
*         STEP(J) = MAX(0.5D0*STEP(J),TIME - T0(J))
*  20 CONTINUE
*
*       Delay saving KS variables in block version until end-point in KSTERM.
      IF (time.GT.tblock) THEN
*       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
      END IF
*
*       Save c.m. type index.
      kstarm(imerge) = kstar(i)
*
*       Terminate inner pair in order to merge components after updating.
      CALL ksterm
*
*       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
*
*       Copy perturber list (NB! JPERT used by KSTERM).
      DO 28 l = 1,nnb
          jpert(l) = jsave(l)
   28 CONTINUE
*
*       Obtain potential energy with respect to inner components.
      jlist(1) = icomp
      jlist(2) = jcomp
      CALL nbpot(2,nnb,pot1)
*
*       Save relative configuration and define old binary as composite body.
      DO 30 k = 1,3
          xrel(k,imerge) = x(k,icomp) - x(k,jcomp)
          vrel(k,imerge) = x0dot(k,icomp) - x0dot(k,jcomp)
          x(k,icomp) = (body(icomp)*x(k,icomp) + body(jcomp)*x(k,jcomp))
     &                                      /(body(icomp) + body(jcomp))
          x0dot(k,icomp) = (body(icomp)*x0dot(k,icomp) +
     &                      body(jcomp)*x0dot(k,jcomp))/
     &                                       (body(icomp) + body(jcomp))
*       Initialize primary velocity of JCOMP1 (needed in KSREG & KSINIT).
          x0dot(k,jcomp1) = xdot(k,jcomp1)
   30 CONTINUE
*
*       Form new c.m. body & associated ghost of zero mass and neighbours.
      body(icomp) = body(icomp) + body(jcomp)
      body(jcomp) = 0.0d0
      list(1,jcomp) = 0
*
*       Initialize integration variables to prevent spurious predictions.
      DO 40 k = 1,3
          x0dot(k,jcomp) = 0.0d0
          xdot(k,jcomp) = 0.0d0
          f(k,jcomp) = 0.0d0
          fdot(k,jcomp) = 0.0d0
          d2(k,jcomp) = 0.0d0
          d3(k,jcomp) = 0.0d0
          d2r(k,jcomp) = 0.0d0
          d3r(k,jcomp) = 0.0d0
   40 CONTINUE
*
*       Set large value of T0 which avoids integration of ghost particle.
      t0(jcomp) = 1.0e+06
*       Set large X0 & X to avoid neighbour problems (no escape removal).
      x0(1,jcomp) = 1.0e+06
      x(1,jcomp) = 1.0e+06
*
*       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
*       Note that inner binary correction is included with IPAIR.
      CALL nbpot(2,nnb,pot3)
      jlist(1) = jcomp1
      CALL nbpot(1,nnb,pot4)
*
*       Update 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
*
*       Include interaction of inner c.m. & neighbours to give net effect.
   50 jlist(1) = icomp
      CALL nbpot(1,nnb,pot2)
*
*       Form square of c.m. velocity correction due to tidal effects.
*     VI2 = X0DOT(1,ICOMP)**2 + X0DOT(2,ICOMP)**2 + X0DOT(3,ICOMP)**2
      dphi = (pot2 - pot1) + (pot4 - pot3)
*     CORR = 1.0 + 2.0*DPHI/(BODY(ICOMP)*VI2)
*     IF (CORR.LE.0.0D0) CORR = 0.0
*
*       Adjust c.m. velocity by net tidal energy correction.
*     DO 60 K = 1,3
*         X0DOT(K,ICOMP) = SQRT(CORR)*X0DOT(K,ICOMP)
*  60 CONTINUE
*
*       Remove ghost from local neighbour lists (any others done in REGINT).
      jlist(1) = jcomp
      CALL nbrem(icomp,1,nnb)
*
*       Also remove ghost from list of ICOMP (use NTOT as dummy here).
      jpert(1) = icomp
      CALL nbrem(ntot,1,1)
*
*       Check for sufficient perturbers (RMAX = apocentre set in IMPACT).
      IF (sqrt(cmsep2)*rmax.GT.2.0*rs(icomp).OR.nnb.LE.3) THEN
          fac = nnbmax/nnb
          fac = min(fac,4.0d0)
          rsi = fac**0.3333*rs(icomp)
          IF (rsi.GT.rs(i)) THEN
              CALL nblist(icomp,rsi)
          END IF
      END IF
*
*       Perform KS regularization of hierarchical system (ICOMP & JCOMP1).
      jcomp = jcomp1
      CALL ksreg
*
*       Define negative c.m. name for identification & termination.
      name(ntot) = nzero - name(ntot)
*
*       Set c.m. & ghost names for merger identification (include escape).
      namem(imerge) = name(ntot)
      nameg(imerge) = name(jcomp1)
*       Subtract 1 for hyperbolic flybys to count stable systems only.
      IF (h(npairs).GT.0.0) nmerg = nmerg - 1

*       Copy stability limit for termination test A(1 - E) < R0 in KSINT.
      r0(npairs) = pcrit 
*       Ensure that flag denotes new rather than primordial binary.
      list(2,jcomp) = 0
*
*       Add merger energy (or: in clude DPHI and activate X0DOT correction).
      emerge = emerge + zmu*hm(imerge) + dphi
*
*       Set phase indicator = -1 for new time-step list in routine INTGRT.
      iphase = -1
*
  100 RETURN
*
      END