xtpert.f

Go to the documentation of this file.

      SUBROUTINE xtpert(ACC,CHTIME)
*
*
*       External perturbations on chain members.
*       ----------------------------------------
*
      include 'common6.h'
      REAL*8  m,mass,mc,mij,mkk
      parameter(nmx=10,nmx3=3*nmx,nmxm=nmx*(nmx-1)/2)
      common/chain1/  xch(nmx3),vch(nmx3),m(nmx),
     &                zz(nmx3),wc(nmx3),mc(nmx),
     &                xj(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/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
      REAL*8  acc(nmx3)
*
*
*       Predict current coordinates of perturbers & chain members.
      time0 = time
      isub = isys(5)
      time = chtime + t0s(isub)
      CALL xcpred(1)
*
*       Initialize the external perturbations.
      nk = 3*nn
      DO 1 k = 1,nk
          acc(k) = 0.0d0
    1 CONTINUE
*
*       Consider each perturber in turn.
      npc = listc(1) + 1
      DO 20 l = 2,npc
          j = listc(l)
*       Sum perturber contributions over each chain component.
          ik = -3
          ks = 0
          DO 10 i = 1,nn
              ik = ik + 3
              a1 = x(1,j) - xc(1,i)
              a2 = x(2,j) - xc(2,i)
              a3 = x(3,j) - xc(3,i)
              rij2 = a1*a1 + a2*a2 + a3*a3
              IF (j.LE.n) go to 5
              jpair = j - n
*       Check c.m. approximation (only resolve KS components once).
              IF (rij2.GT.cmsep2*r(jpair)**2) go to 5
              IF (ks.EQ.0) THEN
                  CALL ksres(jpair,j1,j2,rij2)
                  ks = 1
              END IF
              kdum = j1
              k = kdum
    4         a1 = x(1,k) - xc(1,i)
              a2 = x(2,k) - xc(2,i)
              a3 = x(3,k) - xc(3,i)
              rij2 = a1*a1 + a2*a2 + a3*a3
              a6 = body(k)/(rij2*sqrt(rij2))
              acc(ik+1) = acc(ik+1) + a1*a6
              acc(ik+2) = acc(ik+2) + a2*a6
              acc(ik+3) = acc(ik+3) + a3*a6
*       See whether the second component has been included.
              IF (k.EQ.kdum) THEN
                  k = k + 1
                  go to 4
              END IF
              go to 10
*
    5         a6 = body(j)/(rij2*sqrt(rij2))
              acc(ik+1) = acc(ik+1) + a1*a6
              acc(ik+2) = acc(ik+2) + a2*a6
              acc(ik+3) = acc(ik+3) + a3*a6
   10     CONTINUE
   20 CONTINUE
*
*       Restore the global time.
      time = time0
*
      RETURN
      END