ccc  massfunc.f    Time-stamp: <2004-06-01 08:21:00 hamana>
ccc  If this does not work well, try the compile option that
ccc  sets default size of REAL to 8 bytes, such like
ccc  "-r8" or  "-real_size 64"
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev,omek
      common /matter_Pk/ gamfin,atps
      common /defs/ pi,tps,ch0,rn,onethr
      common /st/ p_st,a_st,st_norm

      parameter (numin=-25)
      parameter (numax=25)
      parameter (delnu=0.1)

      character*99 fname

ccc useful numbers
      pi=4.0*atan(1.0)
      tps=2.0*pi**2
      ch0=3000.0
      rn=8.0/ch0
      onethr=1.0/3.0
      e2ten=2.30259 !ln(x)=ln(10)*log(x), ln(10)=e2ten
      zero=0.0
      rhocr=2.7755e11 ! M_sun h^2 /Mpc^3
      
ccc flags
      igrfrag=-1

ccc 
      write(*,*)'### massfunc.f'
      write(*,*)'### version-1 [20/06/2003]'
      write(*,*)'###'
      write(*,*)'### Use this at your own risk' 
      write(*,*)'### Bug report to Takashi Hamana'
      write(*,*)'###'
      write(*,*)'### If this does not work well'
      write(*,*)'### try a compile option that'
      write(*,*)'### defines REAL as DOUBLE PRECISION (REAL*8)'
      write(*,*)'### such like "-r8" or "-real_size 64"'
      write(*,*)'###'
      write(*,*)'### This uses BBKS CDM transfer function'
      write(*,*)'### with Sugiyama95 fitting function for Gamma'
      write(*,*)'###'
      write(*,*)'### NOTE **Valid only for Flat cosmology**'
      write(*,*)'###        Omega_m+Lambda should be 1'
      write(*,*)'#########################################'

ccc mass function
 901  write(*,*)'[1] ST99 [2] PS'
      read(*,*)imodel
      if (imodel.eq.2) then
         p_st=0.0  ! PS parameter 
         a_st=1.0  ! PS parameter
      elseif (imodel.eq.1) then
         p_st=0.3    ! ST parameter 
         a_st=0.707  ! ST parameter 
      else
         goto 901
      endif

ccc parameter input
 902  write(*,*)'cosmological parameters'
      write(*,*)'[1] LCDM' 
      write(*,*)'    Omega_m=0.3'
      write(*,*)'    Lambda=0.7'
      write(*,*)'    Omega_b=0.0'
      write(*,*)'    h=0.7'
      write(*,*)'    sigma_8=0.9'
      write(*,*)'[2] SCDM' 
      write(*,*)'    Omega_m=1.0'
      write(*,*)'    Lambda=0.0'
      write(*,*)'    Omega_b=0.0'
      write(*,*)'    h=0.5'
      write(*,*)'    sigma_8=0.5'
      write(*,*)'[3] else (flat model only)'
      read(*,*)icosmo
      if (icosmo.eq.1) then
         hubp=0.7
         omem=0.3
         omev=0.7
         omeb=0.0
         sig8=0.9
      elseif (icosmo.eq.2) then
         hubp=0.5
         omem=1.0
         omev=0.0
         omeb=0.0
         sig8=0.5
      elseif (icosmo.eq.3) then
         write(*,*)'Omega_m ?'
         read(*,*)omem
         write(*,*)'Lambda ?'
         read(*,*)omev
         write(*,*)'Omage_b ?'
         read(*,*)omeb
         write(*,*)'h ?'
         read(*,*)hubp
         write(*,*)'sigma_8'
         read(*,*)sig8
      else
         goto 902
      endif
      omebh2=omeb*hubp**2
      omek=omem+omev-1.0
      rho0=omem*rhocr

ccc realization redshift 
      write(*,*)'redshift ?'
      read(*,*)zs
      as=1.0/(1.0+zs)

ccc Gamma factor
      gamfin=exp(omeb+sqrt(2.0*hubp)*omeb/omem)/
     $     (omem*hubp*ch0)

ccc normalization
      call normPk(anrm,sig8)
      atps=anrm/tps
      
ccc output file name
      write(*,*)'output file name ?'
      read(*,'(a99)')fname

ccc construct mass function
      call int_st(st_norm)

ccc construct sigma-R relation at z=0
      call int_sigmar

ccc growth rate at zs
      g=gr(as,igrfrag)
      ga2=(g*as)**2

ccc Virial overdensity delta_vir at zs
      call deltas(as,deltac,deltavir)
      deltac2=deltac**2
      deltavir3=deltavir**onethr

ccc r_ast at zs
      call search_rast(rast,ga2,deltac2)

      open(10,file=fname,status='unknown')
      write(10,*)'# M; dn(M)/dM; M*dn(M)/dM; dn(M)/dlog10M'
      do 30 inu=numin,numax     ! loop for nu
cccccc for nu=nu(inu)
         anu=10.0**(inu*delnu)  ! nu
         sigma2_nu=deltac2/anu  ! sigma^2
         rnu=sigma2r(sigma2_nu,ga2) ! R

         amss=4.0*pi*rho0*rnu**3/3.0
         coef=st_norm*rho0/amss

         stmf=coef*fnu_nu(anu)/anu

         anu1=10.0**((inu+0.1)*delnu)
         sigma2_nu=deltac2/anu1
         rnu=sigma2r(sigma2_nu,ga2)
         amass1=4.0*pi*rho0*rnu**3/3.0
         
         anu2=10.0**((inu-0.1)*delnu)
         sigma2_nu=deltac2/anu2
         rnu=sigma2r(sigma2_nu,ga2)
         amass2=4.0*pi*rho0*rnu**3/3.0
         
         dnudm=(anu1-anu2)/(amass1-amass2)
         stmf0=stmf*dnudm
         stmfm=stmf*dnudm*amss
         dnudlog10m=(anu1-anu2)/(log10(amass1)-log10(amass2))
         stmflog10=stmf*dnudlog10m
         
         write(10,100)amss,stmf0,stmfm,stmflog10
 30   continue
 100  format(4(1pe13.5))
      close(10)
      end


ccc Sheth & Tormen mass function ! NOT normalized
      function fnu_nu(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /st/ p,a,st_norm
      y=a*x
      fnu_nu=(1.0+1.0/y**p)*sqrt(y)*exp(-0.5*y)
      return
      end

ccc normalization of Sheth & Tormen mass function    
      subroutine int_st(st_norm)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      external fnu_nu_log
      amin=-80.0
      amax=80.0
      call qromb(fnu_nu_log,amin,amax,ss)
      st_norm=1.0/ss
      return
      end
      function fnu_nu_log(z)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /st/ p,a,st_norm
      x=exp(z)
      y=a*x
      fnu_nu_log=(1.0+y**(-p))*sqrt(y)*exp(-0.5*y)
      return
      end

ccc delta_c and delta_vir 
ccc VALID only FLAT model
      subroutine deltas(a,deltac,deltavir)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev,omek
      common /defs/ pi,tps,ch0,rn,onethr
      x=a*(1.0/omem-1.0)**(1.0/3.0)
      deltac=3.0*(12.0*pi)**(2.0/3.0)/20.0*
     $     (1.0-0.00123*log10(1.0+x**3))
      deltavir=18.0*pi**2*(1.0+0.4093*x**2.71572)
      return
      end


      subroutine search_rast(rast,ga2,deltac2)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /sigma_r_relation/ sigmar0,r0,drdsig1,drdsig2
      parameter (irmin=-350)
      parameter (irmax=250)
      parameter (delr=0.02)
      parameter (irall=irmax-irmin+1)
      dimension sigmar0(irall)
      dimension r0(irall)
      sig0=log(deltac2/ga2)
      call hunt(sigmar0,irall,sig0,jlo)
      d=(r0(jlo+1)-r0(jlo))/(sigmar0(jlo+1)-sigmar0(jlo))
      r1=d*(sig0-sigmar0(jlo))+r0(jlo)
      rast=exp(r1)
      return
      end
ccc sigma_R - R relation
      subroutine int_sigmar
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /sigma_r_relation/ sigmar0,r0,drdsig1,drdsig2
      parameter (irmin=-350)
      parameter (irmax=250)
      parameter (delr=0.02)
      parameter (irall=irmax-irmin+1)
      dimension sigmar0(irall)
      dimension r0(irall)
      do 10 i=1,irall
         ii=irmax-i+1
         x=10.0**(delr*ii)
         r0(i)=log(x)
         sigmar0(i)=log(sigma_r(x)) ! Notice sigma_r(x)=sigma_R^2
 10   continue
      drdsig1=(r0(2)-r0(1))/(sigmar0(2)-sigmar0(1))
      drdsig2=(r0(irall)-r0(irall-1))/
     $     (sigmar0(irall)-sigmar0(irall-1))
      return
      end
ccc
      function sigma2r(sig,ga2)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /sigma_r_relation/ sigmar0,r0,drdsig1,drdsig2
      parameter (irmin=-350)
      parameter (irmax=250)
      parameter (delr=0.02)
      parameter (irall=irmax-irmin+1)
      dimension sigmar0(irall)
      dimension r0(irall)
      sig0=log(sig/ga2)
      if (sig0.le.sigmar0(1)) then
         y=r0(1)+drdsig1*(sig0-sigmar0(1))
      elseif (sig0.ge.sigmar0(irall)) then
         y=r0(irall)+drdsig2*(sig0-sigmar0(irall))
      else
         call hunt(sigmar0,irall,sig0,jlo)
         d=(r0(jlo+1)-r0(jlo))/(sigmar0(jlo+1)-sigmar0(jlo))
         y=d*(sig0-sigmar0(jlo))+r0(jlo)
      endif
      sigma2r=exp(y)
      return
      end


ccccccccccccccccccccccccccccccccccccccccccccc
      function sigma_r(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /defs/ pi,tps,ch0,rn,onethr
      common /sigmar/rsmooth
      external dsigma_r
      wmin=-30.0
      wmax=40.0
      rsmooth=x/ch0
      call qromb(dsigma_r,wmin,wmax,ss)
      sigma_r=ss
      return
      end   
ccccccccccccccccccccccccccccccccccccccccccccc
      function dsigma_r(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /sigmar/rsmooth
      w=exp(x)
      rk=rsmooth*w
      dsigma_r=delta_lin0(w)*filtTH(rk)**2
      return
      end 
ccccccccccccccccccccccccccccccccccccccccccccc
      function delta_lin0(w)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /matter_Pk/ gfin,atps
      q=w*gfin
      delta_lin0=atps*trnf2(q)*w**4
      return
      end 

cccccccccc normalization of power spectrum 
      subroutine normPk(a_norm,sig8)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      external delta
      common /defs/ pi,tps,ch0,rn,onethr
      wmin=-30.0
      wmax=40.0
      call qromb(delta,wmin,wmax,seki)
      y=seki/tps
      a_norm=sig8*sig8/y
      return
      end
ccccccccc
      function delta(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /defs/ pi,tps,ch0,rn,onethr
      common /matter_Pk/ gfin,atps
      w=exp(x)
      q=w*gfin
      r=w*rn
      filt=filtTH(r)
      tranf=trnf2(q)
      delta=w**4*filt*filt*tranf
      return
      end
ccccccccc
      function filtTH(r)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      filtTH=3.0*(sin(r)-r*cos(r))/r**3
      return
      end
ccccccccc 
      function trnf2(q)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /matter_Pk/ gamfin,atps
      qq=q*q
      trnf2=0.182628*log(1.0+2.34*q)**2/qq
     $     /sqrt(1.0+3.89*q+259.21*qq
     $     +162.77*qq*q+2027.17*qq*qq)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      function gr(a,igrfrag)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev,omek
      if (igrfrag.eq.-1) then
         pre=1.0
         et=1.0/e(pre)
         ot=omem*et
         vt=omev*et
         grnorm=0.4*(ot**(0.571428571)-vt+
     $        (1.0+0.5*ot)*(1.0+vt*0.014285714))/ot
         igrfrag=1
      endif
      et=1.0/e(a)
      ot=omem*a*et
      vt=omev*et*a**4
      graw=2.5*ot/(ot**(0.571428571)-vt+
     $     (1.0+0.5*ot)*(1.0+vt*0.014285714))
      gr=graw*grnorm
      return
      end
ccc      
      function e(a)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev,omek
      a2=a*a
      e=a*omem-a2*omek+a2*a2*omev
      return
      end

ccc QROMB package
      SUBROUTINE qromb(func,a,b,ss)
      INTEGER JMAX,JMAXP,K,KM
      REAL a,b,func,ss,EPS
      EXTERNAL func
      PARAMETER (EPS=1.e-6, JMAX=40, JMAXP=JMAX+1, K=10, KM=K-1)
CU    USES polint,trapzd
      INTEGER j
      REAL dss,h(JMAXP),s(JMAXP)
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd(func,a,b,s(j),j)
        if (j.ge.K) then
          call polint(h(j-KM),s(j-KM),K,0.,ss,dss)
          if (abs(dss).le.EPS*abs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
11    continue
      pause 'too many steps in qromb'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE polint(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=20)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)pause 'failure in polint'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE trapzd(func,a,b,s,n)
      INTEGER n
      REAL a,b,s,func
      EXTERNAL func
      INTEGER it,j
      REAL del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.

      SUBROUTINE hunt(xx,n,x,jlo)
      INTEGER jlo,n
      REAL x,xx(n)
      INTEGER inc,jhi,jm
      LOGICAL ascnd
      ascnd=xx(n).gt.xx(1)
      if(jlo.le.0.or.jlo.gt.n)then
        jlo=0
        jhi=n+1
        goto 3
      endif
      inc=1
      if(x.ge.xx(jlo).eqv.ascnd)then
1       jhi=jlo+inc
        if(jhi.gt.n)then
          jhi=n+1
        else if(x.ge.xx(jhi).eqv.ascnd)then
          jlo=jhi
          inc=inc+inc
          goto 1
        endif
      else
        jhi=jlo
2       jlo=jhi-inc
        if(jlo.lt.1)then
          jlo=0
        else if(x.lt.xx(jlo).eqv.ascnd)then
          jhi=jlo
          inc=inc+inc
          goto 2
        endif
      endif
3     if(jhi-jlo.eq.1)return
      jm=(jhi+jlo)/2
      if(x.gt.xx(jm).eqv.ascnd)then
        jlo=jm
      else
        jhi=jm
      endif
      goto 3
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.

      FUNCTION rtflsp(func,x1,x2,xacc)
      INTEGER MAXIT
      REAL rtflsp,x1,x2,xacc,func
      EXTERNAL func
      PARAMETER (MAXIT=1000)
      INTEGER j
      REAL del,dx,f,fh,fl,swap,xh,xl
      fl=func(x1)
      fh=func(x2)
      if(fl*fh.gt.0.) pause 'root must be bracketed in rtflsp'
      if(fl.lt.0.)then
        xl=x1
        xh=x2
      else
        xl=x2
        xh=x1
        swap=fl
        fl=fh
        fh=swap
      endif
      dx=xh-xl
      do 11 j=1,MAXIT
        rtflsp=xl+dx*fl/(fl-fh)
        f=func(rtflsp)
        if(f.lt.0.) then
          del=xl-rtflsp
          xl=rtflsp
          fl=f
        else
          del=xh-rtflsp
          xh=rtflsp
          fh=f
        endif
        dx=xh-xl
        if(abs(del).lt.xacc.or.f.eq.0.)return
11    continue
      pause 'rtflsp exceed maximum iterations'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.

CCC QROMB2 
      SUBROUTINE qromb2(func,a,b,ss)
      INTEGER JMAX,JMAXP,K,KM
      REAL a,b,func,ss,EPS
      EXTERNAL func
      PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
CU    USES polint2,trapzd2
      INTEGER j
      REAL dss,h(JMAXP),s(JMAXP)
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd2(func,a,b,s(j),j)
        if (j.ge.K) then
          call polint2(h(j-KM),s(j-KM),K,0.,ss,dss)
          if (abs(dss).le.EPS*abs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
11    continue
      pause 'too many steps in qromb2'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE polint2(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)pause 'failure in polint2'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE trapzd2(func,a,b,s,n)
      INTEGER n
      REAL a,b,s,func
      EXTERNAL func
      INTEGER it,j
      REAL del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.


CCC QROMB3
      SUBROUTINE qromb3(func,z,a,b,ss)
      INTEGER JMAX,JMAXP,K,KM
      REAL a,b,func,ss,EPS,z
      EXTERNAL func
      PARAMETER (EPS=1.e-5, JMAX=100, JMAXP=JMAX+1, K=5, KM=K-1)
CU    USES polint3,trapzd3
      INTEGER j
      REAL dss,h(JMAXP),s(JMAXP)
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd3(func,z,a,b,s(j),j)
        if (j.ge.K) then
          call polint3(h(j-KM),s(j-KM),K,0.,ss,dss)
          if (abs(dss).le.EPS*abs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
11    continue
      pause 'too many steps in qromb3'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE polint3(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)pause 'failure in polint3'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE trapzd3(func,z,a,b,s,n)
      INTEGER n
      REAL a,b,s,func,z
      EXTERNAL func
      INTEGER it,j
      REAL del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(z,a)+func(z,b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(z,x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.