
******************************************************************
* GETMASS - given a maximum & minimum mass range, randomly
* select a star mass such that a large ensemble of the output
* will follow the Kroupa (2001; MNRAS 322, 231) IMF.
* 
* INPUT:
* masmax = Maximum mass probed for random masses
* masmin = Minimum "    "      "   "      "
* OUTPUT:
* mass = output mass
* Eric Mamajek 11/9/01
******************************************************************
      program getmass

      implicit double precision (a-z)
      integer idum,choice,n,no,nb,na,nf,ng,nk,nm

      open(unit = 10,file="imf.out")
      open(unit = 20,file="imf_O.out")
      open(unit = 21,file="imf_B.out")
      open(unit = 22,file="imf_A.out")
      open(unit = 23,file="imf_F.out")
      open(unit = 24,file="imf_G.out")
      open(unit = 25,file="imf_K.out")
      open(unit = 26,file="imf_M.out")
      open(unit = 27,file="imf_substellar.out")

* INITIALIZE
      n  = 0
      no = 0
      nb = 0
      na = 0
      nf = 0
      ng = 0
      nk = 0
      nm = 0
      sum   = 0.00000000d0
      sum2  = 0.0000000d0

* ZAMS mass distribution
      momax = 150.0d0   ! O stars
      momin =  13.2d0   ! OB star boundary mass
      mbmin =   2.4d0   ! BA star boundary mass
      mamin =   1.55d0  ! AF star boundary mass
      mfmin =   1.12d0  ! FG star boundary mass
      mgmin =   0.82d0  ! G s
      mkmin =   0.52d0  ! K stars
      mmmin =   0.08d0  ! M stars

      mbmax = momin
      mamax = mbmin
      mfmax = mamin
      mgmax = mfmin 
      mkmax = mgmin
      mmmax = mkmin

* Power law indices from Kroupa 2001, where xi ~ m^-a
      a0 = 0.3    ! 0.01 < m/Mo < 0.08
      a1 = 1.3    ! 0.08 < m/Mo < 0.50
      a2 = 2.3    ! 0.50 < m/Mo 

      m0 = 0.01
      m1 = 0.08   ! divisions between power laws
      m2 = 0.50

      norm0 = 0.25118864  ! = (m0)**a0
      norm1 = 0.02009509  ! = m1*norm0
      norm2 = 0.01004755  ! = m2*norm1

      write(*,*)'Create fake stellar sample using Kroupa(2001) IMF'
      write(*,*)'[0] create # stars for N solar masses w/ IMF'
      write(*,*)'[else] create THIS number of stars from IMF...'
      read(*,*)choice

      write(*,*)'Recommended MIN mass: 0.08 Msun'
      write(*,*)'Recommended MAX mass: 150  Msun'
      write(*,*)'MIN mass, MAX mass [Msun], negative integer seed'
      read(*,*)masmin,masmax,idum

*      masmin = 0.08d0  ! defaults
*      masmax = 150.0d0 ! defaults
*      idum   = -7      ! defaults

      range = masmax - masmin
      done  = 0

      if (choice .eq. 0) then 
        write(*,*)'Solar masses of stars?'
        read(*,*)sum 
        do while (sum .gt. 0) 
          x = masmin + range*ran3(idum) ! PICK RANDOM MASS IN RANGE
          if (x .ge. m2) then !HIGH MASS END >1Mo
            y = norm2*x**(-1*a2) !SALPETER LAW
          elseif (x .ge. m1 .and. x .lt. m2) then
            y = norm1*x**(-1*a1)
          else                !BELOW M1 = 0.08 Mo
            y = norm0*x**(-1*a2)
          endif
          yran = rann2(idum)
          if (yran .lt. y) then
            mass = x
            write(10,999)mass

            if (mass .gt. momin .and. mass .le. momax) then
               write(20,999)mass
               no = no + 1
            elseif (mass .gt. mbmin .and. mass .le. mbmax) then
               write(21,999)mass
               nb = nb + 1
            elseif (mass .gt. mamin .and. mass .le. mamax) then
               write(22,999)mass
               na = na + 1
            elseif (mass .gt. mfmin .and. mass .le. mfmax) then
               write(23,999)mass
               nf = nf + 1
            elseif (mass .gt. mgmin .and. mass .le. mgmax) then
               write(24,999)mass
               ng = ng + 1
            elseif (mass .gt. mkmin .and. mass .le. mkmax) then
               write(25,999)mass
               nk = nk + 1
            elseif (mass .gt. mmmin .and. mass .le. mmmax) then
               write(26,999)mass
               nm = nm + 1
            else
               write(27,999)mass
            endif

            sum  = sum  - mass
            sum2 = sum2 + mass
            n = n + 1
          endif
        enddo
        meanmass = sum2/dble(n)
        write(*,*)' N = ',n,' stars '
        write(*,*)' M = ',sum2,' M/Mo'
        write(*,*)'<M>= ',meanmass, ' M/Mo '
      else ! just make the number of stars that I entered... and quit asking questions
        nmax = choice 
        i = 0
        do 80 while (i .le. nmax) 
          i = i + 1
          do 50 while (done .ne. 1) 
            x = masmin + range*ran3(idum) ! PICK RANDOM MASS IN RANGE
            if (x .ge. m2) then !HIGH MASS END >1Mo
               y = norm2*x**(-1*a2) !SALPETER LAW
            elseif (x .ge. m1 .and. x .lt. m2) then
               y = norm1*x**(-1*a1)
            else                !BELOW M1 = 0.08 Mo
               y = norm0*x**(-1*a2)
            endif
            yran = rann2(idum)
            if (yran .lt. y) then
               done = 1
               mass = x
               write(10,999)mass

               if (mass .gt. momin .and. mass .le. momax) then
                  write(20,999)mass
                  no = no + 1
               elseif (mass .gt. mbmin .and. mass .le. mbmax) then
                  write(21,999)mass
                  nb = nb + 1
               elseif (mass .gt. mamin .and. mass .le. mamax) then
                  write(22,999)mass
                  na = na + 1
               elseif (mass .gt. mfmin .and. mass .le. mfmax) then
                  write(23,999)mass
                  nf = nf + 1
               elseif (mass .gt. mgmin .and. mass .le. mgmax) then
                  write(24,999)mass
                  ng = ng + 1
               elseif (mass .gt. mkmin .and. mass .le. mkmax) then
                  write(25,999)mass
                  nk = nk + 1
               elseif (mass .gt. mmmin .and. mass .le. mmmax) then
                  write(26,999)mass
                  nm = nm + 1
               else
                  write(27,999)mass
               endif
               sum = sum + mass
               n   = n + 1
            endif
 50       end do 
          done = 0
 80     end do
        meanmass = sum/dble(n)
        write(*,*)' N  = ',n,' stars '
        write(*,*)' M  = ',sum,' M/Mo'
        write(*,*)'<M> = ',meanmass, ' M/Mo '
        write(*,*)'N(O) = ',no,' +- ',dsqrt(dble(no))
        write(*,*)'N(B) = ',nb,' +- ',dsqrt(dble(nb))
        write(*,*)'N(A) = ',na,' +- ',dsqrt(dble(na))
        write(*,*)'N(F) = ',nf,' +- ',dsqrt(dble(nf))
        write(*,*)'N(G) = ',ng,' +- ',dsqrt(dble(ng))
        write(*,*)'N(K) = ',nk,' +- ',dsqrt(dble(nk))
        write(*,*)'N(M) = ',nm,' +- ',dsqrt(dble(nm))
      endif
      write(*,*)'see output files: imf*.out'
      write(*,*)'type >wc imf*.out'
 999  format(f7.3)
      close(10)
      close(20)
      close(21)
      close(22)
      close(23)
      close(24)
      close(25)
      close(26)
      close(27)
      end



******************************************************************
* RANN2 - random number function
******************************************************************
      function rann2(idum)
      integer idum,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
      real rann2,am,eps,rnmx
      parameter (im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1,
     *  ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211,
     *  ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7,rnmx=1.-eps)
      integer idum2,j,k,iv(ntab),iy
      save iv,iy,idum2
      data idum2/123456789/, iv/ntab*0/, iy/0/
      if (idum .le. 0) then
         idum = max(-idum,1)
         idum2 = idum
         do j = ntab+8,1,-1
            k = idum/iq1
            idum = ia1*(idum-k*iq1) - k*ir1
            if (idum .lt. 0) idum = idum + im1
            if (j .le. ntab) iv(j) = idum
         enddo
         iy = iv(1)
      endif
      k = idum/iq1
      idum = ia1*(idum-k*iq1) - k*ir1
      if (idum .lt. 0) idum = idum + im1
      k = idum2/iq2
      idum2 = ia2*(idum2-k*iq2) - k*ir2
      if (idum2 .lt. 0) idum2 = idum2 + im2
      j = 1 + iy/ndiv
      iy = iv(j) - idum2
      iv(j) = idum
      if (iy .lt. 1) iy = iy + imm1
      rann2 = min(am*iy,rnmx)
      return
      end

******************************************************************
* RAN3 - random number function
******************************************************************

      FUNCTION ran3(idum)
      INTEGER idum
      INTEGER MBIG,MSEED,MZ
C     REAL MBIG,MSEED,MZ
      REAL ran3,FAC
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)
C     REAL mj,mk,ma(55)
      SAVE iff,inext,inextp,ma
      DATA iff /0/
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum)
        mj=mod(mj,MBIG)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3=mj*FAC
      return
      END
