absorb.f

Go to the documentation of this file.

      SUBROUTINE absorb(ISUB)
*
*
*       Absorption of chain member(s).
*       -----------------------------
*
      include 'common6.h'
      parameter(nmx=10,nmx3=3*nmx,nmx4=4*nmx,nmxm=nmx*(nmx-1)/2)
      REAL*8  m,mass,mc,mij,mkk,xcm(3),vcm(3)
      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)
*
*
*       Define discrete time for new polynomial (DT2 < 0 is OK).
      dt2 = t0s(isub) + timec - tprev
      dt8 = (tblock - tprev)/8.0d0
*       Avoid zero since TBLOCK = TPREV during the first block-step.
      IF (dt8.EQ.0.0d0) dt8 = step(ich)/8.0d0
*
*       Adopt the nearest truncated step (at most 8 subdivisions).
      IF (dt2.GT.0.0d0) THEN
          CALL stepk(dt2,dtn2)
          dtn = nint(dtn2/dt8)*dt8
      ELSE
*       Choose negative step if pericentre time < TPREV.
          dt2 = -dt2
          CALL stepk(dt2,dtn2)
          dtn = -nint(dtn2/dt8)*dt8
      END IF
*
*       Update time for new polynomial initializations (but check T - T0).
      time = tprev + dtn
*
*       Avoid prediction skip by XVPRED in case TIME - T0 = 0.0.
      IF (time - t0(ich).EQ.0.0d0) time = time + dt8/16.0d0
*
*       Re-define initial epoch for consistency (ignore phase error).
      t0s(isub) = time - timec
*
      IF (kz(30).GT.2) THEN
          WRITE (6,1)  time+toff, dt2, dt8
    1     FORMAT (' ABSORB:    TIME DT2 DT8 ',f12.6,1p,2e10.2)
      END IF
*
*       Increase membership of chain (JCLOSE: single body or KS pair).
      nch0 = nch
      CALL setsys
*
*       Improve coordinates & velocities of c.m. body to order F3DOT.
      CALL xvpred(ich,-1)
*
      sum = 0.0
      DO 5 k = 1,3
          xcm(k) = 0.0
          vcm(k) = 0.0
    5 CONTINUE
*
*       Accumulate mass-weighted moments of absorbed particle(s).
      DO 15 l = nch0+1,nch
          j = jlist(l)
          sum = sum + body(j)
          DO 10 k = 1,3
              xcm(k) = xcm(k) + body(j)*x(k,j)
              vcm(k) = vcm(k) + body(j)*xdot(k,j)
   10     CONTINUE
   15 CONTINUE
*
*       Form combined c.m. of old chain and new perturber(s).
      DO 20 k = 1,3
          xcm(k) = (body(ich)*x(k,ich) + xcm(k))/(body(ich) + sum)
          vcm(k) = (body(ich)*xdot(k,ich) + vcm(k))/(body(ich) + sum)
   20 CONTINUE
*
*       Define new relative coordinates & velocities and add to chain.
      lk = 3*nch0
      DO 30 l = nch0+1,nch
          j = jlist(l)
          SIZE(l) = radius(j)
          istar(l) = kstar(j)
          rij2 = 0.0
          DO 25 k = 1,3
              lk = lk + 1
              x4(k,l) = x(k,j) - xcm(k)
              xdot4(k,l) = xdot(k,j) - vcm(k)
              xch(lk) = x4(k,l)
              vch(lk) = xdot4(k,l)
              rij2 = rij2 + (x4(k,l) - x4(k,l-1))**2
   25     CONTINUE
*       Initialize new inverse distance(s) (some value needed in chpert.f).
          rinv(l-1) = 1.0/sqrt(rij2)
   30 CONTINUE
*
*       Re-define old chain variables with respect to new c.m.
      lk = 0
      DO 40 l = 1,nch0
          DO 35 k = 1,3
              lk = lk + 1
              xch(lk) = xch(lk) - (xcm(k) - x(k,ich))
              vch(lk) = vch(lk) - (vcm(k) - xdot(k,ich))
   35     CONTINUE
   40 CONTINUE
*
*       Create ghost particle(s) and remove from neighbour lists.
      DO 50 l = nch0+1,nch
          j = jlist(l)
          CALL ghost(j)
   50 CONTINUE
*
*       Update total mass and initialize new c.m. body variables.
      body(ich) = body(ich) + sum
      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
*
*       Remove ghost particle(s) from neighbour list of #ICH.
      jpert(1) = ich
      jlast = ntot + 1
      DO 70 l = nch0+1,nch
          jlist(1) = jlist(l)
          CALL nbrem(jlast,1,1)
   70 CONTINUE
*
*       Perform re-initialization of c.m. polynomials & perturber list.
      CALL reinit(isub)
*
*       Check optional output for centre of mass condition.
      IF (kz(30).GT.2) THEN
          DO 80 k = 1,6
              cm(k) = 0.0
   80     CONTINUE
*
          lk = 0
          DO 90 l = 1,nch
              DO 85 k = 1,3
                  lk = lk + 1
                  cm(k) = cm(k) + bodyc(l)*xch(lk)
                  cm(k+3) = cm(k+3) + bodyc(l)*vch(lk)
   85         CONTINUE
   90     CONTINUE
*
          DO 95 k = 1,6
              cm(k) = cm(k)/cm(7)
   95     CONTINUE
*
          WRITE (6,99)  (cm(k),k=1,6)
   99     FORMAT (' ABSORB:   CM ',1p,6e9.1)
      END IF
*
      RETURN
*
      END