Nbody6
imf.f
Go to the documentation of this file.
1  SUBROUTINE imf(BODY10,BODYN)
2 *
3 *
4 * Initial mass function.
5 * ----------------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Generate Scalo IMF (Eggleton, Fitchett & Tout, Ap.J. 347, 998).
11  iter = 1
12 * Set upper limit of mass spectrum & corresponding mass fraction.
13  bodym = body10
14  xm = 0.998
15 * Find initial value of XM for upper & lower mass by iteration.
16  1 y0 = (1.0 - xm)**0.75 + 0.032*(1.0 - xm)**0.25
17  y2 = (0.75 + 0.008/sqrt(1.0 - xm))/
18  & ((1.0 - xm) + 0.032*sqrt(1.0 - xm))
19  y1 = bodym - 0.19*xm/y0
20  ypr = -0.19*(1.0 + xm*y2)/y0
21 * Try next guess of Newton-Raphson iteration.
22  xm = xm - y1/ypr
23 * Check upper & lower limits.
24  IF (xm.GT.1.0) xm = 0.99999
25  IF (xm.LT.0.0) xm = 0.001
26 * Set current mass and check the convergence.
27  bodys = 0.19*xm/y0
28  IF (abs(bodys - bodym).GT.1.0d-06*bodym) go to 1
29 *
30 * Save upper value of XM and perform second iteration.
31  IF (iter.EQ.1) THEN
32  x1 = xm
33  bodym = bodyn
34  xm = 0.4
35  iter = 2
36  go to 1
37  END IF
38 *
40  dx = (x1 - xm)/float(n - 1)
41  zmass = 0.0d0
42  DO 10 i = 1,n
43  xi = x1 - dx*float(i-1)
44  zm0 = (1.0 - xi)**0.75 + 0.032*(1.0 - xi)**0.25
45  body(i) = 0.19*xi/zm0
46  zmass = zmass + body(i)
47  10 CONTINUE
48 *
49  WRITE (6,20) body(1), body(n), zmass/float(n)
50  20 FORMAT (/,12x,'REALISTIC MASS FUNCTION:',' BODY(1) =',1pe9.2,
51  & ' BODY(N) =',e9.2,' <M> =',e9.2)
52 *
53 * Replace input value by actual mean mass in solar units.
54  zmbar = zmass/float(n)
55 *
56  RETURN
57 *
58  END