Nbody6
lagr2.f
Go to the documentation of this file.
1  SUBROUTINE lagr2(C)
2 *
3 *
4 * Mass distribution for two mass groups.
5 * --------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 r2,zmi(nmax)
9  common/work1/ r2(nmax)
10  parameter(lx=11)
11  REAL*8 c(3),flagr(lx),rlagr(lx)
12 * DATA FLAGR/-1.9,-1.7,-1.5,-1.3,-1.1,-.9,-.7,-.5,-.3,-.1/
13 * DATA FLAGR/0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,
14 * & 0.75,0.9/
15  DATA flagr/0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.625,0.75,0.99/
16 *
17 *
18 * Set square radii of single particles & c.m. bodies (NAME <= NZERO/5).
19  nam1 = nzero/5
20 * Apply the mass test in same proportion for two Plummer spheres.
21  IF (kz(5).EQ.2) THEN
22  nam1 = (nzero - n1)/5
23  nam2 = n1 + (nzero - n1)/5
24  END IF
25  iter = 0
26  nblack = 0
27  DO 5 i = 1,n
28  IF (kstar(i).EQ.14) nblack = nblack + 1
29  5 CONTINUE
30 *
31  1 np = 0
32  zm1 = 0.0
33  DO 10 i = ifirst,ntot
34  IF (kz(5).EQ.1) THEN
35  IF (name(i).GT.nam1.AND.name(i).LE.nzero) go to 10
36  ELSE IF (kz(5).EQ.2) THEN
37  IF ((name(i).GT.nam1.AND.name(i).LE.n1).OR.
38  & (name(i).GT.nam2.AND.name(i).LE.nzero)) go to 10
39  END IF
40  zm1 = zm1 + body(i)
41  np = np + 1
42  r2(np) = (x(1,i) - c(1))**2 + (x(2,i) - c(2))**2 +
43  & (x(3,i) - c(3))**2
44  jlist(np) = i
45  ilist(np) = i
46  zmi(np) = body(i)
47  10 CONTINUE
48 *
49 * Improve choice of NAM1 using relative deviation (max 10 tries).
50  dm = (zm1 - 0.5*zmass)/zmass
51  IF (abs(dm).GT.0.001.AND.iter.LE.10) THEN
52  iter = iter + 1
53  IF (abs(dm).GT.0.05) dm = dm*0.05/abs(dm)
54  nm = dm*float(n-npairs)
55  IF (nm.NE.0) THEN
56  nam1 = nam1 - nm
57  IF (kz(5).EQ.2) THEN
58  nam2 = nam2 - nm*(nzero - n1)/(5*nam2)
59  END IF
60  go to 1
61  END IF
62  END IF
63 *
64 * Obtain sum of inverse separations to evaluate mass segregation.
65  pot1 = 0.0
66  DO 14 l = 1,np-1
67  i = jlist(l)
68  DO 12 ll = l+1,np
69  j = jlist(ll)
70  rij2 = 0.0
71  DO 11 k = 1,3
72  rij2 = rij2 + (x(k,i) - x(k,j))**2
73  11 CONTINUE
74  pot1 = pot1 + 1.0/sqrt(rij2)
75  12 CONTINUE
76  14 CONTINUE
77 * Define mean geometric distance.
78  rp1 = float(np)*(np - 1)/(2.0*pot1)
79 *
80 * Sort square distances with respect to the centre C.
81  CALL sort1(np,r2,jlist)
82
83 * Determine Lagrangian radii for specified mass fractions.
84  DO 20 il = 1,lx
85  zm = 0.0
86  zmh = flagr(il)*zm1
87  i = 0
88  15 i = i + 1
89  im = jlist(i)
90  zm = zm + body(im)
91  IF (zm.LT.zmh) go to 15
92  rlagr(il) = sqrt(r2(i))
93  20 CONTINUE
94 *
95 * Obtain half-mass radius separately.
96  zm = 0.0
97  zmh = 0.5*zm1
98  i = 0
99  25 i = i + 1
100  im = jlist(i)
101  zm = zm + body(im)
102  IF (zm.LT.zmh) go to 25
103
104  rh1 = sqrt(r2(i))
105  np1 = np
106 *
107  WRITE (31,30) time+toff, (log10(rlagr(k)),k=1,lx)
108  30 FORMAT (' ',f8.1,13f7.3)
109  CALL flush(31)
110 *
111 * Sort masses from first group in increasing order.
112  CALL sort1(np,zmi,ilist)
113  zmx1 = zmi(1)
114  zmx0 = 0.5*(zmx1 + body1)
115  itry = 0
116  32 np0 = 0
117  lx0 = 0
118 * Count dominant masses above half minimum + maximum single stars.
119  DO 35 l = 1,np
120  IF (zmi(l).GT.zmx0) THEN
121  IF (lx0.EQ.0) lx0 = l
122  np0 = np0 + 1
123  END IF
124  35 CONTINUE
125 * Ensure at least 1% in the most massive bin or use NBLACK > 2.
126  itry = itry + 1
127  npcrit = sqrt(float(n - npairs))
128  npcrit = min(npcrit,10)
129  IF (nblack.GT.2) npcrit = nblack
130  IF (np0.LT.npcrit.AND.itry.LT.20) THEN
131  zmx0 = 0.9*zmx0
132  go to 32
133  END IF
134  IF (np0.GT.npcrit.AND.itry.LT.20) THEN
135  zmx0 = 1.1*zmx0
136  go to 32
137  END IF
138 *
139 * Form sum of inverse separations to evaluate mass segregation.
140  pot0 = 0.0
141  lx0 = max(lx0,1)
142  DO 44 l = lx0,np-1
143  i = ilist(l)
144  DO 42 ll = l+1,np
145  j = ilist(ll)
146  rij2 = 0.0
147  DO 41 k = 1,3
148  rij2 = rij2 + (x(k,i) - x(k,j))**2
149  41 CONTINUE
150  pot0 = pot0 + 1.0/sqrt(rij2)
151  42 CONTINUE
152  44 CONTINUE
153 * Define mean geometric distance.
154  IF (np0.LE.1) pot0 = 1.0
155  rp0 = float(np0)*(np0 - 1)/(2.0*pot0)
156 *
157 * Search for single black holes and include KS binary.
158  n14 = 0
159  n14x = 0
160  semi = 0
161  potbh = 0.0
162  zmbh = 0.0
163  DO 45 i = ifirst,ntot
164  IF (kstar(i).EQ.14.AND.body(i).GT.0.0d0) THEN
165  n14 = n14 + 1
166  ilist(n14) = i
167  zmbh = zmbh + body(i)*smu
168  ELSE IF (i.GT.n) THEN
169 * Note: the following lines must be used for KS binary.
170  j1 = 2*(i - n) - 1
171  IF (kstar(j1).EQ.14.AND.kstar(j1+1).EQ.14) THEN
172 * Add only one KS component to avoid small RIJ2 but use extra count.
173  n14 = n14 + 1
174  n14x = n14x + 1
175  ilist(n14) = i
176  zmbh = zmbh + body(i)*smu
177  ip = i - n
178  semi = -0.5*body(i)/h(ip)
179  END IF
180  END IF
181  45 CONTINUE
182 *
183 * Form sum of inverse separations to evaluate mass segregation.
184  pot12 = 0.0
185  DO 48 l = 1,n14-1
186  i = ilist(l)
187  DO 47 ll = l+1,n14
188  j = ilist(ll)
189  rij2 = 0.0
190  DO 46 k = 1,3
191  rij2 = rij2 + (x(k,i) - x(k,j))**2
192  46 CONTINUE
193  potbh = potbh + 1.0/sqrt(rij2)
194 * Note internal terms not included at present.
195  IF (i.LT.ifirst) pot12 = pot12 + 1.0/sqrt(rij2)
196  47 CONTINUE
197  48 CONTINUE
198 * Define mean geometric distance (modified by binary 10/2/11).
199  n14 = n14 + n14x
200  IF (n14.GT.1) rbh = float(n14)*(n14 - 1)/(2.0*potbh)
201 * Omit internal binary contribution for secondary definition.
202  IF (pot12.GT.0.0) THEN
203  rx = float(n14-1)*(n14 - 2)/(2.0*(potbh-pot12))
204  ELSE
205  rx = 0.0
206  END IF
207 *
208 * Search for single neutron stars and include KS binary.
209  n13 = 0
210  semi2 = 0
211  pot13 = 0.0
212  r13 = 0.0
213  zm13 = 0.0
214  DO 145 i = ifirst,ntot
215  IF (kstar(i).EQ.13) THEN
216  n13 = n13 + 1
217  ilist(n13) = i
218  zm13 = zm13 + body(i)*smu
219  ELSE IF (i.GT.n) THEN
220  j1 = 2*(i - n) - 1
221  IF (kstar(j1).EQ.13.AND.kstar(j1+1).EQ.13) THEN
222  n13 = n13 + 1
223  ilist(n13) = i
224  zm13 = zm13 + body(i)*smu
225  ip = i - n
226  semi2 = -0.5*body(i)/h(ip)
227  END IF
228  END IF
229  145 CONTINUE
230 *
231 * Form sum of inverse separations to evaluate mass segregation.
232  DO 148 l = 1,n13-1
233  i = ilist(l)
234  DO 147 ll = l+1,n13
235  j = ilist(ll)
236  rij2 = 0.0
237  DO 146 k = 1,3
238  rij2 = rij2 + (x(k,i) - x(k,j))**2
239  146 CONTINUE
240  pot13 = pot13 + 1.0/sqrt(rij2)
241  147 CONTINUE
242  148 CONTINUE
243 * Define mean geometric distance (modified by binary 10/2/11).
244  IF (n13.GT.1) r13 = float(n13)*(n13 - 1)/(2.0*pot13)
245 *
246 * Treat the low-mass particles in the same way (exclude c.m. bodies).
247  np = 0
248  zm2 = 0.0
249  DO 50 i = ifirst,n
250  IF (kz(5).EQ.1) THEN
251  IF (name(i).LE.nam1) go to 50
252  ELSE IF (kz(5).EQ.2) THEN
253  IF (name(i).LE.nam1.OR.
254  & (name(i).GT.n1.AND.name(i).LE.nam2)) go to 50
255  END IF
256  zm2 = zm2 + body(i)
257  np = np + 1
258  r2(np) = (x(1,i) - c(1))**2 + (x(2,i) - c(2))**2 +
259  & (x(3,i) - c(3))**2
260  jlist(np) = i
261  50 CONTINUE
262 *
263 * Treat second mass group in the same way.
264  pot2 = 0.0
265  DO 54 l = 1,np-1
266  i = jlist(l)
267  DO 52 ll = l+1,np
268  j = jlist(ll)
269  rij2 = 0.0
270  DO 51 k = 1,3
271  rij2 = rij2 + (x(k,i) - x(k,j))**2
272  51 CONTINUE
273  pot2 = pot2 + 1.0/sqrt(rij2)
274  52 CONTINUE
275  54 CONTINUE
276  rp2 = float(np)*(np - 1)/(2.0*pot2)
277 *
278 * Sort square distances with respect to the centre C.
279  CALL sort1(np,r2,jlist)
280
281  DO 60 il = 1,lx
282  zm = 0.0
283  zmh = flagr(il)*zm2
284  i = 0
285  55 i = i + 1
286  im = jlist(i)
287  zm = zm + body(im)
288  IF (zm.LT.zmh) go to 55
289  rlagr(il) = sqrt(r2(i))
290  60 CONTINUE
291 *
292 * Obtain half-mass radius separately.
293  zm = 0.0
294  zmh = 0.5*zm2
295  i = 0
296  65 i = i + 1
297  im = jlist(i)
298  zm = zm + body(im)
299  IF (zm.LT.zmh) go to 65
300 *
301  rh2 = sqrt(r2(i))
302  np2 = np
303 *
304  WRITE (32,30) time+toff, (log10(rlagr(k)),k=1,lx)
305  CALL flush(32)
306 *
307  WRITE (6,70) time+toff, np0, np1, np2, rh1, rh2, rp0, rp1, rp2
308  70 FORMAT(/,' MASS GROUPS: T NP0 NP1 NP2 RM1 RM2 RP0 RP1 RP2 ',
309  & f7.1,i5,i6,i7,1x,5f7.3)
310 *
311  IF (n14.GT.1.OR.n13.GT.1) THEN
312  n14 = max(n14,1)
313  WRITE (6,80) time+toff, tphys, n14, zmbh/float(n14), rbh,
314  & semi, n13, r13, rx
315  80 FORMAT (/,' BH/NS SUBSYSTEM: T TPHYS NBH <MBH> RBH A ',
316  & 'NS R13 RX ',
317  & 2f7.1,i4,f6.1,f7.3,1p,e10.2,0p,i5,2f7.3)
318  END IF
319 *
320  RETURN
321 *
322  END