imf.f

Go to the documentation of this file.

      SUBROUTINE imf(BODY10,BODYN)
*
*
*       Initial mass function.
*       ----------------------
*
      include 'common6.h'
*
*
*       Generate Scalo IMF (Eggleton, Fitchett & Tout, Ap.J. 347, 998).
      iter = 1
*       Set upper limit of mass spectrum & corresponding mass fraction.
      bodym = body10
      xm = 0.998
*       Find initial value of XM for upper & lower mass by iteration.
    1 y0 = (1.0 - xm)**0.75 + 0.032*(1.0 - xm)**0.25
      y2 = (0.75 + 0.008/sqrt(1.0 - xm))/
     &                              ((1.0 - xm) + 0.032*sqrt(1.0 - xm))
      y1 = bodym - 0.19*xm/y0
      ypr = -0.19*(1.0 + xm*y2)/y0
*       Try next guess of Newton-Raphson iteration.
      xm = xm - y1/ypr
*       Check upper & lower limits.
      IF (xm.GT.1.0) xm = 0.99999
      IF (xm.LT.0.0) xm = 0.001
*       Set current mass and check the convergence.
      bodys = 0.19*xm/y0
      IF (abs(bodys - bodym).GT.1.0d-06*bodym) go to 1
*
*       Save upper value of XM and perform second iteration.
      IF (iter.EQ.1) THEN
          x1 = xm
          bodym = bodyn
          xm = 0.4
          iter = 2
          go to 1
      END IF
*
*       Set incremental fraction and assign individual masses.
      dx = (x1 - xm)/float(n - 1)
      zmass = 0.0d0
      DO 10 i = 1,n
          xi = x1 - dx*float(i-1)
          zm0 = (1.0 - xi)**0.75 + 0.032*(1.0 - xi)**0.25
          body(i) = 0.19*xi/zm0
          zmass = zmass + body(i)
   10 CONTINUE
*
      WRITE (6,20)  body(1), body(n), zmass/float(n)
   20 FORMAT (/,12x,'REALISTIC MASS FUNCTION:','   BODY(1) =',1pe9.2,
     &                                '  BODY(N) =',e9.2,'  <M> =',e9.2)
*
*       Replace input value by actual mean mass in solar units.
      zmbar = zmass/float(n)
*
      RETURN
*
      END