Nbody6
 All Files Functions Variables
imfbd.f
Go to the documentation of this file.
1 C
2 C Note: this routine contains the correction from Carsten Weidner
3 C (Nov.2004), that for LOM above a certain value the integral
4 C was not computed correctly.
5 C
6 C
7 C===========================================================
8  REAL*8 FUNCTION imfbd(XX,LOM,UPM)
9 *
10 *
11 * 5-part-power-law IMF extending into BD regime.
12 * ---------------------------------------------
13 * (P.Kroupa, Feb.99, Heidelberg)
14 *
15 * On input of a random deviate XX and the lower (LM) and upper (UM)
16 * mass limits (in Msun), IMFBD returns a mass in Msun.
17 * (For equal masses, LOM=UPM and IMFBD = UPM.)
18 * LOM and UPM must be kept unchanged.
19 *
20 *
21  IMPLICIT NONE
22  REAL*8 ml,mh,m0,m1,m2,mu
23  REAL*8 alpha0,alpha1,alpha2,alpha3,alpha4
24  REAL*8 k,xh,x0,x1,x2,xu,mup
25  REAL*8 mtot,massh,mass0,mass1,mass2,massu
26  REAL*8 xx,lom,upm,lm,um,zm
27  REAL*8 beta,intgr,mintgr,mgen,xin
28  EXTERNAL intgr,mintgr,mgen
29  LOGICAL firstbd
30 *
31  SAVE firstbd,xh,x0,x1,x2,xu,k,lm,um
32  DATA firstbd /.true./
33 *
34 * -------------------------------------------------------------------
35 * MF (Kroupa, Tout & Gilmore, 1993, MNRAS, 262, 545)
36 * extended into BD regime (masses in Msun).
37 * (Note: ALPHA0 is between ML and MH, etc.)
38 *
39 c DATA ML, MH, M0, M1, M2, MU/
40 c & 0.01D0, 0.08D0, 0.5D0, 1.0D0, 8.0D0, 100.0D0/
41 c DATA ALPHA0, ALPHA1, ALPHA2, ALPHA3, ALPHA4/
42 c & 0.3D0, 1.3D0, 2.2D0, 2.7D0, 2.7D0/
43 *
44 *
45 C Standard of average IMF:
46 *
47  DATA ml, mh, m0, m1, m2, mu/
48  & 0.01d0, 0.08d0, 0.5d0, 1.0d0, 8.0d0, 100.0d0/
49  DATA alpha0, alpha1, alpha2, alpha3, alpha4/
50  & 0.3d0, 1.3d0, 2.3d0, 2.3d0, 2.3d0/
51 *
52 *
53 *
54 C DATA ML, MH, M0, M1, M2, MU/
55 C & 0.01D0, 0.08D0, 0.5D0, 1.0D0, 8.0D0, 150.0D0/
56 C DATA ALPHA0, ALPHA1, ALPHA2, ALPHA3, ALPHA4/
57 C & 0.3D0, 1.3D0, 2.3D0, 2.35D0, 2.35D0/
58 *
59 *
60 *.....................................................................
61 * Salpeter MF:
62 *
63 C DATA ML, MH, M0, M1, M2, MU/
64 C & 0.01D0, 0.08D0, 0.08D0, 1.0D0, 8.0D0, 100.0D0/
65 C DATA ALPHA0, ALPHA1, ALPHA2, ALPHA3, ALPHA4/
66 C & 0.3D0, 0.3D0, 1.2D0, 2.3D0, 2.3D0/
67 *
68 *.....................................................................
69 * Flat MF:
70 *
71 c DATA ML, MH, M0, M1, M2, MU/
72 c & 0.01D0, 0.08D0, 0.5D0, 1.0D0, 5.0D0, 100.0D0/
73 c DATA ALPHA0, ALPHA1, ALPHA2, ALPHA3, ALPHA4/
74 c & 1.D0, 1.D0, 1.D0, 1.D0, 1.D0/
75 * --------------------------------------------------------------------
76 *
77  IF (.NOT.firstbd.AND.um.EQ.lm)THEN
78  zm = um
79  goto 500
80  END IF
81 *
82  IF (firstbd) THEN
83  lm = lom
84  um = upm
85 *
86 * Remind user of encoded MF:
87 *
88  WRITE(6,*)
89  WRITE(6,'(a)')' The 5-part power-law MF:'
90  WRITE(6,'(a)')' ML MH M0 M1 M2 MU'
91  WRITE(6,'(6F8.3)')ml,mh,m0,m1,m2,mu
92  WRITE(6,'(a)')' ALPHA0 ALPHA1 ALPHA2 ALPHA3 ALPHA4'
93  WRITE(6,'(3x,6F8.3)')alpha0,alpha1,alpha2,alpha3,alpha4
94  WRITE(6,*)
95 *
96 * Check MF:
97 *
98  IF (ml.GT.mh) THEN
99  WRITE(6,*)
100  WRITE(6,*)' FATAL: ML > MH'
101  WRITE(6,*)' STOP'
102  WRITE(6,*)
103  stop
104  END IF
105  IF (mh.GT.m0) THEN
106  WRITE(6,*)
107  WRITE(6,*)' FATAL: MH > M0'
108  WRITE(6,*)' STOP'
109  WRITE(6,*)
110  stop
111  END IF
112  IF (m0.GT.m1) THEN
113  WRITE(6,*)
114  WRITE(6,*)' FATAL: M0 > M1'
115  WRITE(6,*)' STOP'
116  WRITE(6,*)
117  stop
118  END IF
119  IF (m1.GT.m2) THEN
120  WRITE(6,*)
121  WRITE(6,*)' FATAL: M1 > M2'
122  WRITE(6,*)' STOP'
123  WRITE(6,*)
124  stop
125  END IF
126  IF (m2.GT.mu) THEN
127  WRITE(6,*)
128  WRITE(6,*)' FATAL: M2 > MU'
129  WRITE(6,*)' STOP'
130  WRITE(6,*)
131  stop
132  END IF
133 *
134 * Check for physical input mass range.
135 *
136  IF (um.EQ.lm) THEN
137  WRITE(6,*)
138  WRITE(6,*)' WARNING: UM (= ',um,' Msun) = LM= ',
139  & lm,' Msun'
140  WRITE(6,*)
141  firstbd = .false.
142  zm = um
143  goto 500
144  END IF
145 *
146  IF (um.LT.lm) THEN
147  WRITE(6,*)
148  WRITE(6,*)' FATAL: UM (= ',um,' Msun) < LM= ',
149  & lm,' Msun'
150  WRITE(6,*)' STOP'
151  WRITE(6,*)
152  stop
153  END IF
154  IF (lm.GT.mu) THEN
155  WRITE(6,*)
156  WRITE(6,*)' FATAL: BODYN (LM= ',lm,' Msun) > MU'
157  WRITE(6,*)' STOP'
158  WRITE(6,*)
159  stop
160  END IF
161  IF (lm.LT.ml) THEN
162  WRITE(6,*)
163  WRITE(6,*)' WARNING: BODYN (LM= ',lm,' Msun) < ML'
164  WRITE(6,*)' BODYN = ',ml, ' Msun adopted.'
165  WRITE(6,*)
166  lm = ml
167  END IF
168  IF (um.GT.mu) THEN
169  WRITE(6,*)
170  WRITE(6,*)' WARNING: BODY1 (UM= ',um,' Msun) > MU= ',
171  & mu,' Msun'
172  WRITE(6,*)' BODY1 = ',mu, ' Msun adopted.'
173  WRITE(6,*)
174  um = mu
175  END IF
176 *
177 * Determine normalisation K for unit probablity (LM<=M<=UM).
178 *
179  k = 0.0d0
180  xh = 0.0d0
181  x0 = 0.0d0
182  x1 = 0.0d0
183  x2 = 0.0d0
184  xu = 0.0d0
185 *
186  IF (lm.GE.ml.AND.lm.LE.mh) THEN
187 *
188  IF (um.LE.mh) THEN
189  mup = um
190  ELSE IF (um.GT.mh) THEN
191  mup = mh
192  END IF
193  beta = mh**alpha0
194  xh = intgr(lm,mup,alpha0,beta)
195 *
196  IF (um.GT.mh) THEN
197  IF (um.LE.m0) THEN
198  mup = um
199  ELSE IF (um.GT.m0) THEN
200  mup = m0
201  END IF
202  beta = mh**alpha1
203  x0 = intgr(mh,mup,alpha1,beta)
204  END IF
205 *
206  IF (um.GT.m0) THEN
207  IF (um.LE.m1) THEN
208  mup = um
209  ELSE IF (um.GT.m1) THEN
210  mup = m1
211  END IF
212  beta = m0**alpha2 * (mh/m0)**alpha1
213  x1 = intgr(m0,mup,alpha2,beta)
214  END IF
215 *
216  IF (um.GT.m1) THEN
217  IF (um.LE.m2) THEN
218  mup = um
219  ELSE IF (um.GT.m2) THEN
220  mup = m2
221  END IF
222  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
223  x2 = intgr(m1,mup,alpha3,beta)
224  END IF
225 *
226  IF (um.GT.m2) THEN
227  mup = um
228  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
229  + (m1/m2)**alpha3
230  xu = intgr(m2,mup,alpha4,beta)
231  END IF
232 *
233  ELSE IF (lm.GT.mh.AND.lm.LE.m0) THEN
234 *
235  IF (um.LE.m0) THEN
236  mup = um
237  ELSE IF (um.GT.m0) THEN
238  mup = m0
239  END IF
240  beta = mh**alpha1
241  x0 = intgr(lm,mup,alpha1,beta)
242 *
243  IF (um.GT.m0) THEN
244  IF (um.LE.m1) THEN
245  mup = um
246  ELSE IF (um.GT.m1) THEN
247  mup = m1
248  END IF
249  beta = m0**alpha2 * (mh/m0)**alpha1
250  x1 = intgr(m0,mup,alpha2,beta)
251  END IF
252 *
253  IF (um.GT.m1) THEN
254  IF (um.LE.m2) THEN
255  mup = um
256  ELSE IF (um.GT.m2) THEN
257  mup = m2
258  END IF
259  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
260  x2 = intgr(m1,mup,alpha3,beta)
261  END IF
262 *
263  IF (um.GT.m2) THEN
264  mup = um
265  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
266  + (m1/m2)**alpha3
267  xu = intgr(m2,mup,alpha4,beta)
268  END IF
269 *
270  ELSE IF (lm.GT.m0.AND.lm.LE.m1) THEN
271 *
272  IF (um.LE.m1) THEN
273  mup = um
274  ELSE IF (um.GT.m1) THEN
275  mup = m1
276  END IF
277  beta = m0**alpha2 * (mh/m0)**alpha1
278  x1 = intgr(lm,mup,alpha2,beta)
279 *
280  IF (um.GT.m1) THEN
281  IF (um.LE.m2) THEN
282  mup = um
283  ELSE IF (um.GT.m2) THEN
284  mup = m2
285  END IF
286  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
287  x2 = intgr(m1,mup,alpha3,beta)
288  END IF
289 *
290  IF (um.GT.m2) THEN
291  mup = um
292  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
293  + (m1/m2)**alpha3
294  xu = intgr(m2,mup,alpha4,beta)
295  END IF
296 *
297  ELSE IF (lm.GT.m1.AND.lm.LE.m2) THEN
298 *
299  IF (um.LE.m2) THEN
300  mup = um
301  ELSE IF (um.GT.m2) THEN
302  mup = m2
303  END IF
304  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
305  x2 = intgr(lm,mup,alpha3,beta)
306 *
307  IF (um.GT.m2) THEN
308  mup = um
309  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
310  + (m1/m2)**alpha3
311  xu = intgr(m2,mup,alpha4,beta)
312  END IF
313 *
314  ELSE IF (lm.GT.m2) THEN
315 *
316  mup = um
317  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
318  + (m1/m2)**alpha3
319  xu = intgr(lm,mup,alpha4,beta)
320 *
321  END IF
322 *
323 *
324  k = xh + x0 + x1 + x2 + xu
325  xh = xh/k
326  x0 = x0/k
327  x1 = x1/k
328  x2 = x2/k
329  xu = xu/k
330 *
331 * Write number fractions to log file.
332 *
333  WRITE(6,'(a,F7.3,a,F7.3,a)')
334  & ' Number fractions for LM=',lm,', Msun, UM=',um,' Msun:'
335  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') ml,' -',mh,' Msun:', xh
336  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') mh,' -',m0,' Msun:', x0
337  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') m0,' -',m1,' Msun:', x1
338  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') m1,' -',m2,' Msun:', x2
339  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') m2,' -',mu,' Msun:', xu
340  WRITE(6,*)
341 *
342 * Determine average mass and mass fractions (for log-file).
343 *
344  mtot = 0.0d0
345  massh = 0.0d0
346  mass0 = 0.0d0
347  mass1 = 0.0d0
348  mass2 = 0.0d0
349  massu = 0.0d0
350 *
351  IF (lm.GE.ml.AND.lm.LE.mh) THEN
352 *
353  IF (um.LE.mh) THEN
354  mup = um
355  ELSE IF (um.GT.mh) THEN
356  mup = mh
357  END IF
358  beta = mh**alpha0
359  massh = mintgr(lm,mup,alpha0,beta)
360 *
361  IF (um.GT.mh) THEN
362  IF (um.LE.m0) THEN
363  mup = um
364  ELSE IF (um.GT.m0) THEN
365  mup = m0
366  END IF
367  beta = mh**alpha1
368  mass0 = mintgr(mh,mup,alpha1,beta)
369  END IF
370 *
371  IF (um.GT.m0) THEN
372  IF (um.LE.m1) THEN
373  mup = um
374  ELSE IF (um.GT.m1) THEN
375  mup = m1
376  END IF
377  beta = m0**alpha2 * (mh/m0)**alpha1
378  mass1 = mintgr(m0,mup,alpha2,beta)
379  END IF
380 *
381  IF (um.GT.m1) THEN
382  IF (um.LE.m2) THEN
383  mup = um
384  ELSE IF (um.GT.m2) THEN
385  mup = m2
386  END IF
387  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
388  mass2 = mintgr(m1,mup,alpha3,beta)
389  END IF
390 *
391  IF (um.GT.m2) THEN
392  mup = um
393  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
394  + (m1/m2)**alpha3
395  massu = mintgr(m2,mup,alpha4,beta)
396  END IF
397 *
398  ELSE IF (lm.GT.mh.AND.lm.LE.m0) THEN
399 *
400  IF (um.LE.m0) THEN
401  mup = um
402  ELSE IF (um.GT.m0) THEN
403  mup = m0
404  END IF
405  beta = mh**alpha1
406  mass0 = mintgr(lm,mup,alpha1,beta)
407 *
408  IF (um.GT.m0) THEN
409  IF (um.LE.m1) THEN
410  mup = um
411  ELSE IF (um.GT.m1) THEN
412  mup = m1
413  END IF
414  beta = m0**alpha2 * (mh/m0)**alpha1
415  mass1 = mintgr(m0,mup,alpha2,beta)
416  END IF
417 *
418  IF (um.GT.m1) THEN
419  IF (um.LE.m2) THEN
420  mup = um
421  ELSE IF (um.GT.m2) THEN
422  mup = m2
423  END IF
424  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
425  mass2 = mintgr(m1,mup,alpha3,beta)
426  END IF
427 *
428  IF (um.GT.m2) THEN
429  mup = um
430  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
431  + (m1/m2)**alpha3
432  massu = mintgr(m2,mup,alpha4,beta)
433  END IF
434 *
435  ELSE IF (lm.GT.m0.AND.lm.LE.m1) THEN
436 *
437  IF (um.LE.m1) THEN
438  mup = um
439  ELSE IF (um.GT.m1) THEN
440  mup = m1
441  END IF
442  beta = m0**alpha2 * (mh/m0)**alpha1
443  mass1 = mintgr(lm,mup,alpha2,beta)
444 *
445  IF (um.GT.m1) THEN
446  IF (um.LE.m2) THEN
447  mup = um
448  ELSE IF (um.GT.m2) THEN
449  mup = m2
450  END IF
451  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
452  mass2 = mintgr(m1,mup,alpha3,beta)
453  END IF
454 *
455  IF (um.GT.m2) THEN
456  mup = um
457  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
458  + (m1/m2)**alpha3
459  massu = mintgr(m2,mup,alpha4,beta)
460  END IF
461 *
462  ELSE IF (lm.GT.m1.AND.lm.LE.m2) THEN
463 *
464  IF (um.LE.m2) THEN
465  mup = um
466  ELSE IF (um.GT.m2) THEN
467  mup = m2
468  END IF
469  beta = m1**alpha3 * (mh/m0)**alpha1 * (m0/m1)**alpha2
470  mass2 = mintgr(lm,mup,alpha3,beta)
471 *
472  IF (um.GT.m2) THEN
473  mup = um
474  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
475  + (m1/m2)**alpha3
476  massu = mintgr(m2,mup,alpha4,beta)
477  END IF
478 *
479  ELSE IF (lm.GT.m2) THEN
480 *
481  mup = um
482  beta = m2**alpha4 * (mh/m0)**alpha1 * (m0/m1)**alpha2 *
483  + (m1/m2)**alpha3
484  massu = mintgr(lm,mup,alpha4,beta)
485 *
486  END IF
487 *
488 *
489  mtot = massh + mass0 + mass1 + mass2 + massu
490  massh = massh/mtot
491  mass0 = mass0/mtot
492  mass1 = mass1/mtot
493  mass2 = mass2/mtot
494  massu = massu/mtot
495 *
496 * Write mass fractions to log file.
497 *
498  WRITE(6,'(a,F7.3,a,F7.3,a)')
499  & ' Mass fractions for LM=',lm,', Msun, UM=',um,' Msun:'
500  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') ml,' -',mh,' Msun:', massh
501  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') mh,' -',m0,' Msun:', mass0
502  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') m0,' -',m1,' Msun:', mass1
503  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') m1,' -',m2,' Msun:', mass2
504  WRITE(6,'(F7.3,a,F7.3,a,F7.4)') m2,' -',mu,' Msun:', massu
505  WRITE(6,'(a,F7.3,a)')' Expected average stellar mass: ',
506  & mtot/k,' Msun'
507  WRITE(6,*)
508 *
509  firstbd = .false.
510  END IF
511 *
512 * Generate stellar mass from random deviate.
513 *
514  IF (lm.GE.ml.AND.lm.LE.mh) THEN
515 *
516  IF (xx.LE.xh) THEN
517  xin = xx
518  beta = k/mh**alpha0
519  zm = mgen(xin,lm,alpha0,beta)
520  ELSE IF (xx.GT.xh.AND.xx.LE.(xh+x0)) THEN
521  xin = xx-xh
522  beta = k/mh**alpha1
523  zm = mgen(xin,mh,alpha1,beta)
524  ELSE IF (xx.GT.(xh+x0).AND.xx.LE.(xh+x0+x1)) THEN
525  xin = xx-xh-x0
526  beta = k/m0**alpha2 * (m0/mh)**alpha1
527  zm = mgen(xin,m0,alpha2,beta)
528  ELSE IF (xx.GT.(xh+x0+x1).AND.xx.LE.(xh+x0+x1+x2)) THEN
529  xin = xx-xh-x0-x1
530  beta = k/m1**alpha3 * (m1/m0)**alpha2 * (m0/mh)**alpha1
531  zm = mgen(xin,m1,alpha3,beta)
532  ELSE IF (xx.GT.(xh+x0+x1+x2)) THEN
533  xin = xx-xh-x0-x1-x2
534  beta = k/m2**alpha4 * (m2/m1)**alpha3 *
535  + (m1/m0)**alpha2 * (m0/mh)**alpha1
536  zm = mgen(xin,m2,alpha4,beta)
537  END IF
538 *
539  ELSE IF (lm.GT.mh.AND.lm.LE.m0) THEN
540 *
541  IF (xx.GT.xh.AND.xx.LE.(xh+x0)) THEN
542  xin = xx-xh
543  beta = k/mh**alpha1
544  zm = mgen(xin,lm,alpha1,beta)
545  ELSE IF (xx.GT.(xh+x0).AND.xx.LE.(xh+x0+x1)) THEN
546  xin = xx-xh-x0
547  beta = k/m0**alpha2 * (m0/mh)**alpha1
548  zm = mgen(xin,m0,alpha2,beta)
549  ELSE IF (xx.GT.(xh+x0+x1).AND.xx.LE.(xh+x0+x1+x2)) THEN
550  xin = xx-xh-x0-x1
551  beta = k/m1**alpha3 * (m1/m0)**alpha2 * (m0/mh)**alpha1
552  zm = mgen(xin,m1,alpha3,beta)
553  ELSE IF (xx.GT.(xh+x0+x1+x2)) THEN
554  xin = xx-xh-x0-x1-x2
555  beta = k/m2**alpha4 * (m2/m1)**alpha3 *
556  + (m1/m0)**alpha2 * (m0/mh)**alpha1
557  zm = mgen(xin,m2,alpha4,beta)
558  END IF
559 *
560  ELSE IF (lm.GT.m0.AND.lm.LE.m1) THEN
561 *
562  IF (xx.GT.(xh+x0).AND.xx.LE.(xh+x0+x1)) THEN
563  xin = xx-xh-x0
564  beta = k/m0**alpha2 * (m0/mh)**alpha1
565  zm = mgen(xin,lm,alpha2,beta)
566  ELSE IF (xx.GT.(xh+x0+x1).AND.xx.LE.(xh+x0+x1+x2)) THEN
567  xin = xx-xh-x0-x1
568  beta = k/m1**alpha3 * (m1/m0)**alpha2 * (m0/mh)**alpha1
569  zm = mgen(xin,m1,alpha3,beta)
570  ELSE IF (xx.GT.(xh+x0+x1+x2)) THEN
571  xin = xx-xh-x0-x1-x2
572  beta = k/m2**alpha4 * (m2/m1)**alpha3 *
573  + (m1/m0)**alpha2 * (m0/mh)**alpha1
574  zm = mgen(xin,m2,alpha4,beta)
575  END IF
576 *
577  ELSE IF (lm.GT.m1.AND.lm.LE.m2) THEN
578 *
579  IF (xx.GT.(xh+x0+x1).AND.xx.LE.(xh+x0+x1+x2)) THEN
580  xin = xx-xh-x0-x1
581  beta = k/m1**alpha3 * (m1/m0)**alpha2 * (m0/mh)**alpha1
582  zm = mgen(xin,lm,alpha3,beta)
583  ELSE IF (xx.GT.(xh+x0+x1+x2)) THEN
584  xin = xx-xh-x0-x1-x2
585  beta = k/m2**alpha4 * (m2/m1)**alpha3 *
586  + (m1/m0)**alpha2 * (m0/mh)**alpha1
587  zm = mgen(xin,m2,alpha4,beta)
588  END IF
589 *
590  ELSE IF (lm.GT.m2) THEN
591 *
592  IF (xx.GT.(xh+x0+x1+x2)) THEN
593  xin = xx-xh-x0-x1-x2
594  beta = k/m2**alpha4 * (m2/m1)**alpha3 *
595  + (m1/m0)**alpha2 * (m0/mh)**alpha1
596  zm = mgen(xin,lm,alpha4,beta)
597  END IF
598 *
599  END IF
600 *
601  500 imfbd = zm
602 *
603  RETURN
604  END
605 C===========================================================
606  real*8 FUNCTION intgr(MA,MB,ALPHA,BETA)
607 *
608 * Integral of mass function over mass interval MA to MB
609 * -----------------------------------------------------
610 *
611  IMPLICIT NONE
612  REAL*8 alpha,beta,ma,mb
613 *
614 *
615  IF (alpha.NE.1.d0) THEN
616  intgr = (mb**(1.d0-alpha)-ma**(1.d0-alpha)) * beta/(1.d0-alpha)
617  ELSE IF (alpha.EQ.1.d0) THEN
618  intgr = log(mb/ma) * beta
619  END IF
620 *
621  RETURN
622  END
623 C===========================================================
624  real*8 FUNCTION mintgr(MA,MB,ALPHA,BETA)
625 *
626 * Total mass over mass interval MA to MB
627 * --------------------------------------
628 *
629  IMPLICIT NONE
630  REAL*8 alpha,beta,ma,mb
631 *
632 *
633  IF (alpha.NE.2.d0) THEN
634  mintgr = (mb**(2.d0-alpha)-ma**(2.d0-alpha)) * beta/(2.d0-alpha)
635  ELSE IF (alpha.EQ.2.d0) THEN
636  mintgr = log(mb/ma) * beta
637  END IF
638 *
639  RETURN
640  END
641 C===========================================================
642  real*8 FUNCTION mgen(XIN,MA,ALPHA,BETA)
643 *
644 * Generate mass
645 * -------------
646 *
647  IMPLICIT NONE
648  REAL*8 alpha,beta,xin,ma
649 *
650 *
651  IF (alpha.NE.1.d0) THEN
652  mgen = (1.d0-alpha) * xin * beta
653  mgen = (mgen + ma**(1.d0-alpha))**(1.d0/(1.d0-alpha))
654  ELSE IF (alpha.EQ.1.d0) THEN
655  mgen = ma * dexp(xin * beta)
656  END IF
657 *
658  RETURN
659  END
660 C===========================================================