reset2.f

Go to the documentation of this file.

      SUBROUTINE reset2
*
*
*      Termination 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
      REAL*8  xx(3,3),vv(3,3)
      SAVE istab
      DATA istab /0/
*
*
*       Set index of disrupted pair and save output parameters.
      ipair = kspair
      i = n + ipair
      e1 = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(i)
      g1 = gamma(ipair)
      r1 = r(ipair)
      semi1 = -0.5*body(i)/h(ipair)
      ecc2 = (1.0 - r1/semi1)**2 + tdot2(ipair)**2/(body(i)*semi1)
      ecc = sqrt(ecc2)
*
*       Locate current position in the merger table.
      imerge = 0
      DO 1 k = 1,nmerge
          IF (namem(k).EQ.name(i)) imerge = k
    1 CONTINUE
*
*       Produce info on quintuplet or sextuplet (NAMEG has outer c.m. name).
      IF (nameg(imerge).GT.nzero) THEN
          ri2 = 0.0
          DO 2 k = 1,3
              ri2 = ri2 + (x(k,i) - rdens(k))**2
    2     CONTINUE
          ri = sqrt(ri2)
          pmin = semi1*(1.0 - ecc)
          eb = 0.5*body(2*ipair-1)*body(2*ipair)/semi1
          eb = eb*float(n-npairs)/zkin
          zm = body(i)*smu
          which1 = ' QUINTUPLET'
*       Determine previous merger index with NAME greater by 2*NZERO.
          jm = 1
          DO 3 k = 1,nmerge
              IF (namem(k).EQ.name(i) + 2*nzero) jm = k
    3     CONTINUE
*       Use previous merger index to identify ghost from earlier level.
          DO 4 j = ifirst,ntot
              IF (name(j).EQ.nameg(jm)) jg = j
    4     CONTINUE
          IF (jg.GT.n) which1 = ' SEXTUPLET '
          WRITE (6,5)  which1, time+toff, zm, name(2*ipair-1),
     &                 name(2*ipair), nameg(imerge), ri, ecc, eb, semi1,
     &                 pmin, g1
    5     FORMAT (' END',a11,'   T MT NM1 NM2 NM3 RI E1 EB1 A1 PM G1 ',
     &                           f9.2,f6.2,3i6,2f6.2,f6.1,1p,3e10.2)
      END IF
*
*       Check presence of [[B,S],S] quartet or [[B,B],S] quintuplet.
      IF (nameg(imerge).GT.0.AND.nameg(imerge).LE.nzero) THEN
          i1 = 2*ipair - 1
          i2 = 2*ipair
          CALL findj(i1,jg,im)
          zm = body(i)*smu
          IF (name(i2).LE.nzero) THEN
              WRITE (6,9)  time+toff, zm, name(i1), name(i2), name(jg),
     &                     ecc, semi1, r1, g1
    9         FORMAT (/,' END QUARTET    T MT NM1 NMG NM3 E1 A1 R1 G1 ',
     &                                   f9.2,f6.2,3i6,f6.2,1p,3e10.2)
          ELSE
              WRITE (6,11)  time+toff, zm, name(i1), name(i2), name(jg),
     &                      ecc, semi1, r1, g1
   11         FORMAT (/,' END QUINTUP2    T MT NM1 NMG NM3 E1 A1 R1 G1',
     &                                    f10.2,f6.2,3i6,f6.2,1p,3e10.2)
          END IF
      END IF
*
*       Include diagnostics for double triple ([[B,S],[B,S]]).
      IF (nameg(imerge).LT.0) THEN
          i1 = 2*ipair - 1
          j1 = 2*ipair
*       Obtain current merger index (note NAMEM(JM) = NAMEG(JG)).
          CALL findj(j1,jj,jg)
          im = 1
          jm = 1
*       Determine original merger indices of the two inner binaries.
          DO 16 k = 1,nmerge
              IF (namem(k).EQ.name(i) + 2*nzero) im = k
              IF (namem(k).EQ.nameg(jg)) jm = k
   16     CONTINUE
          ai = -0.5*(cm(1,im) + cm(2,im))/hm(im)
          aj = -0.5*(cm(1,jm) + cm(2,jm))/hm(jm)
          pmin = semi1*(1.0 - ecc)
          zm = body(i)*smu
          WRITE (6,18)  time+toff, zm, name(i1), nameg(im), -namem(jm),
     &                  nameg(jm), ecc, ai, aj, r(ipair), semi1, pmin,
     &                  g1
   18     FORMAT (/,' END HITRIP    T MT NM E1 AI AJ R1 A1 PM G1 ',
     &                              f9.2,f6.2,4i6,f6.2,1p,6e10.2)
      END IF
*
*       Check optional diagnostics for hierarchy.
      IF ((kz(18).EQ.1.OR.kz(18).EQ.3).AND.kstar(imerge).LE.20) THEN
          CALL hiarch(ipair)
      END IF
*
*       Save neighbours for correction procedure.
      nnb = list(1,i) + 1
      DO 6 l = 2,nnb
          j = list(l,i)
          jpert(l) = j
    6 CONTINUE
*
*       Ensure that c.m. coordinates are known to highest order.
      CALL xvpred(i,0)
*
*       Predict neighbour coordinates & velocities (XDOT used by FPOLY1).
      DO 7 l = 2,nnb
          j = jpert(l)
          CALL xvpred(j,0)
    7 CONTINUE
*
*       Obtain current coordinates & velocities and specify KS components.
      CALL resolv(ipair,2)
      icomp = 2*ipair - 1
      jcomp = icomp + 1
*
*       Initialize mass, coordinates and velocities for new c.m. body.
      body(i) = body(icomp)
      DO 8 k = 1,3
          x(k,i) = x(k,icomp)
          xdot(k,i) = xdot(k,icomp)
          x0dot(k,i) = xdot(k,icomp)
    8 CONTINUE
*
*       Add outer component to perturber list.
      jpert(1) = jcomp
*
*       Sum first part of potential energy correction due to tidal effect.
      jlist(1) = icomp
      CALL nbpot(1,nnb,pot1)
*
*       Find correct location of ghost particle using identification name.
      icm = i
      jcomp1 = jcomp
      DO 10 i = 1,ntot
          IF (body(i).EQ.0.0d0.AND.name(i).EQ.nameg(imerge)) jcomp1 = i
   10 CONTINUE
*
*       Regularize two-body configuration if JCOMP1 cannot be identified.
      IF (jcomp.EQ.jcomp1) THEN
          WRITE (6,12)  imerge, nameg(imerge), jcomp
   12     FORMAT (/,5x,'WARNING!    RESET2    JCOMP NOT IDENTIFIED ',
     &                       '   IM =',i3,'  NAMEG =',i6,'  JCOMP =',i6)
          go to 100
      END IF
*
*       Initialize basic variables for ghost and new c.m. (JCOMP -> JCOMP1).
      j1 = jcomp1
      j = jcomp
   13 t0(j1) = time
      body(j1) = body(j)
      DO 14 k = 1,3
          x(k,j1) = x(k,j)
          x0(k,j1) = x(k,j)
          xdot(k,j1) = xdot(k,j)
          x0dot(k,j1) = xdot(k,j)
   14 CONTINUE
      IF (j.EQ.jcomp) THEN
          j1 = icm
          j = icomp
          go to 13
      END IF
*
*       Restore masses, coordinates & velocities of inner binary.
      body(icomp) = cm(1,imerge)
      body(jcomp) = cm(2,imerge)
      zm = -body(icomp)/(body(icomp) + body(jcomp))
*
*       Begin with second component since ICOMP holds new c.m. variables.
      i = jcomp
      DO 20 kcomp = 1,2
          DO 15 k = 1,3
              x(k,i) = x(k,icomp) + zm*xrel(k,imerge)
              x0dot(k,i) = x0dot(k,icomp) + zm*vrel(k,imerge)
              xdot(k,i) = x0dot(k,i)
*       Note that XDOT is only needed for improved polynomials of JCOMP.
   15     CONTINUE
          i = icomp
          zm = body(jcomp)/(body(icomp) + body(jcomp))
   20 CONTINUE
*
*       Copy KS variables for inner binary (small TDOT2 near apo/peri).
      i1 = 2*ipair - 1
      t0(i1) = time
      list(1,i1) = 1
      h(ipair) = hm(imerge)
      r(ipair) = 0.0d0
      tdot2(ipair) = 0.0d0
      DO 30 k = 1,4
          u(k,ipair) = um(k,imerge)
          u0(k,ipair) = u(k,ipair)
          udot(k,ipair) = umdot(k,imerge)
          r(ipair) = r(ipair) + u(k,ipair)**2
   30 CONTINUE
*
*       Add #JCOMP1 to all neighbour lists containing ICM.
      jlist(1) = jcomp1
      CALL nbrest(icm,1,nnb)
*
*       Add #ICM to neighbour list of #JCOMP1.
      jlist(1) = jcomp1
      jpert(1) = icm
      CALL nbrest(icm,1,1)
*
*       Form new neighbour list for #ICM (NB! NBREST(JCOMP1,1,1) skips).
      rs0 = rs(icm)
      CALL nblist(icm,rs0)
*
*       Get new list for #JCOMP1 (NNB = 1 from MERGE2/NBREM for ghost c.m.).
      rs0 = rs(jcomp1)
      CALL nblist(jcomp1,rs0)
*
*       Initialize force polynomial for outer component using resolved c.m.
      CALL fpoly1(jcomp1,jcomp1,0)
      CALL fpoly2(jcomp1,jcomp1,0)
*
*       Initialize c.m. polynomials and activate inner binary.
      CALL ksin2(3)
*
*       Restore original name of inner hierarchy (c.m. NAME set in MERGE2).
      name(icm) = name(icm) + 2*nzero
*
*       Locate position of inner binary in merger table (IM = 1 for safety).
      im = 1
      DO 35 k = 1,nmerge
          IF (namem(k).EQ.name(icm)) im = k
   35 CONTINUE
*
*       Determine inclination between inner relative motion and outer orbit.
      rv = 0.0
      DO 40 k = 1,3
          xx(k,1) = xrel(k,im)
          xx(k,2) = 0.0
          xx(k,3) = x(k,jcomp)
          vv(k,1) = vrel(k,im)
          vv(k,2) = 0.0
          vv(k,3) = xdot(k,jcomp)
          rv = rv + xrel(k,im)*vrel(k,im)
   40 CONTINUE
      CALL inclin(xx,vv,x(1,icm),xdot(1,icm),angle)
*
*       Evaluate stability parameter from current elements (minimum ECC).
      semi0 = -0.5*body(icomp)/hm(im)
      semi = -0.5*body(icm)/h(ipair)
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(icm)*semi)
      ecc1 = sqrt(ecc2)
*       Form inner eccentricity from two-body elements.
      rin = sqrt(xrel(1,im)**2 + xrel(2,im)**2 + xrel(3,im)**2)
      ecc2 = (1.0 - rin/semi0)**2 + rv**2/(body(icomp)*semi0)
      ecc = sqrt(ecc2)
*
*       Evaluate the general stability function.
      IF (ecc1.LT.1.0) THEN
          nst = nstab(semi0,semi,ecc,ecc1,angle,cm(1,imerge),
     &                                    cm(2,imerge),body(jcomp))
          IF (nst.EQ.0) THEN
              pcrit = 0.98*semi*(1.0 - ecc1)
              pcr = stability(cm(1,imerge),cm(2,imerge),body(jcomp),
     &                                          ecc,ecc1,angle)*semi0
              istab = istab + 1
              IF (istab.LT.50) THEN
                  WRITE (6,41)  ecc, ecc1, semi0, semi, pcrit, pcr
   41             FORMAT (' STABLE2    E E1 A A1 PCR PC99 ',
     &                                 2f7.3,1p,4e10.2)
              END IF
          ELSE
              pcrit = 1.02*semi*(1.0 - ecc1)
          END IF
      ELSE
          pcrit = stability(cm(1,imerge),cm(2,imerge),body(jcomp),
     &                                      ecc,ecc1,angle)*semi0
      END IF
*
*       Set critical pericentre distance for stability check.
      r0(ipair) = pcrit
*
*       Form new force polynomials for perturbers inside 10*R1 (skip J > N).
      nnb1 = list(1,icm) + 1
      rs2 = (10.0*r1)**2
      DO 46 l = 2,nnb1
          j = list(l,icm)
          rij2 = 0.0
          DO 42 k = 1,3
              rij2 = rij2 + (x(k,icm) - x(k,j))**2
   42     CONTINUE
          IF (rij2.LT.rs2.AND.j.NE.jcomp1.AND.j.LE.n) THEN
              t0(j) = time
              DO 45 k = 1,3
                  x0dot(k,j) = xdot(k,j)
   45         CONTINUE
              CALL fpoly1(j,j,0)
              CALL fpoly2(j,j,0)
          END IF
   46 CONTINUE
*
*       Rename perturber list for routine NBPOT.
      jpert(1) = jcomp1
*
*       Restore Roche stage indicator (any ghost c.m. is OK).
      kstar(icm) = kstarm(imerge)
*
*       See whether the outer component is a single or composite particle.
      pot3 = 0.0d0
      pot4 = 0.0d0
      IF (jcomp1.LE.n) go to 50
*
*       Restore component masses for outer binary.
      jpair = jcomp1 - n
      body(2*jpair-1) = cm(3,imerge)
      body(2*jpair) = cm(4,imerge)
*
*       Update look-up times & radii and check possible Roche condition.
      IF (kz(19).GE.3) THEN
          IF (kstar(jcomp1).GT.0.AND.kstar(jcomp1).LE.10) THEN
              CALL trflow(jpair,dtr)
              tev(jcomp1) = time + dtr
              tmdot = min(tev(jcomp1),tmdot)
              tmdot = min(tev(2*jpair),tmdot)
          END IF
      END IF
*
*       Obtain coordinates & velocities of unperturbed binary components.
      CALL resolv(jpair,1)
*
*       Select new perturbers and initialize polynomials for KS motion.
      CALL kslist(jpair)
      CALL kspoly(jpair,1)
*
*       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.
      eb2 = body(2*jpair-1)*body(2*jpair)*h(jpair)/body(jcomp1)
      emerge = emerge - eb2
*
      e2 = e1/eb2
      eb2 = eb2/be(3)
      dp = pot4 - pot3
      IF (kz(15).GT.1) THEN
          WRITE (6,48)  jpair, h(jpair), body(2*jpair-1),
     &                  body(2*jpair), e2, eb2, r(jpair), gamma(jpair),
     &                  dp
   48     FORMAT (' END OUTER MERGER',i5,'  H =',f7.1,'  M =',2f7.4,
     &            '  E1 =',f6.3,'  EB2 =',f6.3,'  RB2 =',1pe8.1,
     &            '  G2 =',e8.1,'  DP =',e8.1)
      END IF
*
*       Include interaction of body #ICOMP & JCOMP with perturbers.
   50 jlist(1) = icomp
      jlist(2) = jcomp
      CALL nbpot(2,nnb,pot2)
*
*       Form square of c.m. velocity correction due to tidal effects.
*     VI2 = X0DOT(1,ICM)**2 + X0DOT(2,ICM)**2 + X0DOT(3,ICM)**2
      dphi = (pot2 - pot1) + (pot4 - pot3)
*     CORR = 1.0 + 2.0*DPHI/(BODY(ICM)*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,ICM) = SQRT(CORR)*X0DOT(K,ICM)
*  60 CONTINUE
*
*       Modify the merger energy to maintain conservation.
      eb = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(icm)
      emerge = emerge - eb + dphi
*
      e1 = e1/eb
      eb = eb/be(3)
      IF (kz(15).GT.1) THEN
          WRITE (6,65)  imerge, time+toff, body(2*ipair-1),
     &                  body(2*ipair), r1, semi1, eb, e1,
     &                  gamma(ipair), g1, nnb-1
   65     FORMAT (' END MERGE2',i3,'  T =',f8.2,'  M =',2f7.4,
     &            '  R1 =',1pe8.1,'  A1 =',e8.1,'  EB =',0pf6.3,
     &            '  E1 =',f6.3,'  GB =',1pe8.1,'  G =',0pf6.3,
     &            '  NB =',i3)
      END IF
*
*       Check Roche look-up time.
      IF (kstar(icm).GT.0.AND.kstar(icm).LE.10) THEN
          k = icm - n
          CALL trflow(k,dtr)
          tev(icm) = min(tev(icm),time + dtr)
          tmdot = min(tev(icm),tmdot)
      END IF
*
*       Reduce merger counter and update tables (unless last or only pair).
   70 nmerge = nmerge - 1
      DO 80 l = imerge,nmerge
          l1 = l + 1
          hm(l) = hm(l1)
          nameg(l) = nameg(l1)
          namem(l) = namem(l1)
          kstarm(l) = kstarm(l1)
          DO 74 k = 1,3
              xrel(k,l) = xrel(k,l1)
              vrel(k,l) = vrel(k,l1)
   74     CONTINUE
          DO 75 k = 1,4
              cm(k,l) = cm(k,l1)
              um(k,l) = um(k,l1)
              umdot(k,l) = umdot(k,l1)
   75     CONTINUE
   80 CONTINUE
*
*       Examine merger list for possible escapers (retain up to 3 levels).
      DO 90 l = 1,nmerge
          DO 85 j = 1,npairs
              IF (namem(l).EQ.name(n+j).OR.
     &            namem(l).EQ.name(n+j) + 2*nzero.OR.
     &            namem(l).EQ.name(n+j) + 4*nzero) go to 90
   85     CONTINUE
*       Remove tables for any merger not identified.
          imerge = l
          go to 70
   90 CONTINUE
*
*       Set IPHASE < 0 for new time-step list in INTGRT.
      iphase = -1
*
  100 RETURN
*
      END