Nbody6
 All Files Functions Variables
imf2.f
Go to the documentation of this file.
1  SUBROUTINE imf2(BODY10,BODYN)
2 *
3 *
4 * Mass function for binaries & single stars.
5 * ------------------------------------------
6 * BODY10,BODYN are in M_sun.
7 *
8  include 'common6.h'
9  INTEGER kcm(nmax)
10  REAL*8 ran2,imfbd
11  REAL*8 lm,um,bcm
12  EXTERNAL imfbd
13  common/work1/ bcm(nmax)
14  REAL*8 copyx(3,nmax),copyv(3,nmax)
15  parameter(maxm=0)
16 *
17 *========================= if KZ(20)=2:
18 * KTG3 MF with alpha1=1.3 :
19  DATA g1,g2,g3,g4 /0.19,1.55,0.050,0.6/
20 * KTG3 MF with alpha1=1.1 :
21 c DATA G1,G2,G3,G4 /0.28,1.14,0.010,0.1/
22 *=========================
23 *
24 *
25 * Change PARAMETER to MAXM = 1 for taking BODY10 as massive star.
26  bdymax = 0.0d0
27 *
28 * Generate initial mass function (N-NBIN0 singles & 2*NBIN0 binaries).
29  kdum = idum1
30  zmass = 0.0d0
31  DO 10 i = 1,n+nbin0
32  5 xx = ran2(kdum)
33 *
34 * Choose between Kroupa et al. (M.N. 262, 545) & Eggleton (Book).
35  IF (kz(20).EQ.2.OR.kz(20).EQ.4) THEN
36  zm = 0.08 + (g1*xx**g2 + g3*xx**g4)/(1.0 - xx)**0.58
37  ELSE IF (kz(20).EQ.3.OR.kz(20).EQ.5) THEN
38  zm = 0.3*xx/(1.0 - xx)**0.55
39 * Allow discrimination between correlated & uncorrelated binary masses.
40  ELSE IF (kz(20).EQ.6.OR.kz(20).EQ.7) THEN
41  lm = bodyn
42  um = body10
43  zm = imfbd(xx,lm,um)
44  END IF
45 * Include possibility of setting non-MS types.
46  kstar(i) = 1
47 *
48 * See whether the mass falls within the specified range.
49  IF (zm.GE.bodyn.AND.zm.LE.body10) THEN
50  body(i) = zm
51  zmass = zmass + body(i)
52  ELSE
53  go to 5
54  END IF
55 *
56 * Find index of most massive star (C.O. 20/12/10).
57  IF (body(i).GT.bdymax) THEN
58  bdymax = body(i)
59  imaxx = i
60  END IF
61  10 CONTINUE
62 *
63 * Set maximum mass from upper limit BODY10 and correct total mass.
64  IF (maxm.GT.0) THEN
65  zmass = zmass - body(imaxx)
66  body(imaxx) = body10
67  zmass = zmass + body(imaxx)
68  END IF
69 *
70 * See whether to skip mass function for binaries.
71  IF (nbin0.EQ.0) go to 50
72 *
73 * Merge binary components in temporary variable for sorting.
74  DO 20 i = 1,nbin0
75  bcm(i) = body(2*i-1) + body(2*i)
76  jlist(i) = i
77  kcm(2*i-1) = kstar(2*i-1)
78  kcm(2*i) = kstar(2*i)
79  20 CONTINUE
80 *
81 * Sort total binary masses in increasing order.
82  IF (nbin0.GT.1) THEN
83  CALL sort1(nbin0,bcm,jlist)
84  END IF
85 *
86 * Save scaled binary masses in decreasing order for routine BINPOP.
87  DO 30 i = 1,nbin0
88  jb = jlist(nbin0-i+1)
89  body0(2*i-1) = max(body(2*jb-1),body(2*jb))/zmass
90  body0(2*i) = min(body(2*jb-1),body(2*jb))/zmass
91 *( C.O. 15.11.07)
92 * Adopt correlation f(m1/m2) also for Brown Dwarf IMF (kz(20) = 7).
93 * Set up realistic primordial binary population according to
94 * a) Kouwenhoven et al. 2007, A&A, 474, 77
95 * b) Kobulnicky & Fryer 2007, ApJ, 670, 747
96  IF ((kz(20).eq.4).OR.
97  & (kz(20).eq.5).OR.
98  & (kz(20).eq.7)) THEN
99 * Adopt correlation (m1/m2)' = (m1/m2)**0.4 & constant sum (Eggleton).
100  zmb = body0(2*i-1) + body0(2*i)
101  ratio = body0(2*i-1)/body0(2*i)
102  body0(2*i) = zmb/(1.0 + ratio**0.4)
103  body0(2*i-1) = zmb - body0(2*i)
104  END IF
105  kstar(2*i-1) = kcm(2*jb-1)
106  kstar(2*i) = kcm(2*jb)
107  30 CONTINUE
108 *
109 * Merge binary components into single stars for scaling purposes.
110  zmb = 0.0
111  DO 40 i = 1,nbin0
112  body(nbin0-i+1) = bcm(i)
113  zmb = zmb + bcm(i)
114  40 CONTINUE
115 *
116  WRITE (6,45) nbin0, body(1), body(nbin0), zmb/float(nbin0)
117  45 FORMAT (//,12x,'BINARY STAR IMF: NB =',i6,
118  & ' RANGE =',1p,2e10.2,' <MB> =',e9.2)
119 *
120 * Move the single stars up to form compact array of N members.
121  50 IF (n.LE.nbin0) go to 90
122  ns = 0
123  DO 60 l = 1,n-nbin0
124  body(nbin0+l) = body(2*nbin0+l)
125  ns = ns + 1
126  bcm(ns) = body(nbin0+l)
127  kcm(ns) = kstar(2*nbin0+l)
128 *TEST> C.O. 20.12.2010
129 * Create local copies of phase space vector.
130  DO 55 k = 1,3
131  copyx(k,ns) = x(k,nbin0 + l)
132  copyv(k,ns) = xdot(k,nbin0 + l)
133  55 CONTINUE
134  jlist(ns) = ns
135  60 CONTINUE
136 *
137 * Sort masses of single stars in increasing order.
138  CALL sort1(ns,bcm,jlist)
139 *
140 * Copy masses of single stars to COMMON in decreasing order.
141  zms = 0.0
142  DO 70 i = 1,ns
143  body(n-i+1) = bcm(i)
144  zms = zms + bcm(i)
145  kstar(n-i+1+nbin0) = kcm(jlist(i))
146 *TEST> C.O. 20.12.2010
147 * Recover the corresponding coordinates and velocities.
148  j = jlist(i)
149  DO 65 k = 1,3
150  x(k,n-i+1) = copyx(k,j)
151  xdot(k,n-i+1) = copyv(k,j)
152  65 CONTINUE
153  70 CONTINUE
154 *
155  WRITE (6,80) n-nbin0, body(nbin0+1), body(n), zms/float(n-nbin0)
156  80 FORMAT (/,12x,'SINGLE STAR IMF: NS =',i6,' RANGE =',1p,2e10.2,
157  & ' <MS> =',e9.2)
158 *
159 * Replace input value by actual mean mass in solar units.
160  90 zmbar = zmass/float(n)
161 *
162 * Save random number sequence in IDUM1 for future use.
163  idum1 = kdum
164 *
165  RETURN
166 *
167  END