chpert.f

Go to the documentation of this file.

      SUBROUTINE chpert(GAMX)
*
*
*       Final chain perturber selection.
*       --------------------------------
*
      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/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      REAL*8  acc(nmx3),fp(nmx3)
*
*
*       Skip on zero perturbers.
      gamx = 0.0
      IF (listc(1).EQ.0) go to 40
*
*       Determine global values of XC first time (also called in REINIT).
      IF (timec.EQ.0.0d0) THEN
          CALL xcpred(0)
      END IF
*
*       Initialize the external perturbations.
      nk = 3*nn
      DO 1 k = 1,nk
          acc(k) = 0.0d0
    1 CONTINUE
*
*       Copy provisional perturber list for evaluation loop.
      npc = listc(1) + 1
      DO 2 l = 2,npc
          jpert(l) = listc(l)
    2 CONTINUE
*
*       Consider each provisional perturber in turn.
      np = 1
      DO 20 l = 2,npc
          j = jpert(l)
          ik = -3
          itime = 0
*       Sum perturber contributions over each chain component.
          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
              a6 = body(j)/(rij2*sqrt(rij2))
              fp(ik+1) = a1*a6
              fp(ik+2) = a2*a6
              fp(ik+3) = a3*a6
              acc(ik+1) = acc(ik+1) + fp(ik+1)
              acc(ik+2) = acc(ik+2) + fp(ik+2)
              acc(ik+3) = acc(ik+3) + fp(ik+3)
*
*       Use standard perturbation test (GMIN) for acceptance.
              IF (i.GT.1) THEN
                  df2 = 0.0
                  DO 5 k = 1,3
                      df = fp(ik+k) - fp(ik+k-3)
                      df2 = df2 + df**2
    5             CONTINUE
                  gam = sqrt(df2)/((bodyc(i-1) + bodyc(i))*rinv(i-1)**2)
*       Add accepted perturber to LISTC first time only.
                  IF (gam.GT.gmin.AND.itime.EQ.0) THEN
                      itime = 1
                      np = np + 1
                      listc(np) = j
                  END IF
              END IF
   10     CONTINUE
   20 CONTINUE
*
*       Save new membership.
      listc(1) = np - 1
*
*       Evaluate the total relative perturbations and save maximum.
      DO 30 i = 1,nn-1
          li = 3*(i - 1)
          df2 = 0.0
          DO 25 k = 1,3
              df = acc(k+li+3) - acc(k+li)
              df2 = df2 + df**2
   25     CONTINUE
*       Note that rejected perturbers are included in each final value.
          gam = sqrt(df2)/((bodyc(i) + bodyc(i+1))*rinv(i)**2)
          gamx = max(gam,gamx)
   30 CONTINUE
*
   40 RETURN
      END