Nbody6
chlist.f
Go to the documentation of this file.
1  SUBROUTINE chlist(I)
2 *
3 *
4 * Perturber list for chain regularization.
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  & xi(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/cpert/ rgrav,gpert,ipert,npert
17 *
18 *
19 * Prepare perturber selection using maximum of RSUM & RGRAV.
20  rpert = max(rsum,rgrav)
21 * Restrict effective size to RMIN (energy may be near zero).
22  rpert = min(rpert,rmin)
23  rcrit2 = 2.0*rpert**2/body(i)
24  rcrit3 = rcrit2*rpert/gmin
25 * Base fast search on maximum binary mass (2*BODY1).
26  rcrit2 = 2.0*rcrit2*body1*cmsep2
27  pmax = 0.0
28 *
29 * Select new perturbers from chain c.m. neighbour list.
30  nnb1 = 1
31  nnb2 = list(1,i) + 1
32  DO 10 l = 2,nnb2
33  j = list(l,i)
34  w1 = x(1,j) - x(1,i)
35  w2 = x(2,j) - x(2,i)
36  w3 = x(3,j) - x(3,i)
37  rsep2 = w1*w1 + w2*w2 + w3*w3
38 * Include any merged c.m. bodies in the fast test.
39  IF (rsep2.LT.rcrit2.OR.name(j).LT.0) THEN
40  rij3 = rsep2*sqrt(rsep2)
41 * Estimate perturbed distance from tidal limit approximation.
42  IF (rij3.LT.body(j)*rcrit3) THEN
43  nnb1 = nnb1 + 1
44  listc(nnb1) = j
45  pmax = max(body(j)/rij3,pmax)
46  END IF
47  END IF
48  10 CONTINUE
49 *
50 * Save perturber membership.
51  listc(1) = nnb1 - 1
52 *
53 * Form current perturbation from nearest body (effective value RPERT).
54  gpert = 2.0*pmax*rpert**3/body(ich)
55 *
56 * Obtain actual perturbation (note GAMX zero first time).
57  CALL chpert(gamx)
58 *
59 * WRITE (6,20) ICH, NNB2-1, LISTC(1), RSUM, RPERT, GPERT, GAMX
60 * 20 FORMAT (' CHLIST: ICH NB NP RSUM RPERT GPERT GX ',3I4,1P,4E9.1)
61 *
62  IF (gamx.GT.0.0d0) THEN
63  gpert = gamx
64  END IF
65  npert = listc(1)
66
67 * Ensure consistent perturbation indicator after re-determination.
68  IF (gpert.GT.1.0d-06.AND.npert.GT.0) THEN
69  ipert = 1
70  ELSE
71  ipert = 0
72  END IF
73 *
74  RETURN
75 *
76  END