Nbody6
 All Files Functions Variables
core.f
Go to the documentation of this file.
1  SUBROUTINE core
2 *
3 *
4 * Density centre & core radius.
5 * -----------------------------
6 *
7  include 'common6.h'
8  REAL*8 rlist(lmax),rho,rhos
9  common/work1/ rho(nmax)
10 *
11 *
12 * Set generous cutoff.
13  rcore2 = rscale**2
14  IF (time.GT.0.0d0) rcore2 = max(9.0d0*rc**2,rcore2)
15 *
16 * Select > N/5 central particles.
17  5 nc = 0
18  DO 10 i = ifirst,ntot
19  ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
20  & (x(3,i) - rdens(3))**2
21  IF (ri2.LT.rcore2) THEN
22  nc = nc + 1
23  jlist(nc) = i
24  END IF
25  10 CONTINUE
26 *
27  IF (nc.LT.n/5) THEN
28  rcore2 = 1.5*rcore2
29  go to 5
30  END IF
31 *
32 * Obtain individual densities.
33  rhos = 0.0
34  DO 50 l = 1,nc
35  i = jlist(l)
36  nnb = list(1,i)
37 * Skip and set density to zero if there are no neighbours.
38  IF (nnb.EQ.0) THEN
39  rho(i) = 0.0
40  go to 50
41  END IF
42 * Estimate D6 to limit the candidates (additional factor of 2).
43  d6 = 2.0*(6.0/float(nnb))**0.66*rs(i)**2
44  xi = x(1,i)
45  yi = x(2,i)
46  zi = x(3,i)
47  20 n6 = 0
48 *
49  DO 25 lj = 1,nnb
50  j = list(lj+1,i)
51  rij2 = (xi - x(1,j))**2 + (yi - x(2,j))**2 +
52  & (zi - x(3,j))**2
53  IF (rij2.LT.d6) THEN
54  n6 = n6 + 1
55  ilist(n6) = j
56  rlist(n6) = rij2
57  END IF
58  25 CONTINUE
59 *
60  IF (n6.LT.6) THEN
61 * Make another iteration in rare cases.
62  d6 = 1.5*d6
63  IF (n6.LT.nnb) go to 20
64  END IF
65 *
66 * Sort list of square distances.
67  DO 40 ii = 1,n6
68  DO 35 jj = ii+1,n6
69  IF (rlist(jj).LT.rlist(ii)) THEN
70  rdum = rlist(ii)
71  idum = ilist(ii)
72  rlist(ii) = rlist(jj)
73  ilist(ii) = ilist(jj)
74  rlist(jj) = rdum
75  ilist(jj) = idum
76  END IF
77  35 CONTINUE
78  40 CONTINUE
79 *
80  i6 = min(n6,6)
81 * Calculate total mass of five nearest neighbours.
82  xmass = 0.0d0
83  DO 45 lj = 1,i6-1
84  xmass = xmass + body(ilist(lj))
85  45 CONTINUE
86 *
87  rho(i) = xmass/(rlist(i6)*sqrt(rlist(i6)))
88 * Define particle number density (not used).
89 * RHON = RHO(I)/ZMASS
90 * Assign zero weight if not enough neighbours.
91  IF (n6.LT.6) rho(i) = 0.0d0
92  rhos = max(rhos,rho(i))
93  50 CONTINUE
94 *
95 * Determine the density centre.
96  DO 60 k = 1,3
97  rdens(k) = 0.0d0
98  60 CONTINUE
99 *
100  rho1 = 0.0d0
101  DO 70 l = 1,nc
102  i = jlist(l)
103  DO 65 k = 1,3
104  rdens(k) = rdens(k) + rho(i)*x(k,i)
105  65 CONTINUE
106  rho1 = rho1 + rho(i)
107  70 CONTINUE
108 *
109 * Set current density centre based on improved determination.
110  rho1 = max(rho1,zmass/rscale**3)
111  DO 75 k = 1,3
112  rdens(k) = rdens(k)/rho1
113 * Ignore density centre for Plummer binary models except at end.
114  IF (kz(5).EQ.2.AND.ttot.LT.tcrit) THEN
115  rdens(k) = 0.0
116  END IF
117  75 CONTINUE
118 *
119 * Obtain density radius & averaged density.
120  rc = 0.0d0
121  rho2 = 0.0d0
122  DO 80 l = 1,nc
123  i = jlist(l)
124  xid = x(1,i) - rdens(1)
125  yid = x(2,i) - rdens(2)
126  zid = x(3,i) - rdens(3)
127  rid2 = xid**2 + yid**2 + zid**2
128  rc = rc + rho(i)**2*rid2
129  rho2 = rho2 + rho(i)**2
130  80 CONTINUE
131 *
132 * Form core radius and average & maximum density (scaled in OUTPUT).
133  IF (rho2.GT.0.0d0) rc = sqrt(rc/rho2)
134  rhod = rho2/rho1
135  rhod = (3.0/12.566)*rhod
136  rhom = (3.0/12.566)*rhos
137 * NB! Scaling factor 5 of Casertano & Hut eq. (V.3) replaced by XMASS.
138  IF (rc.EQ.0.0d0) THEN
139  rc = rscale
140  nc = n/2
141  END IF
142 *
143 * Set square core radius and its inverse for routine REGINT.
144  rc2 = rc**2
145  rc2in = 1.0/rc2
146 *
147 * Sum particles & mass inside the core radius and set rms velocity.
148  nc1 = 0
149  vc = 0.d00
150  zmc = 0.d00
151  DO 85 l = 1,nc
152  i = jlist(l)
153  xid = x(1,i) - rdens(1)
154  yid = x(2,i) - rdens(2)
155  zid = x(3,i) - rdens(3)
156  rid2 = xid**2 + yid**2 + zid**2
157  IF (rid2.LT.rc2) THEN
158  nc1 = nc1 + 1
159  vc = vc + xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
160  zmc = zmc + body(i)
161  END IF
162  85 CONTINUE
163 *
164 * Set core membership & rms velocity.
165  nc = max(nc1,2)
166  vc = sqrt(vc/float(nc))
167 *
168  RETURN
169 *
170  END