chfirr.f

Go to the documentation of this file.

      SUBROUTINE chfirr(I,IR,XI,XIDOT,FIRR,FD)
*
*
*       Irregular force & derivative on chain c.m.
*       ------------------------------------------
*
      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)
      REAL*8  xi(3),xidot(3),dx(3),dv(3),firr(3),fd(3),fp(3*ncmax),
     &        fpsum(3),fpd(3*ncmax),fpdsum(3)
      SAVE nskip
      DATA nskip /0/
*
*
*       Update perturber list at frequent intervals.
      IF (ir.GT.0) THEN
          CALL chlist(ich)
*       Note XC & UC for chain members predicted at new block-time.
      ELSE
          nskip = nskip + 1
          IF (nskip.GE.10) THEN
              CALL chlist(ich)
              nskip = 0
          END IF
      END IF
*
*       Subtract all close perturber forces acting on c.m. body #ICH.
      rpert2 = cmsep2*(0.5*rsum)**2
      npc = listc(1) + 1
      jdum = 0
      DO 30 l = 2,npc
          j = listc(l)
          dx(1) = x(1,j) - xi(1)
          dx(2) = x(2,j) - xi(2)
          dx(3) = x(3,j) - xi(3)
          rij2 = dx(1)**2 + dx(2)**2 + dx(3)**2
*       Skip correcting all interactions outside RPERT2.
          IF (rij2.GT.rpert2) go to 30
          IF (j.LE.n) go to 5
          jp = j - n
          j1 = 2*jp - 1
*       See whether c.m. approximation applies.
          IF (rij2.GT.cmsep2*r(jp)**2.OR.list(1,j1).EQ.0) go to 5
*       Resolve components of pair #JP.
          CALL ksres2(jp,j1,j2,rij2)
*       Note danger of low-order prediction (high-order KSRES2 in INTGRT).
          j = j1
          jdum = j1
*       Sum over single body or individual KS components.
    5     dr2 = 0.0
          drdv = 0.0
          DO 10 k = 1,3
              dx(k) = x(k,j) - xi(k)
              dv(k) = xdot(k,j) - xidot(k)
          dr2 = dr2 + dx(k)**2
          drdv = drdv + dx(k)*dv(k)
   10     CONTINUE
          dr2i = 1.0/dr2
          dr3i = body(j)*dr2i*sqrt(dr2i)
      drdv = 3.0*drdv*dr2i
*
*       Subtract force and first derivative.
          DO 20 k = 1,3
              firr(k) = firr(k) - dx(k)*dr3i
              fd(k) = fd(k) - (dv(k) - dx(k)*drdv)*dr3i
   20     CONTINUE
          IF (j.EQ.jdum) THEN
              j = j + 1
              go to 5
          END IF
   30 CONTINUE
*
*       Initialize the perturbing force components and vectorial sum.
      DO 35 k = 1,3*nch
          fp(k) = 0.0d0
          fpd(k) = 0.0d0
   35 CONTINUE
      DO 40 k = 1,3
          fpsum(k) = 0.0d0
          fpdsum(k) = 0.0d0
   40 CONTINUE
*
*       Sum each perturber contribution over the chain components.
      kdum = 0
      im1 = 0
      DO 70 ll = 2,npc
          k = listc(ll)
          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
*
*       Skip perturbers consistent with c.m. approximation test above.
          IF (rij2.GT.rpert2) go to 70
*       Decide appropriate summation (components or c.m. approximation).
          IF (k.GT.n) THEN
*       Resolve the binary inside 20*R (OK for tidal effect on chain).
              IF (rij2.LT.400.0*r(k-n)**2) THEN
                  kdum = 2*(k - n) - 1
                  k = kdum
                  go to 55
              END IF
          END IF
*
*       Obtain perturbation on each component.
          DO 50 im = 1,nch
              im1 = 3*(im - 1)
   46         dr2 = 0.0
              drdv = 0.0
              DO 48 l = 1,3
                  dx(l) = x(l,k) - xc(l,im)
                  dv(l) = xdot(l,k) - uc(l,im)
              dr2 = dr2 + dx(l)**2
              drdv = drdv + dx(l)*dv(l)
   48         CONTINUE
              dr2i = 1.0/dr2
              dr3i = body(k)*dr2i*sqrt(dr2i)
          drdv = 3.0*drdv*dr2i
*
*       Accumulate perturbing force & derivative.
              DO 49 l = 1,3
                  fp(im1+l) = fp(im1+l) + dx(l)*dr3i
                  fpd(im1+l) = fpd(im1+l) + (dv(l) - dx(l)*drdv)*dr3i
   49         CONTINUE
*
              IF (k.EQ.kdum) THEN
                  k = k + 1
                  go to 46
              ELSE
                  IF (kdum.GT.0) k = kdum
              END IF
   50     CONTINUE
*
*       Reset dummy index after use (otherwise bug with two KS pairs).
          kdum = 0
          go to 70
*
*       Sum over individual components of pair #J using c.m. approximation.
   55     dr2 = 0.0
          drdv = 0.0
          DO 60 l = 1,3
              dx(l) = x(l,k) - xi(l)
              dv(l) = xdot(l,k) - xidot(l)
          dr2 = dr2 + dx(l)**2
          drdv = drdv + dx(l)*dv(l)
   60     CONTINUE
          dr2i = 1.0/dr2
          dr3i = body(k)*dr2i*sqrt(dr2i)
      drdv = 3.0*drdv*dr2i
*
*       Add force & first derivative.
          DO 65 l = 1,3
              firr(l) = firr(l) + dx(l)*dr3i
              fd(l) = fd(l) + (dv(l) - dx(l)*drdv)*dr3i
   65     CONTINUE
          IF (k.EQ.kdum) THEN
              k = k + 1
              go to 55
          END IF
   70 CONTINUE
*
*       Add any perturbations on components to the c.m. force & derivative.
      IF (im1.GT.0) THEN
*       Sum individual perturbations.
          DO 80 l = 1,nch
              l1 = 3*(l - 1)
              DO 75 k = 1,3
                  fpsum(k) = fpsum(k) + bodyc(l)*fp(l1+k)
                  fpdsum(k) = fpdsum(k) + bodyc(l)*fpd(l1+k)
   75         CONTINUE
   80     CONTINUE
*
*       Include mass-weighted contributions.
          bodyin = 1.0/body(i)
          DO 90 k = 1,3
              firr(k) = fpsum(k)*bodyin + firr(k)
              fd(k) = fpdsum(k)*bodyin + fd(k)
   90     CONTINUE
      END IF
*
      RETURN
*
      END