program massfn
      implicit double precision (a-h,o-z)
      omega=0.25d0
      omegal=0.7d0
      omegab=0.05d0
c      write(6,*) 'h,sigma_8,n'
c      read (*,*) h,sig8,tn
      h=0.7d0
      sig8=0.8d0
      tn=1.d0
      call ResetPowerSpectrum(omega,omegal,omegab,h,sig8,tn)
      z=-1.d0
CC     Number of grid points in redshift:
      jmax=21
CC     Number of grid points in mass (enter twice, as integer in `imax' 
CC                                   and double-precision in `fimax'):
      imax=4000 
      fimax=4000.d0
CC     Minimum mass considered (in solar masses):
      tMmin=1.d6
CC     Maximum mass considered:
      tMmax=1.d15
CC     Maximum redshift considered:
      zmax=20.d0
c
      deltaz=(zmax+1.d0)/jmax
      ftM=(tMmax/tMmin)**(1.d0/fimax)
      do j=1,jmax
      z=z+deltaz
      tM=tMmin
      do i=1,imax
      tM=tM*ftM
      x=1.326d0/(1.d0+z)
      sigM=sigmaOfM(tM,dsdM,1)*GrowthFac(z)
      sigminx=1.047d0+0.599d0*(datan(7.386d0*(x-0.526d0))/
     &    3.14d0+0.5d0)
      sigmint=1.047d0+0.599d0*(datan(7.386d0*(1.393d0-0.526d0))/
     &    3.14d0+0.5d0)
      bone=sigminx/sigmint
      cminx=3.681d0+1.352d0*(datan(6.948d0*(x-0.424d0))/
     &    3.14d0+0.5d0)
      cmint=3.681d0+1.352d0*(datan(6.948d0*(1.393d0-0.424d0))/
     &    3.14d0+0.5d0)
      bzero=cminx/cmint
      sigp=sigM*bone
      csigp=2.881d0*((sigp/1.257d0)**1.022d0+1.d0)*
     &    dexp(0.06d0/(sigp*sigp))
      c=bzero*csigp
         write(12,'(d20.6,d20.6,d20.6)') z,dlog10(tM),c
      enddo
      enddo
      end

      function dndlM(z,tM)
c
c     Find dn/dlogM in 1/(Mpc^3), at mass tM (solar).
c     n is the comoving halo density.
c
      implicit double precision (a-h,o-z)
      double precision N_nu
      common/GLOBALVARIABLES/omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu,
     $     beta_c,sound_horizon,theta_cmb,omega,omegal,z_equality
      pi=3.1415926535897932385d0
      omz=fomega(omega,omegal,z)
      dcritz=del_crit0(omz,omegal)/GrowthFac(z)
      sigM=sigmaOfM(tM,dsdM,1)
      dlsdlM=tM*dsdM/sigM
      tdn=dsqrt(2.d0/pi)*dcritz*dabs(dlsdlM)*dexp(
     *     -0.707d0*dcritz*dcritz/(2.d0*sigM*sigM))/(tM*sigM)
c     *     -dcritz*dcritz/(2.d0*sigM*sigM))/(tM*sigM)
      tdn=tdn*2.78d11*omhh
c Sheth Tormen correction factor
      tdn=tdn*0.2709*(1.d0+(sigM*sigM/(0.707d0*dcritz*dcritz))**0.3d0)
      dndlM=tdn
      return
      end

      double precision function TFintegrator(f,a,b,tol)
c
c     This calls an integration routine, and integrates the function
c     f between a and b, to relative accuracy tol. The routine used here
c     is rombint, but this can easily be changed.
c
      implicit double precision (a-h,k,o-z)
      double precision f
      external f
      TFintegrator=rombint(f,a,b,tol)
      return
      end

      subroutine Growth(z,DD0)
c
c     From Eisenstein & Hu:
c     DD0 = D1(z)/D1(z=0) inverse of growth from z to 0     
c
      implicit double precision (a-h,k,o-z)
      double precision N_nu

      common/GLOBALVARIABLES/omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu,
     *     beta_c,sound_horizon,theta_cmb,omega,omegal,z_equality

c      q = k*theta_cmb**2/omhh
c
c      y_fs = 17.2*f_nu*(1.+0.488*f_nu**(-7./6.))*(N_nu*q/f_nu)**2

      oz=omega*(1.d0+z)**3 /(omegal+(1.d0-omegal-omega)*(1.d0+z)**2
     $     +omega*(1.d0+z)**3)
      olz=omegal /(omegal+(1.d0-omegal-omega)*(1.d0+z)**2+omega*(1.d0+z)
     $     **3)

      D = (1.d0+z_equality)/(1.d0+z)*5.d0*oz/2.d0*(oz**(4.d0/7.d0)-olz
     $     +(1.d0+oz/2.d0)*(1.d0+olz/70.d0))**(-1.d0)

      DD0= D/((1.d0+z_equality)*5.d0*omega/2.d0*(omega**(4.d0/7.d0)
     $     -omegal+(1.d0+omega/2.d0)*(1.d0+omegal/70.d0))**(-1.d0))

c
c     With neutrinos:
c
c      p_cb = -(5.-sqrt(1.+24*(1.-f_nu)))/4.
c      DD_cb = (1.+(D/(1.+y_fs))**(0.7))**(-p_cb/0.7)*D**(p_cb)
c
c      DD_cbnu = ((1.-f_nu)**(-0.7/p_cb)
c     *     +(D/(1.+y_fs))**(0.7))**(-p_cb/0.7)*D**(p_cb)

      return
      end 

      function GrowthFac(z)
      implicit double precision (a-h,o-z)
      double precision GrowthFac
      call Growth(z,t)
      GrowthFac=t
      return
      end

      function del_crit0(om0,omLam)
      implicit double precision (a-h,o-z)
      double precision omLam,del_crit0   
c
c     This is the very minor dependence on cosmology of
c     the critical collapse overdensity 1.69 .
c
      pi=3.1415926535897932385d0
      t1=.15d0*(12.d0*pi)**(2.d0/3.d0)
      if (om0 .gt. (1.d0-1.d-5)) then
         del_crit0=t1
      elseif (omLam .lt. 1.d-5) then
         del_crit0=t1*(om0**.0185d0)
      else
         del_crit0=t1*(om0**.0055d0)
      endif
      return
      end

      function fomega(om0,omLam,z)
      implicit double precision (a-h,o-z)
      double precision fomega
      f=omLam/(1.d0+z)**3
      fomega=om0/(om0+f+(1.d0-om0-omLam)/(1.d0+z))
      return
      end

      subroutine ResetPowerSpectrum(omegaIn,omegalIn,omegab,h,sig8,tn)
c
c     Sets transfer function parameters; also sets the normalization by
c     computing sigma_8.
c
      implicit double precision (a-h,k,o-z)
      double precision N_nu,sigmatop,sigmatop2
      common /GLOBALVARIABLES/ omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu,
     *     beta_c,sound_horizon,theta_cmb,omega,omegal,z_equality
      common /PASSSIGMA/ scale,tilt,zEHu,ipower
      common /normalize/ snorm
      external sigmatop,sigmatop2
      omega=omegaIn
      omegal=omegalIn
      omeganu=0.d0
      T_cmb=2.726d0
      N_nu=0.d0
      zEHU=0.d0
c     Translate Parameters into GLOBALVARIABLES form
      f_nu = omeganu/omega
      f_baryon = omegab/omega
      theta_cmb=T_cmb/2.7d0
      omhh = omega*h*h
      call TFset_parameters
      tilt = tn
c     No gravity waves:
      ipower=1
c     
c     Normalize to Sigma_8: Mass fluctuation in a top hat of 8./h Mpc at
c     redshift z=0; total mass: CDM+baryon
c     
      scale = 8.d0/h      
c      ipower=0
      tol= 1.d-7
      sigma_8= dsqrt(TFintegrator(sigmatop2,1.d-6,0.001d0/scale,tol)+
     $     TFintegrator(sigmatop,dlog(0.001d0/scale),dlog(0.1d0/scale)
     $     ,tol)+ TFintegrator(sigmatop,dlog(0.1d0/scale),dlog(1.d0
     $     /scale),tol)+ TFintegrator(sigmatop,dlog(1.d0/scale),dlog(10
     $     .d0/scale),tol)+ TFintegrator(sigmatop,dlog(10.d0/scale)
     $     ,dlog(100.d0/scale),tol))
      snorm=sig8/sigma_8
      return
      end

      double precision function sigmaOfM(tM,dsdM,ideriv)
c
c     tM is the halo mass in solar mass. Returns sigma(tM) and
c     calculates (only if ideriv=1) dsdM=d sigma(tM)/d tM
c
      implicit double precision (a-h,k,o-z)
      parameter (pi=3.1415926535897932385d0)
      double precision N_nu,klimits(7),sigmatop,sigmatop2,dsigmatop
     $     ,dsigmatop2
      common /GLOBALVARIABLES/ omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu,
     $     beta_c,sound_horizon,theta_cmb,omega,omegal,z_equality
      common /PASSSIGMA/ scale,tilt,zEHU,ipower
      common /normalize/ snorm
      external sigmatop,sigmatop2,dsigmatop,dsigmatop2
c     scale is in comoving Mpc:
      scale=(tM*3.d0/(4.d0*pi*2.77545d11*omhh))**(1.d0/3.d0)
      tol=2.d-6
      klimits(1)=dmin2(1.d-3/scale,.2d0)
      klim=2.d2/scale
      khi=klim
      klhi=dlog(khi)
      if (klhi .gt. klimits(1)) klim=dmin2(klim,khi)
      nk=int(dlog10(klim/klimits(1))+.5d0)
      if (nk .lt. 5) then 
         nk=5
      elseif (nk .gt. 7) then
         nk=7
      endif
      kfac=dlog(klim/klimits(1))/(nk-1.d0)
      klimits(1)=dlog(klimits(1))
      do i=2,nk
         klimits(i)=klimits(i-1)+kfac
      enddo
      klow=1.d-9
      if (khi .lt. klow) then
         sigmaOfM=0.d0
      elseif (klhi .lt. klimits(1)) then
         sigmaOfM=TFintegrator(sigmatop2,klow,khi,tol)
      else
         sigmaOfM=TFintegrator(sigmatop2,klow,dexp(klimits(1)),tol)
         do i=2,nk
            sigmaOfM=sigmaOfM+TFintegrator(sigmatop,klimits(i-1
     $           ),klimits(i),tol)
         enddo
      endif
      sigmaOfM=snorm*dsqrt(sigmaOfM)
      if (ideriv .eq. 1) then
         if (khi .lt. klow) then
            dsigmaOfM=0.d0
         elseif (klhi .lt. klimits(1)) then
            dsigmaOfM=TFintegrator(dsigmatop2,klow,khi,tol)
         else
            dsigmaOfM=TFintegrator(dsigmatop2,klow,dexp(klimits(1)),tol)
            do i=2,nk
               dsigmaOfM=dsigmaOfM+TFintegrator(dsigmatop,klimits(i-1)
     $              ,klimits(i),tol)
            enddo
         endif
c
c     Boundary term from integration by parts:
c
         dsigmaOfM=dsigmaOfM+sigmatop(klhi)/scale
         dsdM=1.d0/(8.d0*pi*2.77545d11*omhh* scale**2)
         dsdM=snorm**2 *dsdM*dsigmaOfM/sigmaOfM
      endif
      return
      end

      double precision function TF_master(k,dTFdk,ideriv)
c     
c     Uses Eisenstein & Hu online code: 
c     http://www.sns.ias.edu/~whu/transfer/transfer.html
c
c     Modification by Rennan Barkana:
c     Also calculates dTFdk=d TF_master/d k (if ideriv=1)
c
      implicit double precision (a-h,k,o-z)
      double precision N_nu
      common /GLOBALVARIABLES/ omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu,
     *     beta_c,sound_horizon,theta_cmb,omega,omegal,z_equality
c
      q = k*theta_cmb**2/omhh
      gamma_eff=dsqrt(alpha_nu) + (1.d0-dsqrt(alpha_nu))/
     *     (1.d0+(0.43d0*k*sound_horizon)**4)

      q_eff = q/gamma_eff
      TF_master= dlog(dexp(1.d0)+1.84d0*beta_c*dsqrt(alpha_nu)*q_eff)
      if (ideriv .eq. 1) then
         t1=TF_master
         dt1dqe=1.84d0*beta_c*dsqrt(alpha_nu)/dexp(TF_master)
      endif
      TF_master = TF_master/(TF_master + q_eff**2*
     *     (14.4d0 + 325.d0/(1.d0+60.5d0*q_eff**1.11d0)))
c
      if (ideriv .eq. 1) then
         t2=TF_master
         v=t1/t2
         dv=dt1dqe+2.d0*q_eff*(14.4d0 + 325.d0/(1.d0+60.5d0*q_eff**1
     $        .11d0))-325.d0*60.5d0*1.11d0 *q_eff**2.11d0/(1.d0+60.5d0
     $        *q_eff**1.11d0)**2
         dgamedk=-(1.d0-dsqrt(alpha_nu))*4.d0* k**3 *(0.43d0
     $        *sound_horizon)**4/(1.d0+(0.43d0*k*sound_horizon)**4)**2
         dt2dk=(theta_cmb**2/omhh)*((v*dt1dqe-t1*dv)/v**2) *(1.d0
     $        /gamma_eff-k*dgamedk/gamma_eff**2)
         dTFdk = dt2dk
      endif
      if (N_nu .gt. 1.d-10) pause 'Sorry, no neutrinos here.'
c      q_nu = 3.92d0*q/dsqrt(f_nu/N_nu)
c      TF_master = TF_master*
c     *     (1.d0+(1.2d0*f_nu**(0.64d0)*N_nu**(0.3d0+0.6d0*f_nu))/
c     *	   (q_nu**(-1.6d0)+q_nu**(0.8d0)))
      return 
      end
      
      double precision function dmin2(a,b)
      implicit double precision (a-h,k,o-z)
      if (b .lt. a) then
         dmin2=b
      else
         dmin2=a
      endif
      return
      end

      double precision function sigmatop2(k)
      implicit double precision (a-h,k,o-z)
      sigmatop2=sigmatop(dlog(k))/k
      return
      end

      double precision function sigmatop(kl)
c
c     Specialized to z=0!
c      
      implicit double precision (a-h,k,o-z)
      common /PASSSIGMA/ scale,tilt,zEHU,ipower
      
      k = dexp(kl)
      x = scale*k

      if (zEHU .gt. 1.d-10) pause 'Sorry, z=0 here.'
c     
c     For non-zero z:
c      call Growth(z,k,DD_cb,DD_cbnu,DD0)     
      DD0=1.d0
      DD_cbnu=1.d0
      DD_cb=1.d0

      if (ipower.eq.0) then
         sigmatop = k**(3.d0+tilt)*(DD0*DD_cbnu*TF_master(k,dummy,0))**2
     $        *(3.d0*(x*dcos(x) - dsin(x))/x**3)**2
      else

         sigmatop = k**(3.d0+tilt)*(DD0*DD_cb*TF_master(k,dummy,0))**2
     $        *(3.d0*(x*dcos(x) - dsin(x))/x**3)**2
      endif
      return
      end

      double precision function dsigmatop2(k)
      implicit double precision (a-h,k,o-z)
      dsigmatop2=dsigmatop(dlog(k))/k
      return
      end

      double precision function dsigmatop(kl)
      implicit double precision (a-h,k,o-z)
      common /PASSSIGMA/ scale,tilt,zEHU,ipower
c
c     Note: Used integration by parts in order to get a positive
c     integrand.
c      
      k = dexp(kl)
      x = scale*k
      DD0=1.d0
      DD_cbnu=1.d0
      DD_cb=1.d0
      if (ipower.eq.0) then
         t1=TF_master(k,dTFdk,1)
         dsigmatop = -k**(3.d0+tilt)*(DD0*DD_cbnu)**2 *t1*(2.d0*dTFdk*k
     $        +t1*(3.d0+tilt))*(3.d0*(x*dcos(x) - dsin(x))/x**3)**2
     $        /scale
      else
         t1=TF_master(k,dTFdk,1)
         dsigmatop = -k**(3.d0+tilt)*(DD0*DD_cb)**2 *t1*(2.d0*dTFdk*k
     $        +t1*(3.d0+tilt))*(3.d0*(x*dcos(x) - dsin(x))/x**3)**2
     $        /scale
      endif
      return
      end

      subroutine TFset_parameters
      implicit double precision (a-h,k,o-z)
      
      double precision N_nu
      common /GLOBALVARIABLES/ omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu,
     $     beta_c,sound_horizon,theta_cmb,omega,omegal,z_equality

c     Auxiliary variable

      obhh = omhh*f_baryon
c     Main variables

      z_equality = 2.50d4*omhh*theta_cmb**(-4.d0) - 1.D0
      k_equality = 0.0746d0*omhh*theta_cmb**(-2.d0)

      z_drag = 0.313d0*omhh**(-0.419d0)*(1.d0+0.607d0*omhh**(0.674d0))
      z_drag = 1.d0 + z_drag*obhh**(0.238d0*omhh**(0.223d0))
      z_drag = 1291.d0 * omhh**(0.251d0)/
     &     (1d0 + 0.659d0*omhh**(0.828d0)) * z_drag

      y_d = (1.d0+z_equality)/(1.d0+z_drag)

      R_drag = 31.5d0*obhh*theta_cmb**(-4.d0)*1000.d0/(1.d0 + z_drag)
      R_equality = 31.5d0*obhh*theta_cmb**(-4.d0)
     &     *1000.d0/(1.d0 + z_equality)

      sound_horizon = 2.d0/3.d0/k_equality*dsqrt(6.d0/R_equality)*
     &     dlog(( dsqrt(1.d0+R_drag)+dsqrt(R_drag+R_equality) )
     &     /(1.d0+dsqrt(R_equality)))

      p_c  = -(5.d0-dsqrt(1.d0+24.d0*(1.d0-f_nu-f_baryon)))/4.d0
      p_cb = -(5.d0-dsqrt(1.d0+24.d0*(1.d0-f_nu)))/4.d0
      f_c  = 1.d0-f_nu-f_baryon
      f_cb = 1.d0-f_nu
      f_nub= f_nu+f_baryon


      alpha_nu= (f_c/f_cb)* (2.d0*(p_c+p_cb)+5.d0)/(4.d0*p_cb+5.d0)
      alpha_nu= alpha_nu*(1.d0-0.553d0*f_nub+0.126d0*f_nub**3)
      alpha_nu= alpha_nu/(1.d0-0.193d0*dsqrt(f_nu)+0.169d0*f_nu)
      alpha_nu= alpha_nu*(1.d0+y_d)**(p_c-p_cb)
      alpha_nu= alpha_nu*(1.d0+ (p_cb-p_c)/2.d0 *
     *     (1.d0+1.d0/(4.d0*p_c+3.d0)/(4.d0*p_cb+7.d0))/(1.d0+y_d))
      beta_c=1.d0/(1.d0-0.949d0*f_nub)

      return

      end 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function rombint(f,a,b,tol)
c     rombint returns the integral from a to b of using Romberg integration.
c     The method converges provided that f(x) is continuous in (a,b).
c     f must be double precision and must be declared external in the calling
c     routine.  tol indicates the desired relative accuracy in the integral.
c     
c     Taken from CMBFAST; Seljak & Zaldarriaga (1996)
c
      parameter (MAXITER=30,MAXJ=5)
      implicit double precision (a-h,o-z)
      dimension g(MAXJ+1)
      double precision f
      external f
c     
      h=0.5d0*(b-a)
      gmax=h*(f(a)+f(b))
      g(1)=gmax
      nint=1
      error=1.0d20
      i=0
 10   i=i+1
      if (i.gt.MAXITER.or.(i.gt.5.and.abs(error).lt.tol))
     2     go to 40
c     Calculate next trapezoidal rule approximation to integral.
      g0=0.0d0
      do 20 k=1,nint
         g0=g0+f(a+(k+k-1)*h)
 20   continue

      g0=0.5d0*g(1)+h*g0
      h=0.5d0*h
      nint=nint+nint
      jmax=min(i,MAXJ)
      fourj=1.0d0
      do 30 j=1,jmax
c     Use Richardson extrapolation.
         fourj=4.0d0*fourj
         g1=g0+(g0-g(j))/(fourj-1.0d0)
         g(j)=g0
         g0=g1
 30   continue
      if (abs(g0).gt.tol) then
         error=1.0d0-gmax/g0
      else
         error=gmax
      end if
      gmax=g0
      g(jmax+1)=g0
      go to 10
 40   rombint=g0
      if (i.gt.MAXITER.and.abs(error).gt.tol)
     2     write(*,*) 'rombint failed to converge; integral, error=',
     3     rombint,error
      return
      end