Nbody6
setup.f
Go to the documentation of this file.
1  SUBROUTINE setup
2 *
3 *
4 * Generation of initial coordinates & velocities.
5 * -----------------------------------------------
6 *
7  include 'common6.h'
8  parameter(a0=1.0d0)
9  REAL*8 ran2,a(8)
10 *
11 *
12 * Choose between uniform density and Plummer model.
13  kdum = idum1
14  IF (kz(5).GT.0) go to 20
15 *
16 * Set up a uniform spherical system.
17  DO 10 i = 1,n
18  1 a(1) = 0.0d0
19  DO 2 k = 1,3
20  a(k+1) = 2.0*ran2(kdum) - 1.0
21  a(1) = a(1) + a(k+1)**2
22  2 CONTINUE
23  IF (a(1).GT.1.0) go to 1
24 *
25  4 a(5) = 0.0d0
26  DO 5 k = 1,3
27  a(k+5) = 2.0*ran2(kdum) - 1.0
28  a(5) = a(5) + a(k+5)**2
29  5 CONTINUE
30  IF (a(5).GT.1.0) go to 4
31 *
32  DO 8 k = 1,3
33 * X(K,I) = A(1)*A(K+1)
34 * Density proportional to 1/R**2.
35  x(k,i) = a(k+1)
36 * Constant density.
37  xdot(k,i) = a(k+5)
38 * Isotropic velocities (magnitude randomized; no radial dependence).
39  8 CONTINUE
40  10 CONTINUE
41 *
42  go to200
43 *
44 * Initialize centre of mass terms.
45  20 DO 25 k = 1,3
46  cmr(k) = 0.0d0
47  cmrdot(k) = 0.0d0
48  25 CONTINUE
49 *
50 * Treat Plummer and lower options before Jaffe and Zhao models.
51  IF (kz(5).GE.5) go to 52
52 *
53 * Generate initial conditions from Plummer model (A & A 37, 183).
54  DO 40 i = 1,n
55  30 a(1) = ran2(kdum)
56  IF (a(1).LT.1.0d-10) go to 30
57  ri = (a(1)**(-0.6666667) - 1.0)**(-0.5)
58 * Reject distant particles.
59  IF (ri.GT.10.0) go to 30
60 *
61  a(2) = ran2(kdum)
62  a(3) = ran2(kdum)
63  x(3,i) = (1.0 - 2.0*a(2))*ri
64  x(1,i) = sqrt(ri**2 - x(3,i)**2)*cos(twopi*a(3))
65  x(2,i) = sqrt(ri**2 - x(3,i)**2)*sin(twopi*a(3))
66  32 a(4) = ran2(kdum)
67  a(5) = ran2(kdum)
68  a(6) = a(4)**2*(1.0 - a(4)**2)**3.5
69  IF (0.1*a(5).GT.a(6)) go to 32
70 *
71  a(8) = a(4)*sqrt(2.0)/(1.0 + ri**2)**0.25
72  a(6) = ran2(kdum)
73  a(7) = ran2(kdum)
74  xdot(3,i) = (1.0 - 2.0*a(6))*a(8)
75  xdot(1,i) = sqrt(a(8)**2 - xdot(3,i)**2)*cos(twopi*a(7))
76  xdot(2,i) = sqrt(a(8)**2 - xdot(3,i)**2)*sin(twopi*a(7))
77 *
78 * Accumulate c.m. terms.
79  DO 35 k = 1,3
80  cmr(k) = cmr(k) + body(i)*x(k,i)
81  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
82  35 CONTINUE
83  40 CONTINUE
84 *
85 * Scale coordinates & velocities to analytical expectation values.
86  sx = 1.5*twopi/16.0
87  sv = sqrt(zmass/sx)
88  42 DO 50 i = 1,n
89  DO 45 k = 1,3
90  x(k,i) = x(k,i) - cmr(k)/zmass
91  xdot(k,i) = xdot(k,i) - cmrdot(k)/zmass
92  x(k,i) = sx*x(k,i)
93  xdot(k,i) = sv*xdot(k,i)
94  45 CONTINUE
95  IF (kz(22).EQ.1) THEN
96  WRITE (10,48) body(i), (x(k,i),k=1,3), (xdot(k,i),k=1,3)
97  48 FORMAT (1p,7e14.6)
98  END IF
99  50 CONTINUE
100  IF (kz(22).EQ.1) CALL flush(10)
101 *
102 * Terminate standard cases (#5 <=1) or Jaffee/BH (#5 >= 5 label 52).
103  IF (kz(5).LE.1.OR.kz(5).GE.5) go to 200
104 *
105 * Save membership of first system for colour plot (N2 = NZERO - N1).
106  52 n1 = n
107 * Check optional generation of two orbiting Plummer models.
108  IF (kz(5).EQ.2) THEN
109  READ (5,*) apo, ecc, n2, scale
110  n2 = min(n,n2)
111  semi = apo/(1.0 + ecc)
112  semi = min(semi,50.0d0)
113  semi = max(semi,2.0d0)
114  ecc = min(ecc,0.999d0)
115  ecc = max(ecc,0.0d0)
116  zm2 = 0.0
117  kskip = n1/n2
118  DO 54 i = 1,n2
119  j = kskip*i
120  zm2 = zm2 + body(j)
121  54 CONTINUE
122  fac1 = zm2/(zmass + zm2)
123  fac2 = zmass/(zmass + zm2)
124 * Restrict volume ratio to 125 (i.e. unreasonable density contrast).
125  IF (scale.LE.0.2d0) scale = 0.2
126 * Increase total mass.
127  zmass = zmass + zm2
128 * Set apocentre velocity for new combined mass.
129  vap = sqrt(zmass/semi)*sqrt((1.0 - ecc)/(1.0 + ecc))
130  DO 55 i = 1,n
131  IF (i.LE.n2) THEN
132 * Copy members from first system by uniform skipping (N2 <= N1).
133  j = kskip*i
134  body(i+n) = body(j)
135  x(1,i+n) = scale*x(1,j) + fac2*apo
136  x(2,i+n) = scale*x(2,j)
137  x(3,i+n) = scale*x(3,j)
138  xdot(1,i+n) = xdot(1,j)/sqrt(scale)
139  xdot(2,i+n) = xdot(2,j)/sqrt(scale) + fac2*vap
140  xdot(3,i+n) = xdot(3,j)/sqrt(scale)
141  END IF
142  x(1,i) = x(1,i) - fac1*apo
143  xdot(2,i) = xdot(2,i) - fac1*vap
144  55 CONTINUE
145  ELSE IF (kz(5).EQ.3) THEN
146 * Prepare case of accretion disk with massive perturber.
147  READ (5,*) apo, ecc, dmin, scale
148  rin = 0.5
149  rout = 1.0
150  zmass = 1.0
151  body(1) = zmass
152  DO 58 k = 1,3
153  x(k,1) = 0.0
154  xdot(k,1) = 0.0
155  58 CONTINUE
156 * Generate a thin disk population in circular orbits.
157  DO 60 i = 2,n
158  body(i) = 1.0d-03/float(n)
159  semi = rin + (rout - rin)*float(i)/float(n)
160  vcirc = sqrt((body(1) + body(i))/semi)
161  phase = twopi*ran2(kdum)
162  x(1,i) = semi*cos(phase)
163  x(2,i) = semi*sin(phase)
164  x(3,i) = 0.01*(2.0*ran2(kdum) - 1.0)
165  xdot(1,i) = -vcirc*sin(phase)
166  xdot(2,i) = vcirc*cos(phase)
167  xdot(3,i) = 0.01*vcirc*(2.0*ran2(kdum) - 1.0)
168  60 CONTINUE
169 * Define membership of perturber and ensure no external tide.
170  n2 = 1
171  kz(14) = 0
172 * Redefine solar mass unit and astronomical length scale in AU.
173  zmbar = 1.0
174  rbar = 1.0/2.05d+05
175  body(n+1) = scale*body(1)
176 * Set appropriate mass ratios for transforming to new c.m. frame.
177  fac1 = body(n+1)/(zmass + body(n+1))
178  fac2 = zmass/(zmass + body(n+1))
179  zmass = zmass + body(n+1)
180 * Form orbital elements for massive perturber (avoid ECC = 1).
181  IF (abs(ecc - 1.0).GT.1.0d-05) THEN
182  semi = dmin/(1.0 - ecc)
183  ELSE
184  semi = -1.0d+05
185  END IF
186  vm2 = zmass*(2.0/dmin - 1.0/semi)
187  vap2 = zmass*(2.0/apo - 1.0/semi)
188 * Determine initial y-velocity from angular momentum conservation.
189  vy = sqrt(vm2)*dmin/apo
190  vx = sqrt(vap2 - vy**2)
191 * Place perturber on the Y-axis with appropriate velocities.
192  x(1,n+1) = apo*fac2
193  x(2,n+1) = 0.0
194  x(3,n+1) = 0.0
195  xdot(1,n+1) = -vx*fac2
196  xdot(2,n+1) = vy*fac2
197  xdot(3,n+1) = 0.0
198 * Displace the disk members and include negative y-velocity.
199  DO 70 i = 1,n
200  x(1,i) = x(1,i) - fac1*apo
201  xdot(1,i) = xdot(1,i) + fac1*vx
202  xdot(2,i) = xdot(2,i) - fac1*vy
203  70 CONTINUE
204  ELSE IF (kz(5).EQ.4) THEN
205 * Include two massive bodies (ECC > 1: NAME = 1 & 2 free floating).
206  n2 = 0
207  READ (5,*) semi, ecc, zm1, zm2
208  WRITE (6,72) semi, ecc, zm1, zm2
209  72 FORMAT (/,12x,'MASSIVE BODIES A =',1p,e9.1,
210  & ' E =',0p,f6.2,' M1/<M> =',f6.2,' M2/<M> =',f6.2)
211  body(1) = zm1
212  body(2) = zm2
213  IF (ecc.LT.1.0) THEN
214 * Set apocentre velocity for new combined mass (using NAME = 1 & 2).
215  vap = sqrt((zm1 + zm2)/semi)*sqrt((1.0 - ecc)/(1.0 + ecc))
216  fac1 = zm2/(zm1 + zm2)
217  fac2 = zm1/(zm1 + zm2)
218  DO 75 k = 1,3
219  x(k,1) = 0.0
220  x(k,2) = 0.0
221  xdot(k,2) = 0.0
222  75 CONTINUE
223 * Initialize binary with c.m. at rest (elements change in SCALE).
224  x(1,1) = -fac1*semi*(1.0 + ecc)
225  x(1,2) = fac2*semi*(1.0 + ecc)
226  xdot(2,1) = -fac1*vap
227  xdot(2,2) = fac2*vap
228  END IF
229  ELSE IF (kz(5).EQ.5) THEN
230 * Generate Jaffe model with length scale A0.
231  zmh = 1.0/sqrt(float(2*n))
232  DO 100 i = 1,n
233  80 zm = ran2(kdum)
234  ri = zm*a0/(1.0 - zm)
235 * Reject distant particles.
236  IF (ri.GT.10.0.OR.ri.LT.5.0d-03) go to 80
237  a(2) = ran2(kdum)
238  a(3) = ran2(kdum)
239  x(3,i) = (1.0 - 2.0*a(2))*ri
240  x(1,i) = sqrt(ri**2 - x(3,i)**2)*cos(twopi*a(3))
241  x(2,i) = sqrt(ri**2 - x(3,i)**2)*sin(twopi*a(3))
242  vi2 = 1.0/(ri + a0) + 2.0*zmh/ri
243  DO 90 k = 1,3
244  vk2 = vi2*ran2(kdum)/3.0
245  zz = 2.0*ran2(kdum) - 1.0
246  ss = 1.0
247  IF (zz.LT.0.0) ss = -1.0
248  xdot(k,i) = ss*sqrt(vk2)
249  90 CONTINUE
250 * Accumulate c.m. terms.
251  DO 95 k = 1,3
252  cmr(k) = cmr(k) + body(i)*x(k,i)
253  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
254  95 CONTINUE
255  100 CONTINUE
256  sx = 1.0
257  sv = 1.0
258 * Use existing procedure for c.m. correction.
259  go to 42
260  ELSE IF (kz(5).GE.6) THEN
261 * Generate dwarf galaxy model with cusp & BH (Zhao, MN 278, 488).
262  y0 = 5.0*1.25/twopi/sqrt(float(n))
263 * Calculate RH and ZMH by setting mass inside RH to 5*N^{1/2}*<M>.
264  rh = (exp(y0) - 1.0d0)**0.4
265  zmh = 4.0*twopi/5.0*log(1.0d0 + rh**2.5)
266 * Include possibility of using actual mass (read by SCALE).
267  rcut = 0.0
268  IF (kz(24).LT.0) THEN
270  IF (zmh.GT.0.0) THEN
271  y0 = 1.25/twopi*(zmh/float(n))
272  rh = (exp(y0) - 1.0d0)**0.4
273  END IF
274  END IF
275 *
276 * Add experimental choice of favourite value.
277  IF (kz(28).EQ.3) THEN
278  zmh = 1.0/sqrt(float(2*n))
279  y0 = 1.25/twopi*(zmh/float(n))
280  rh = (exp(y0) - 1.0d0)**0.4
281  END IF
282 *
283 * Define BH mass inside radius of influence (unscaled units).
284  reject = 100.0
285  cutm = 4.0*twopi/5.0*log(1.0d0 + reject**2.5)
286 * Specify mass fraction using cut-off value and total in M_sun.
287  zmh = zmh/(float(n)*zmbar)*cutm
288  y0 = 1.25/twopi*zmh
289  rh = (exp(y0) - 1.0d0)**0.4
290  r2 = 100.0
291  jdum = -1
292  ranj = gasdev(jdum)
293  IF (kz(5).EQ.7) THEN
294  rh = 0.006
295  zmh = log(1.0d0 + rh**0.75)
296  reject = 200.0
297  cutm = log(1.0d0 + reject**0.75)
298  r2 = 200.0
299  END IF
300 *
301  icut = 0
302  DO 120 i = 1,n
303  110 zm = ran2(kdum)
304  IF (kz(5).EQ.6) THEN
305  zx = 1.25/twopi*zm
306  zx = zx*cutm
307  zy = exp(zx)
308  ri = (zy - 1.0d0)**0.4
309  ELSE
310  zx = zm*cutm
311  zy = exp(zx)
312  ri = (zy - 1.0d0)**(4.0/3.0)
313  END IF
314 * Reject large and small distances before scaling.
315  IF (ri.GT.reject.OR.ri.LT.1.0d-04) go to 110
316  a(2) = ran2(kdum)
317  a(3) = ran2(kdum)
318  x(3,i) = (1.0 - 2.0*a(2))*ri
319  x(1,i) = sqrt(ri**2 - x(3,i)**2)*cos(twopi*a(3))
320  x(2,i) = sqrt(ri**2 - x(3,i)**2)*sin(twopi*a(3))
321  115 IF (kz(5).EQ.6) THEN
322  CALL vdisp(ri,r2,zmh,vi2)
323  ELSE
324  CALL vdisp2(ri,r2,zmh,vi2)
325  END IF
326  sigma = sqrt(vi2)
327  vx = sigma*gasdev(jdum)
328  vy = sigma*gasdev(jdum)
329  vz = sigma*gasdev(jdum)
330  IF (max(abs(vx),abs(vy),abs(vz)).GT.3.0*sigma) go to 115
331  xdot(1,i) = vx
332  xdot(2,i) = vy
333  xdot(3,i) = vz
334  IF (ri.GT.rcut) go to 120
335  semi = 2.0/ri - (vx**2 + vy**2 + vz**2)/zmh
336  semi = 1.0/semi
337  rrdot = x(1,i)*vx + x(2,i)*vy + x(3,i)*vz
338  ecc2 = (1.0 - ri/semi)**2 + rrdot**2/(semi*zmh)
339  ecc = sqrt(ecc2)
340  IF (semi*(1.0 - ecc).LT.0.5*rcut) THEN
341  icut = icut + 1
342  go to 110
343  END IF
344  120 CONTINUE
345 *
346  WRITE (6,125) rh, zmh, cutm, rcut, icut
347  125 FORMAT (/,12x,'BLACK HOLE RH =',f6.3,' MBH =',f8.4,
348  & ' CUTM =',f7.3,' RCUT =',f9.5,' ICUT =',i5)
349 *
350  sx = 1.0
351  sv = 1.0
352 * Substitute c.m. value (hardly matters with zero velocity).
353  IF (kz(24).LT.0) body(1) = zmh
354  DO 130 k = 1,3
355  x(k,1) = 0.0
356  xdot(k,1) = 0.0
357  cmr(k) = 0.0
358  cmrdot(k) = 0.0
359  130 CONTINUE
360  DO 140 i = 1,n
361  DO 135 k = 1,3
362  cmr(k) = cmr(k) + body(i)*x(k,i)
363  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
364  135 CONTINUE
365  140 CONTINUE
366  go to 42
367  END IF
368 *
369 * Specify new membership.
370  n = n + n2
371  nzero = n
372  ntot = n
373  IF (n.GE.nmax-10) THEN
374  WRITE (6,150) n, nmax
375  150 FORMAT (' DANGER! LIMIT EXCEEDED N =',i6,' NMAX =',i6)
376  stop
377  END IF
378 *
379  IF (kz(5).EQ.2) THEN
380  WRITE (6,155) semi, ecc, n1, n2, scale
381  155 FORMAT (/,12x,'PLUMMER BINARY A =',f6.2,' E =',f6.2,
382  & ' N1 =',i6,' N2 =',i6,' SCALE =',f6.2)
383  ELSE IF (kz(5).EQ.3) THEN
384  WRITE (6,160) apo, ecc, dmin, scale
385  160 FORMAT (/,12x,'MASSIVE PERTURBER APO =',f6.2,' E =',f6.2,
386  & ' DMIN =',f6.2,' MP/M1 =',f6.2)
387  END IF
388 *
389 * Re-initialize centre of mass terms.
390  DO 162 k = 1,3
391  cmr(k) = 0.0d0
392  cmrdot(k) = 0.0d0
393  162 CONTINUE
394  zmass = 0.0
395  DO 170 i = 1,n
396  zmass = zmass + body(i)
397  DO 165 k = 1,3
398  cmr(k) = cmr(k) + body(i)*x(k,i)
399  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
400  165 CONTINUE
401  170 CONTINUE
402  DO 180 i = 1,n
403  DO 175 k = 1,3
404  x(k,i) = x(k,i) - cmr(k)/zmass
405  xdot(k,i) = xdot(k,i) - cmrdot(k)/zmass
406  175 CONTINUE
407  180 CONTINUE
408 *
409 * Save random number sequence in COMMON for future use.
410  200 idum1 = kdum
411 *
412  RETURN
413 *
414  END
415
416  REAL*8 FUNCTION gasdev(IDUM)
417 C
418 C GAUSSIAN DISTRIBUTION (PRESS P. 203).
419 C -------------------------------------
420 C
421  common/randg/ gset,iset
422 CC DATA ISET/0/
423  REAL*8 ran2
424  IF (iset.EQ.0) THEN
425 1 v1=2.*ran2(idum)-1.
426  v2=2.*ran2(idum)-1.
427  r=v1**2+v2**2
428  IF(r.GE.1.)go to 1
429  fac=sqrt(-2.*log(r)/r)
430  gset=v1*fac
431  gasdev=v2*fac
432  iset=1
433  ELSE
434  gasdev=gset
435  iset=0
436  ENDIF
437  RETURN
438  END
439
440  SUBROUTINE vdisp(RI,R2,ZMH,VI2)
441 *
442 * Velocity dispersion of dwarf galaxy model.
443 * ------------------------------------------
444 *
445  IMPLICIT REAL*8 (a-h,o-z)
446 *
447 *
448  one3 = 1.0/3.0d0
449  sum = 0.0
450  rho = 1.0/(sqrt(ri)*(1.0 + ri**2*sqrt(ri)))
451  ncall = 1000 + 100.0/ri**2
452  dr = 2.0*(r2 - ri)/float(ncall)
453  u = ri
454 * Solve by Simpson's rule (NB! DR is twice the interval).
455  DO 30 j = 1,ncall
456  u1 = u
457  CALL vfunc(zmh,u1,f1)
458  u2 = u + 0.5d0*dr
459  CALL vfunc(zmh,u2,f2)
460  u3 = u + dr
461  CALL vfunc(zmh,u3,f3)
462  sum = sum + one3*(f1 + 4.0d0*f2 + f3)*dr/2.0d0
463  u = u + dr
464  IF (u.LT.1.0) THEN
465  dr = 1.4*dr
466  END IF
467  IF (u.GT.100.0) go to 40
468  30 CONTINUE
469 *
470  40 vi2 = sum/rho
471 *
472  RETURN
473  END
474
475  SUBROUTINE vfunc(ZMH,U,F)
476  IMPLICIT REAL*8 (a-h,m,o-z)
477 *
478  u2 = u**2
479  zmu = 0.8*6.2831853*log(1.0d0 + u2*sqrt(u))
480  rho = 1.0/(sqrt(u)*(1.0 + u2*sqrt(u)))
481  f = rho*(zmu + zmh)/u2
482 *
483  RETURN
484  END
485
486  SUBROUTINE vdisp2(RI,R2,ZMH,VI2)
487 *
488 * Velocity dispersion of r^{-9/4} galaxy model.
489 * ---------------------------------------------
490 *
491  IMPLICIT REAL*8 (a-h,o-z)
492 *
493 *
494  one3 = 1.0/3.0d0
495  sum = 0.0
496  rho = 1.0/(ri**2.25*(1.0 + ri**0.75))
497  ncall = 1000 + 100.0/ri**2
498  dr = 2.0*(r2 - ri)/float(ncall)
499  u = ri
500 * Solve by Simpson's rule (NB! DR is twice the interval).
501  DO 30 j = 1,ncall
502  u1 = u
503  CALL vfunc2(zmh,u1,f1)
504  u2 = u + 0.5d0*dr
505  CALL vfunc2(zmh,u2,f2)
506  u3 = u + dr
507  CALL vfunc2(zmh,u3,f3)
508  sum = sum + one3*(f1 + 4.0d0*f2 + f3)*dr/2.0d0
509  u = u + dr
510  IF (u.LT.1.0) THEN
511  dr = 1.4*dr
512  END IF
513  IF (u.GT.200.0) go to 40
514  30 CONTINUE
515 *
516  40 vi2 = sum/rho
517 *
518  RETURN
519  END
520
521  SUBROUTINE vfunc2(ZMH,U,F)
522  IMPLICIT REAL*8 (a-h,m,o-z)
523 *
524  u2 = u**2
525  zmu = log(1.0d0 + u**0.75)
526  rho = 1.0/(u**2.25*(1.0 + u**0.75))
527  f = rho*(zmu + zmh)/u2
528 *
529  RETURN
530  END