Nbody6
search.f
Go to the documentation of this file.
1  SUBROUTINE search(I,IKS)
2 *
3 *
4 * Close encounter search.
5 * -----------------------
6 *
7  include 'common6.h'
8  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
9  & names(ncmax,5),isys(5)
10 *
11 *
12 * Increase counter for regularization attempts.
13  nkstry = nkstry + 1
14 *
15  fmax = 0.0
16  nclose = 0
17 * Find dominant neighbour by selecting all STEP(J) < 2*DTMIN.
18  l = list(1,i) + 2
19  2 l = l - 1
20  IF (l.LT.2) go to 10
21 *
22 * First see whether any c.m. with small step is within 2*RMIN.
23  IF (list(l,i).LE.n) go to 4
24 *
25  j = list(l,i)
26  IF (step(j).GT.smin) go to 2
27  a1 = x(1,j) - x(1,i)
28  a2 = x(2,j) - x(2,i)
29  a3 = x(3,j) - x(3,i)
30  rij2 = a1*a1 + a2*a2 + a3*a3
31  IF (rij2.GT.rmin22) go to 2
32 *
33  fij = body(j)/rij2
34  IF (fmax.LT.fij) fmax = fij
35 * Abandon further search if c.m. force exceeds half the total force.
36  IF (fmax**2.LT.f(1,i)**2 + f(2,i)**2 + f(3,i)**2) THEN
37  nclose = nclose + 1
38  jlist(nclose) = j
39  go to 2
40  ELSE
41  go to 10
42  END IF
43 *
44 * Continue searching single particles with current value of FMAX.
45  4 jcomp = 0
46 *
47  DO 6 k = l,2,-1
48  j = list(k,i)
49  IF (step(j).GT.smin) go to 6
50  a1 = x(1,j) - x(1,i)
51  a2 = x(2,j) - x(2,i)
52  a3 = x(3,j) - x(3,i)
53  rij2 = a1*a1 + a2*a2 + a3*a3
54  IF (rij2.LT.rmin22) THEN
55  nclose = nclose + 1
56  jlist(nclose) = j
57 * Remember index of every single body with small step inside 2*RMIN.
58  fij = (body(i) + body(j))/rij2
59  IF (fij.GT.fmax) THEN
60  fmax = fij
61 * Save square distance and global index of dominant body.
62  rjmin2 = rij2
63  jcomp = j
64  END IF
65  END IF
66  6 CONTINUE
67 *
68 * See whether dominant component is a single particle inside RMIN.
69  IF (jcomp.LT.ifirst.OR.jcomp.GT.n) go to 10
70 * Accept one single candidate inside 2*RMIN (which makes PERT = 0).
71  IF (rjmin2.GT.rmin2.AND.nclose.GT.1) go to 10
72 *
73  rdot = (x(1,i) - x(1,jcomp))*(xdot(1,i) - xdot(1,jcomp)) +
74  & (x(2,i) - x(2,jcomp))*(xdot(2,i) - xdot(2,jcomp)) +
75  & (x(3,i) - x(3,jcomp))*(xdot(3,i) - xdot(3,jcomp))
76 *
77 * Only select approaching particles (include nearly circular case).
78  rijmin = sqrt(rjmin2)
79  IF (rdot.GT.0.02*sqrt((body(i) + body(jcomp))*rijmin)) go to 10
80 *
81 * Ensure a massive neighbour is included in perturbation estimate.
82  bcm = body(i) + body(jcomp)
83  IF (body1.GT.10.0*bcm) THEN
84  jbig = 0
85  big = bcm
86  nnb1 = list(1,i) + 1
87  DO 20 l = 2,nnb1
88  j = list(l,i)
89  IF (body(j).GT.big) THEN
90  jbig = j
91  big = body(j)
92  END IF
93  20 CONTINUE
95  DO 25 l = 1,nclose
96  IF (jlist(l).EQ.jbig) THEN
97  jbig = 0
98  END IF
99  25 CONTINUE
100  IF (jbig.GT.0) THEN
101  nclose = nclose + 1
102  jlist(nclose) = jbig
103  END IF
104  END IF
105 *
106 * Evaluate vectorial perturbation due to the close bodies.
107  CALL fpert(i,jcomp,nclose,pert)
108 *
109 * Accept #I & JCOMP if the relative motion is dominant (GI < 0.25).
110  gi = pert*rjmin2/bcm
111  IF (gi.GT.0.25) THEN
112 * IF (KZ(4).GT.0.AND.TIME-TLASTT.GT.4.44*TCR/FLOAT(N))
113 * & CALL EVOLVE(JCOMP,0)
114  go to 10
115  END IF
116 *
117 * Exclude any c.m. body of compact subsystem (TRIPLE, QUAD or CHAIN).
118  DO 8 isub = 1,nsub
119  namei = names(1,isub)
120  IF (namei.EQ.name(i).OR.namei.EQ.name(jcomp)) go to 10
121  8 CONTINUE
122 *
123 * Also check possible c.m. body of chain regularization (NAME = 0).
124  IF (nch.GT.0) THEN
125  IF (name(i).EQ.0.OR.name(jcomp).EQ.0) go to 10
126  END IF
127 *
128 * Save index and increase indicator to denote new regularization.
129  icomp = i
130  iks = iks + 1
131 *
132  10 RETURN
133 *
134  END