reduce.f

Go to the documentation of this file.

      SUBROUTINE reduce(IESC,JESC,ISUB)
*
*
*       Reduction of chain.
*       -------------------
*
      include 'common6.h'
      parameter(nmx=10,nmx2=2*nmx,nmx3=3*nmx,nmx4=4*nmx,
     &            nmx8=8*nmx,nmxm=nmx*(nmx-1)/2)
      REAL*8  m,mass,mc,mij,mkk,xcm(3),vcm(3),dxc(3),dvc(3),cg(6)
      common/chain1/  xch(nmx3),vch(nmx3),m(nmx),
     &                zz(nmx3),wc(nmx3),mc(nmx),
     &                xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
     &                mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
      common/chainc/  xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
     &                listc(lmax)
      common/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      common/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
      common/ccoll2/  qk(nmx4),pk(nmx4),rik(nmx,nmx),SIZE(nmx),vstar1,
     &                ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
      common/incond/  x4(3,nmx),xdot4(3,nmx)
      REAL*8  firr(3),fd(3)
*
*
*       Define discrete time for new polynomials (DT2 < 0 is OK).
      time0 = time
      dt2 = t0s(isub) + timec - tprev
      dt8 = (tblock - tprev)/8.0d0
      time = tprev + nint(dt2/dt8)*dt8
      time = min(tblock,time)
*       Re-define initial epoch for consistency (ignore phase error).
      t0s(isub) = time - timec
*
      IF (kz(30).GT.2) THEN
          WRITE (6,1)  time0+toff, time+toff, dt2, dt8
    1     FORMAT (' REDUCE:   TIME0 TIME DT2 DT8 ',2f12.6,1p,2e10.2)
      END IF
*
*       Make new chain from remaining members.
    2 lk = 0
      DO 10 l = 1,nch
          IF (l.EQ.iesc) go to 10
          DO 5 k = 1,3
              lk = lk + 1
              xch(lk) = x4(k,l)
              vch(lk) = xdot4(k,l)
    5     CONTINUE
   10 CONTINUE
*
*       Reduce chain membership and mass (global & local COMMON).
      nch = nch - 1
      nn = nch
      mass = mass - m(iesc)
*
*       Improve coordinates & velocities of c.m. body to order F3DOT.
      CALL xvpred(ich,-1)
*
*       Predict X & XDOT of pertubers to order FDOT.
      nb1 = listc(1) + 1
      DO 15 l = 2,nb1
          j = listc(l)
          CALL xvpred(j,-2)
   15 CONTINUE
*
*       Set new c.m. for reduced system and save old c.m. variables.
      DO 20 k = 1,3
          dxc(k) = -bodyc(iesc)*x4(k,iesc)/(body(ich) - bodyc(iesc))
          dvc(k) = -bodyc(iesc)*xdot4(k,iesc)/(body(ich) - bodyc(iesc))
          xcm(k) = x(k,ich) + dxc(k)
          vcm(k) = xdot(k,ich) + dvc(k)
          cm(k) = x(k,ich)
          cm(k+3) = xdot(k,ich)
   20 CONTINUE
*
*       Re-define new chain variables w.r. to modified c.m. (NB! retain X4).
      lk = 0
      DO 30 l = 1,nch
          DO 25 k = 1,3
              lk = lk + 1
              xch(lk) = xch(lk) - dxc(k)
              vch(lk) = vch(lk) - dvc(k)
   25     CONTINUE
   30 CONTINUE
*
*       Save original mass & neighbour radius of c.m. body.
      bodych = body(ich)
      rs0 = rs(ich)
*
*       Search for global index of escaper.
      DO 40 j = ifirst,ntot
          IF (name(j).EQ.namec(iesc)) THEN
              i = j
              IF (body(j).GT.0.0d0) WRITE (6,35)  i, iesc, namec(iesc)
   35         FORMAT (' WARNING!   NON-ZERO GHOST    I IESC NAMEC ',3i5)
              go to 55
          END IF
   40 CONTINUE
*
*       Switch to another reference body since #ICH is escaping (NAME = 0).
      i = ich
      IF (iesc.GT.1) THEN
          new = 1
      ELSE
          new = 2
      END IF
*
*       Identify global index of new reference body (including binary c.m.).
      DO 45 j = ifirst,ntot
          IF (name(j).EQ.namec(new)) THEN
              ich = j
              go to 50
          END IF
   45 CONTINUE
*
*       Include warning if no reference body (this should not occur).
      WRITE (6,48)  iesc, namec(new)
   48 FORMAT (' REDUCE:   DANGER!   NO REFERENCE BODY    IESC NAME',2i5)
      nch = nch + 1
      nn = nch
      mass = mass + m(iesc)
      go to 100
*
*       Copy neighbour list of old c.m. body for routine NBREST.
   50 nnb = list(1,i)
      DO 52 l = 2,nnb+1
          jpert(l-1) = list(l,i)
   52 CONTINUE
*
*       Restore ghost to lists of neighbours (body #I will be skipped).
      jlist(1) = ich
      CALL nbrest(i,1,nnb)
*
      IF (kz(30).GT.1) THEN
          WRITE (6,53)  name0, name(ich), ich
   53     FORMAT (' REDUCE:    SWITCH C.M.    NAME0 NAMECH ICH ',3i5)
      END IF
*
*       Exchange name of reference body and initialize new c.m. name.
      name(i) = name0
      name0 = name(ich)
      name(ich) = 0
*
*       Update total mass and initialize new c.m. body variables.
   55 body(ich) = bodych - bodyc(iesc)
      cm(7) = body(ich)
      t0(ich) = time
      DO 60 k = 1,3
          x(k,ich) = xcm(k)
          x0(k,ich) = xcm(k)
          xdot(k,ich) = vcm(k)
          x0dot(k,ich) = vcm(k)
   60 CONTINUE
*
*       Restore the mass and transform to global coordinates & velocities.
      body(i) = bodyc(iesc)
      t0(i) = time
      DO 65 k = 1,3
          x(k,i) = x4(k,iesc) + cm(k)
          xdot(k,i) = xdot4(k,iesc) + cm(k+3)
          x0(k,i) = x(k,i)
          x0dot(k,i) = xdot(k,i)
   65 CONTINUE
*
*       Remove chain (and clump) mass & reference name of escaper.
      DO 70 l = iesc,nch
          m(l) = m(l+1)
          bodyc(l) = bodyc(l+1)
          namec(l) = namec(l+1)
          SIZE(l) = SIZE(l+1)
          istar(l) = istar(l+1)
          bodys(l,isub) = bodys(l+1,isub)
          names(l,isub) = names(l+1,isub)
   70 CONTINUE
*
*       Ensure current coordinates & velocities for chain components.
      CALL xcpred(1)
*
*       Initialize neighbour list using current radius of (old) #ICH.
      CALL nblist(ich,rs0)
*
*       Set dominant F & FDOT on body #I for #ICH in FPOLY2.
      jlist(1) = ich
      jlist(2) = jclose
      np = 2
      IF (jclose.EQ.0) np = 1
      CALL fclose(i,np)
*
*       Perform re-initialization of c.m. polynomials & perturber list.
      CALL reinit(isub)
*
*       Copy new chain coordinates & velocities to standard variables.
      lk = 0
      DO 80 l = 1,nch
          DO 75 k = 1,3
              lk = lk + 1
              x4(k,l) = xch(lk)
              xdot4(k,l) = vch(lk)
   75     CONTINUE
   80 CONTINUE
*
*       Copy neighbour list of cm. body for routine NBREST.
      nnb = list(1,ich)
      DO 82 l = 2,nnb+1
          jpert(l-1) = list(l,ich)
   82 CONTINUE
*
*       Restore ghost to lists of neighbours (body #ICH will be skipped).
      jlist(1) = i
      CALL nbrest(ich,1,nnb)
*
*       Make new neighbour list.
      CALL nblist(i,rs0)
*
*       Distinguish between single particle and binary (JESC = 0 & > 0).
      IF (jesc.EQ.0) THEN
*       Initialize force polynomials & time-steps (add differential F & FD).
          iphase = 8
          CALL fpoly1(i,i,0)
          DO 85 k = 1,3
              firr(k) = 0.0
              fd(k) = 0.0
   85     CONTINUE
          CALL fchain(i,0,x(1,i),xdot(1,i),firr,fd)
          rij2 = 0.0
          vij2 = 0.0
          rdot = 0.0
          DO 86 k = 1,3
              f(k,i) = f(k,i) + firr(k)
              fdot(k,i) = fdot(k,i) + fd(k)
              rij2 = rij2 + (x(k,i) - x(k,ich))**2
              vij2 = vij2 + (xdot(k,i) - xdot(k,ich))**2
              rdot = rdot + (x(k,i) - x(k,ich))*(xdot(k,i)-xdot(k,ich))
   86     CONTINUE
          CALL fpoly2(i,i,0)
*       Reduce regular steps for large velocities to minimize NNB change.
          j = ich
          DO 87 l = 1,2
              vi2 = xdot(1,j)**2 + xdot(2,j)**2 + xdot(3,j)**2
              IF (vi2.GT.2.0.AND.stepr(j).GE.8.0*step(j)) THEN
                  stepr(j) = stepr(j)/8.0d0
               END IF
               j = i
   87     CONTINUE
*       Check high-velocity ejection for outwards hyperbolic motion.
          hi = 0.5*vij2 - (body(i) + body(ich))/sqrt(rij2)
          IF (hi.GT.0.0.AND.rdot.GT.0.0) THEN
               CALL hivel(i)
          END IF
          iphase = -1
*       Include case of escaping binary (set JESC < 0 for new KS).
      ELSE IF (jesc.GT.0) THEN
          iclose = i
          IF (jesc.GT.iesc) jesc = jesc - 1
          iesc = jesc
          jesc = -1
          go to 2
      ELSE
*       Initialize KS regularization after second reduction.
          icomp = min(iclose,i)
          jcomp = max(iclose,i)
          CALL ksreg
      END IF
*
*       Re-activate any dormant binary.
      IF (i.GT.n.AND.jesc.EQ.0) THEN
          CALL renew(i)
      END IF
*
*       Prepare c.m. check.
      DO 88 k = 1,6
          cg(k) = 0.0
   88 CONTINUE
*
      lk = 0
      DO 95 l = 1,nch
          DO 90 k = 1,3
              lk = lk + 1
              cg(k) = cg(k) + bodyc(l)*xch(lk)
              cg(k+3) = cg(k+3) + bodyc(l)*vch(lk)
   90     CONTINUE
   95 CONTINUE
*
      DO 96 k = 1,6
          cg(k) = cg(k)/body(ich)
   96 CONTINUE
*
      IF (kz(30).GT.2) THEN
          WRITE (6,97)  time+toff, (cg(k),k=1,6)
   97     FORMAT (' REDUCE:   T CG ',f12.5,1p,6e9.1)
      END IF
*
  100 RETURN
*
      END