reset.f

Go to the documentation of this file.

      SUBROUTINE reset
*
*
*       Restore hierarchical configuration.
*       -----------------------------------
*
      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)
*
*
*       Set index of disrupted pair and save output parameters.
      ipair = kspair
      i = n + ipair
*
*       Check special termination of double hierarchy (routine RESET2).
      IF (name(i).LT.-2*nzero) THEN
          CALL reset2
          go to 100
      END IF
*
      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)
*
*       Locate current position in the merger table.
      imerge = 0
      DO 1 k = 1,nmerge
          IF (namem(k).EQ.name(i)) imerge = k
    1 CONTINUE
*
*       Check optional diagnostics for hierarchy.
      IF ((kz(18).EQ.1.OR.kz(18).EQ.3).AND.kstarm(imerge).LE.20) THEN
          CALL hiarch(ipair)
      END IF
*
*       Save neighbours for correction procedure and rename if moved up.
      nnb = list(1,i) + 1
      DO 4 l = 2,nnb
          j = list(l,i)
          IF (j.GT.i) j = j - 1
          IF (j.LE.2*npairs.AND.j.GT.2*ipair) j = j - 2
          jsave(l) = j
    4 CONTINUE
*
*       Ensure that c.m. coordinates are known to highest order.
      CALL xvpred(i,0)
*
*       Restore original c.m. name and terminate outer KS pair.
      name(i) = nzero - name(i)
      CALL ksterm
*
*       Predict neighbour coordinates & velocities (XDOT used by FPOLY1).
      DO 5 l = 2,nnb
          jpert(l) = jsave(l)
          j = jpert(l)
          CALL xvpred(j,0)
    5 CONTINUE
*
*       Add outer component to perturber list and set old dipole range.
      jpert(1) = jcomp
*     RCRIT2 = CMSEP2*(XREL(1,IMERGE)**2 + XREL(2,IMERGE)**2 +
*    &                                                XREL(3,IMERGE)**2)
*
*       Sum first part of potential energy correction due to tidal effect.
      jlist(1) = icomp
      CALL nbpot(1,nnb,pot1)
*
*       Find the nearest neighbour and reduce steps of active perturbers.
      rjmin2 = 100.0
      jmin = n
      DO 8 l = 1,nnb
          j = jpert(l)
          rij2 = (x(1,j) - x(1,icomp))**2 + (x(2,j) - x(2,icomp))**2 +
     &                                      (x(3,j) - x(3,icomp))**2
*       Identify the closest perturber for dominant force calculation.
          IF (rij2.LT.rjmin2.AND.j.NE.jcomp) THEN
              rjmin2 = rij2
              jmin = j
          END IF
*         IF (RIJ2.LT.RCRIT2) THEN
*       Reduce step of inner binary perturbers (c.m. approximation used).
*             STEP(J) = MAX(0.5D0*STEP(J),TIME - T0(J))
*         END IF
    8 CONTINUE
*
*       Find correct location of ghost particle using identification name.
      jcomp1 = jcomp
*       Note that ghost may be saved in an old binary c.m. (search NTOT). 
      DO 10 i = 1,ntot
          IF (body(i).EQ.0.0d0.AND.name(i).EQ.nameg(imerge)) jcomp = i
   10 CONTINUE
*
*       Regularize two-body configuration if JCOMP cannot be identified.
      IF (jcomp.EQ.jcomp1) THEN
          WRITE (6,12)  nameg(imerge)
   12     FORMAT (/,5x,'DANGER!   JCOMP NOT IDENTIFIED IN RESET',
     &                                                  '   NAMEG =',i5)
          stop
      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
*
*       Add #JCOMP to neighbour lists containing ICOMP (KSREG sets c.m.).
      jlist(1) = jcomp
      CALL nbrest(icomp,1,nnb)
*
*       Ensure that any rare case of missing second component is included.
      DO 30 i = 1,ntot
          nnb1 = list(1,i) + 1
          DO 25 l = 2,nnb1
              IF (list(l,i).GT.icomp) go to 30
*
*       Now see whether #JCOMP has already been added.
              DO 22 k = l,nnb1
                  IF (list(k,i).EQ.jcomp) go to 30
   22         CONTINUE
*
*       Specify index #I and add body #JCOMP as above.
              jpert(1) = i
              CALL nbrest(icomp,1,1)
   25     CONTINUE
   30 CONTINUE
*
*       Set dominant F & FDOT on inner components including main intruder.
      jlist(1) = icomp
      jlist(2) = jcomp
      jlist(3) = jcomp1
      jlist(4) = jmin
*
      CALL fclose(icomp,4)
      CALL fclose(jcomp,4)
*
*       Initialize force polynomials for outer component using resolved c.m.
      CALL fpoly1(jcomp1,jcomp1,0)
      CALL fpoly2(jcomp1,jcomp1,0)
*
*       Rename perturber list for routine NBPOT.
      jpert(1) = jcomp
*
*       Copy basic KS variables for inner binary (small TDOT2 near apo/peri).
      jp1 = npairs + 1
      h(jp1) = hm(imerge)
      r(jp1) = 0.0d0
      tdot2(jp1) = 0.0d0
      DO 40 k = 1,4
          u(k,jp1) = um(k,imerge)
          u0(k,jp1) = u(k,jp1)
          udot(k,jp1) = umdot(k,imerge)
          r(jp1) = r(jp1) + u(k,jp1)**2
   40 CONTINUE
*
*       Save ghost index and re-activate inner binary (JCOMP <-> JCOMP1).
      jcomp1 = jcomp
      CALL ksreg
*
*       Restore Roche stage indicator (any ghost c.m. is OK).
      kstar(n+npairs) = 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
*       Reduce evolution time-scale after any delay during ghost stage.
          IF (tev(2*jpair-1).LT.time + 0.01*tcr) tev(2*jpair-1) = time
          IF (tev(2*jpair).LT.time + 0.01*tcr) tev(2*jpair) = time
          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
*       Note that inner binary correction is included with IPAIR.
      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,45)  jpair, h(jpair), body(2*jpair-1),
     &                  body(2*jpair), e2, eb2, r(jpair), gamma(jpair),
     &                  dp
   45     FORMAT (' END QUAD',i4,'  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,NTOT)**2 + X0DOT(2,NTOT)**2 + X0DOT(3,NTOT)**2
      dphi = (pot2 - pot1) + (pot4 - pot3)
*     CORR = 1.0 + 2.0*DPHI/(BODY(NTOT)*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,NTOT) = SQRT(CORR)*X0DOT(K,NTOT)
*  60 CONTINUE
*
*       Modify the merger energy to maintain conservation.
      eb = body(2*npairs-1)*body(2*npairs)*h(npairs)/body(ntot)
*       Note that EMERGE may contain escaped mergers.
      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*npairs-1),
     &                  body(2*npairs), r1, semi1, eb, e1,
     &                  gamma(npairs), g1, nnb
   65     FORMAT (' END MERGER',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(ntot).GE.10.AND.kstar(ntot).LE.20) THEN
*       Reduce evolution time-scale after any delay during ghost stage.
          IF (tev(2*npairs-1).LT.time + 0.01*tcr) tev(2*npairs-1) = time
          CALL trflow(npairs,dtr)
*       Note DTR < 0 is possible for active Roche after termination.
          tev(ntot) = time + max(dtr,0.0d0)
          tmdot = min(tev(ntot),tev(icomp),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
*
*       Include consistency check on merger names.
      DO 96 i = n+1,ntot
          IF (name(i).LT.0) THEN
              DO 94 l = 1,nmerge
                  IF (namem(l).EQ.name(i)) go to 96
   94         CONTINUE
              WRITE (6,95)  i, name(i), (namem(l),l=1,nmerge)
   95         FORMAT (' DANGER!    RESET    I NAMI NAMEM  ',10i6)
              stop
          END IF
   96 CONTINUE
*
*       Reduce merger energy on zero membership for consistency.
      IF (nmerge.EQ.0) THEN
          be(3) = be(3) - emerge
          emerge = 0.0
      END IF
*
*       Set phase indicator < 0 for new time-step list in routine INTGRT.
      iphase = -1
*
  100 RETURN
*
      END