Nbody6
chpert.f
Go to the documentation of this file.
1  SUBROUTINE chpert(GAMX)
2 *
3 *
4 * Final chain perturber selection.
5 * --------------------------------
6 *
7  include 'common6.h'
8  REAL*8 m,mass,mc,mij,mkk
9  parameter(nmx=10,nmx3=3*nmx,nmxm=nmx*(nmx-1)/2)
10  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
11  & zz(nmx3),wc(nmx3),mc(nmx),
12  & xj(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
13  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
14  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
15  & listc(lmax)
16  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
17  REAL*8 acc(nmx3),fp(nmx3)
18 *
19 *
20 * Skip on zero perturbers.
21  gamx = 0.0
22  IF (listc(1).EQ.0) go to 40
23 *
24 * Determine global values of XC first time (also called in REINIT).
25  IF (timec.EQ.0.0d0) THEN
26  CALL xcpred(0)
27  END IF
28 *
29 * Initialize the external perturbations.
30  nk = 3*nn
31  DO 1 k = 1,nk
32  acc(k) = 0.0d0
33  1 CONTINUE
34 *
35 * Copy provisional perturber list for evaluation loop.
36  npc = listc(1) + 1
37  DO 2 l = 2,npc
38  jpert(l) = listc(l)
39  2 CONTINUE
40 *
41 * Consider each provisional perturber in turn.
42  np = 1
43  DO 20 l = 2,npc
44  j = jpert(l)
45  ik = -3
46  itime = 0
47 * Sum perturber contributions over each chain component.
48  DO 10 i = 1,nn
49  ik = ik + 3
50  a1 = x(1,j) - xc(1,i)
51  a2 = x(2,j) - xc(2,i)
52  a3 = x(3,j) - xc(3,i)
53  rij2 = a1*a1 + a2*a2 + a3*a3
54  a6 = body(j)/(rij2*sqrt(rij2))
55  fp(ik+1) = a1*a6
56  fp(ik+2) = a2*a6
57  fp(ik+3) = a3*a6
58  acc(ik+1) = acc(ik+1) + fp(ik+1)
59  acc(ik+2) = acc(ik+2) + fp(ik+2)
60  acc(ik+3) = acc(ik+3) + fp(ik+3)
61 *
62 * Use standard perturbation test (GMIN) for acceptance.
63  IF (i.GT.1) THEN
64  df2 = 0.0
65  DO 5 k = 1,3
66  df = fp(ik+k) - fp(ik+k-3)
67  df2 = df2 + df**2
68  5 CONTINUE
69  gam = sqrt(df2)/((bodyc(i-1) + bodyc(i))*rinv(i-1)**2)
70 * Add accepted perturber to LISTC first time only.
71  IF (gam.GT.gmin.AND.itime.EQ.0) THEN
72  itime = 1
73  np = np + 1
74  listc(np) = j
75  END IF
76  END IF
77  10 CONTINUE
78  20 CONTINUE
79 *
80 * Save new membership.
81  listc(1) = np - 1
82 *
83 * Evaluate the total relative perturbations and save maximum.
84  DO 30 i = 1,nn-1
85  li = 3*(i - 1)
86  df2 = 0.0
87  DO 25 k = 1,3
88  df = acc(k+li+3) - acc(k+li)
89  df2 = df2 + df**2
90  25 CONTINUE
91 * Note that rejected perturbers are included in each final value.
92  gam = sqrt(df2)/((bodyc(i) + bodyc(i+1))*rinv(i)**2)
93  gamx = max(gam,gamx)
94  30 CONTINUE
95 *
96  40 RETURN
97  END