kspert.f

Go to the documentation of this file.

      SUBROUTINE kspert(I1,NNB0,XI,VI,FP,FD)
*
*
*       Perturbation on KS pair.
*       ------------------------
*
      include 'common6.h'
      common/chainc/  xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
     &                listc(lmax)
      REAL*8  xi(6),vi(6),fp(6),fd(6),tf(3),td(3)
*
*
*       Initialize the perturbing force & first derivative.
      DO 10 k = 1,6
          fp(k) = 0.0d0
          fd(k) = 0.0d0
   10 CONTINUE
*
*       Set index of the last single perturber.
      nnb2 = nnb0 + 1
   15 IF (list(nnb2,i1).LE.n) go to 20
      nnb2 = nnb2 - 1
      IF (nnb2.GT.1) go to 15
*       Include special case of only c.m. perturbers.
      go to 30
*
*       Obtain the perturbation from single particles.
   20 DO 25 l = 2,nnb2
          k = list(l,i1)
          a1 = x(1,k) - xi(1)
          a2 = x(2,k) - xi(2)
          a3 = x(3,k) - xi(3)
          rij2 = a1*a1 + a2*a2 + a3*a3
          a6 = body(k)/(rij2*sqrt(rij2))
          fp(1) = fp(1) + a1*a6
          fp(2) = fp(2) + a2*a6
          fp(3) = fp(3) + a3*a6
          v1 = xdot(1,k) - vi(1)
          v2 = xdot(2,k) - vi(2)
          v3 = xdot(3,k) - vi(3)
          a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
          fd(1) = fd(1) + (v1 - a1*a9)*a6
          fd(2) = fd(2) + (v2 - a2*a9)*a6
          fd(3) = fd(3) + (v3 - a3*a9)*a6
*       Perturbation on first component.
*
          a1 = x(1,k) - xi(4)
          a2 = x(2,k) - xi(5)
          a3 = x(3,k) - xi(6)
          rij2 = a1*a1 + a2*a2 + a3*a3
          a6 = body(k)/(rij2*sqrt(rij2))
          fp(4) = fp(4) + a1*a6
          fp(5) = fp(5) + a2*a6
          fp(6) = fp(6) + a3*a6
          v1 = xdot(1,k) - vi(4)
          v2 = xdot(2,k) - vi(5)
          v3 = xdot(3,k) - vi(6)
          a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
          fd(4) = fd(4) + (v1 - a1*a9)*a6
          fd(5) = fd(5) + (v2 - a2*a9)*a6
          fd(6) = fd(6) + (v3 - a3*a9)*a6
*       Perturbation on second component.
   25 CONTINUE
*
*       See whether to include any remaining c.m. perturbers.
      IF (nnb2.GT.nnb0) go to 40
*
   30 kdum = 0
*       Dummy index to enable summation of c.m. or resolved components.
      nnb3 = nnb2 + 1
      DO 35 l = nnb3,nnb0+1
          k = list(l,i1)
          a1 = x(1,k) - xi(1)
          a2 = x(2,k) - xi(2)
          a3 = x(3,k) - xi(3)
          rij2 = a1*a1 + a2*a2 + a3*a3
*       See whether c.m. approximation applies (ignore unperturbed case).
          j = k - n
          IF (rij2.GT.cmsep2*r(j)**2) THEN
              v1 = xdot(1,k) - vi(1)
              v2 = xdot(2,k) - vi(2)
              v3 = xdot(3,k) - vi(3)
              a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
              go to 33
          END IF
*
*       Resolve pair #J and sum over individual components.
          CALL ksres2(j,j1,j2,rij2)
          kdum = j1
          k = kdum
   32     a1 = x(1,k) - xi(1)
          a2 = x(2,k) - xi(2)
          a3 = x(3,k) - xi(3)
          rij2 = a1*a1 + a2*a2 + a3*a3
          v1 = xdot(1,k) - vi(1)
          v2 = xdot(2,k) - vi(2)
          v3 = xdot(3,k) - vi(3)
          a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
   33     a6 = body(k)/(rij2*sqrt(rij2))
          fp(1) = fp(1) + a1*a6
          fp(2) = fp(2) + a2*a6
          fp(3) = fp(3) + a3*a6
          fd(1) = fd(1) + (v1 - a1*a9)*a6
          fd(2) = fd(2) + (v2 - a2*a9)*a6
          fd(3) = fd(3) + (v3 - a3*a9)*a6
*       Perturbation on first component.
*
          a1 = x(1,k) - xi(4)
          a2 = x(2,k) - xi(5)
          a3 = x(3,k) - xi(6)
          rij2 = a1*a1 + a2*a2 + a3*a3
          a6 = body(k)/(rij2*sqrt(rij2))
          fp(4) = fp(4) + a1*a6
          fp(5) = fp(5) + a2*a6
          fp(6) = fp(6) + a3*a6
          v1 = xdot(1,k) - vi(4)
          v2 = xdot(2,k) - vi(5)
          v3 = xdot(3,k) - vi(6)
          a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
          fd(4) = fd(4) + (v1 - a1*a9)*a6
          fd(5) = fd(5) + (v2 - a2*a9)*a6
          fd(6) = fd(6) + (v3 - a3*a9)*a6
*       Perturbation on second component.
*
          IF (k.EQ.kdum) THEN
              k = k + 1
              go to 32
          END IF
   35 CONTINUE
*
*       Check perturbation correction due to regularized chain.
   40 IF (nch.GT.0) THEN
          DO 45 l = 2,nnb2
              j = list(l,i1)
              IF (j.GT.ich) go to 50
              IF (j.EQ.ich) THEN
                  j1 = i1
                  CALL fchain(j1,1,xi(1),vi(1),fp(1),fd(1))
                  j1 = j1 + 1
                  CALL fchain(j1,1,xi(4),vi(4),fp(4),fd(4))
                  go to 50
              END IF
   45     CONTINUE
      END IF 
*
*       Set the relative perturbing force and first derivative.
   50 DO 55 k = 1,3
          fp(k) = fp(k) - fp(k+3)
          fd(k) = fd(k) - fd(k+3)
          tf(k) = 0.0d0
          td(k) = 0.0d0
   55 CONTINUE
*
*       See whether the linearized perturbation should be included.
      IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
          q1 = xi(1) - xi(4)
          q3 = xi(3) - xi(6)
          CALL xtrnlp(q1,q3,tf)
*
*       Use same formalism for the first derivative (omit Coriolis force).
          vx = vi(1) - vi(4)
          vz = vi(3) - vi(6)
          CALL xtrnlp(vx,vz,td)
          DO 60 k = 1,3
             fp(k) = fp(k) + tf(k)
             fd(k) = fd(k) + td(k)
   60     CONTINUE
      END IF
*
*       Check optional Plummer potential.
      IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
          ri2 = ap2
          rrdot = 0.0
*       Form one central distance and scalar product of relative motion.
          DO 65 k = 1,3
              ri2 = ri2 + xi(k)**2
              rrdot = rrdot + (xi(k) - xi(k+3))*(vi(k) - vi(k+3))
   65     CONTINUE
          zf = 1.0/ri2
*       Write current mass inside RI as MP*R3*ZF^{3/2} (Heggie & Hut p.73).
          fmp = mp*zf*sqrt(zf)
          DO 70 k = 1,3
              xrel = xi(k) - xi(k+3)
              vrel = vi(k) - vi(k+3)
              fp(k) = fp(k) - xrel*fmp
              fd(k) = fd(k) - (vrel - 3.0*rrdot*zf*xrel)*fmp
   70     CONTINUE
      END IF
*
      RETURN
*
      END