ccc pk_halomodel.f   Time-stamp: <2018-06-12 09:49:05 hamana>
c 
c   WAHT IS THIS ?
c   This program computes the halo model prediction for the dark matter
c   power-spectrum.
c   See Hamana, Yoshida & Suto, MNRAS, 568, 462 (2002) and references therein
c   for the description of the halo model.
c
c   Compile: gfortran -fdefault-real-8 -O2 pk_halomodel.f -o pk_halomodel.exe
c
c   Usage: ./pk_halomodel.exe -om {Omega_m value} -h {h value} -ob {Omega_b value} -s8 {sigma_8 value} -z {redshift value} -f {output file name}
c   
c   HISTORY:
c   2018/06/12 --- Ver 2.0 by Takashi Hamana
c   2005/02/07 --- Ver 1.0 by Takashi Hamana
c
c   Bug report to Takashi Hamana
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev
      common /pkshape/ gamfin,atps
      common /defs/ pi,tps,ch0,rn,onethr
      common /st/ p_st,a_st,st_norm
      common /fby/ deltac,deltac2,deltavir3,ga2,zs,rast,ak,as
      common /conc/ b0,c0

      external fnunuby
      external fnunur3y2

      parameter (ikmin=-30)
      parameter (ikmax=40)
      parameter (delk=0.1)

      character*88 filename
      character*256 opt,arg

ccc definitions
      pi=4.0*atan(1.0)
      tps=2.0*pi**2
      ch0=3000.0
      rn=8.0/ch0
      onethr=1.0/3.0

ccc ST mass function
      p_st=0.3   ! 0.0 for PS, 0.3 for ST
      a_st=0.707 ! 1.0 for PS, 0.707 for ST
      
ccc M-c relation by Bullock et al
      c0=8.0
      b0=-0.13

ccc read options
      omem=0.28
      omeb=0.04
      hubp=0.7
      sig8=0.8
      zs=0.0
      filename='pk_halomodel.dat'

! read-in
      narg=iargc()
      do i=1,narg
         call getarg(i,opt)
         call getarg(i+1,arg)
         select case (opt)
         case ('-om')
            read(arg,*)omem
         case ('-h')
            read(arg,*)hubp
         case ('-ob')
            read(arg,*)omeb
         case ('-s8')
            read(arg,*)sig8
         case ('-z')
            read(arg,*)zs
         case ('-f')
            filename=arg
         end select
      enddo

      omev=1.0-omem !!! Flat model only
      ok=0.0
      omebh2=omeb*hubp**2
      
ccc flag
      igrflag=-1

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 construct ST mass function
      call int_st(st_norm)

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

ccc at z=zs, redshift -> scale factor
      as=1.0/(1.0+zs)

ccc growth rate
      g=gr(as,igrflag)
      ga2=(g*as)**2

ccc deltas
      call deltas(as,deltac,deltavir)
      deltac2=deltac**2
      deltavir3=deltavir**onethr

ccc r_ast
      call search_rast(rast,ga2,deltac2)

ccc pk
      open(10,file=filename,STATUS= 'UNKNOWN')
      anumin=-80.0
      anumax=30.0
      do 10 i=ikmin,ikmax
         ak=10.0**(i*delk)
         w=ak*ch0
         pk_lin=delta_lin0(w)*ga2

         call qromb(fnunuby,anumin,anumax,ss1)
         pk_hh=pk_lin*ss1**2
         
         call qromb(fnunur3y2,anumin,anumax,ss2)
         pk_p=ss2*ak**3/(1.5*pi)

         pk_hh_p=pk_hh+pk_p

         write(10,100)ak,pk_lin,pk_hh,pk_p,pk_hh_p
 10   continue
      write(10,99)'# col-1: k [h/Mpc]             '
      write(10,99)'# col-2: pk_lin*k^3/2Pi^2      '
      write(10,99)'# col-3: pk_hh*k^3/2Pi^2       '
      write(10,99)'# col-4: pk_p*k^3/2Pi^2        '
      write(10,99)'# col-5: (pk_hh+pk_p)*k^3/2Pi^2'
 99   format(a31)

      write(10,98)'#  Omega_m = ',omem
      write(10,98)'#  Omega_L = ',omev
      write(10,98)'#  Omega_b = ',omeb
      write(10,98)'#        h = ',hubp
      write(10,98)'#  sigma_8 = ',sig8
      write(10,98)'# redshift = ',zs
 98   format(a13,f8.4)

      close(10)
 100  format(5(1pe13.5))
      end

ccc
      function fnunuby(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /fby/ deltac,deltac2,deltavir3,ga2,zs,rast,ak,as
      common /st/ p,a,st_norm
      anu=exp(x)
      call nu2r_rs(anu,deltac2,deltavir3,
     $     ga2,as,rast,sig,rvir,c_nfw,rs_nfw,rsigma)
      fnunuby=st_norm*fnu_nu(anu)*stbias(anu,deltac)*
     $     yfunc(ak,c_nfw,rs_nfw)
c      if (fnunuby.lt.1.0e-33) fnunuby=0.0
      return
      end

ccc
      function fnunur3y2(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /fby/ deltac,deltac2,deltavir3,ga2,zs,rast,ak,as
      common /st/ p,a,st_norm
      anu=exp(x)
      call nu2r_rs(anu,deltac2,deltavir3,
     $     ga2,as,rast,sig,rvir,c_nfw,rs_nfw,rsigma)
      fnunur3y2=st_norm*fnu_nu(anu)*rsigma**3*
     $     yfunc(ak,c_nfw,rs_nfw)**2
      if (fnunur3y2.lt.1.0e-33) fnunur3y2=0.0
      return
      end

ccc  
      function yfunc(ak,c,rs)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      z=ak*rs
      if (z.gt.1.0e-7) then
         call cisi(z,ciz,siz)
         c1=1.0+c
         zc1=z*c1
         call cisi(zc1,cizc1,sizc1)
         yfunc=(cos(z)*(cizc1-ciz)+sin(z)*(sizc1-siz)
     $        -sin(c*z)/zc1)/(log(c1)-c/c1)
      else 
         yfunc=1.0
      endif
      return
      end

ccc Sheth & Tormen bias
      function stbias(x,deltac)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /st/ p,a,st_norm
      y=a*x
      stbias=1.0+(y-1.0+2.0*p/(1.0+y**p))/deltac
      return
      end

ccc Sheth & Tormen mass function ! NOT narmalized
      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

cccccccccc
      subroutine nu2r_rs(anu,deltac2,deltavir3,
     $     ga2,as,rast,sig,rvir,c_nfw,rs_nfw,rsigma)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /defs/ pi,tps,ch0,rn,onethr
      sig=deltac2/anu
      rsigma=sigma2r(sig,ga2)
      rvir=rsigma/deltavir3
      c_nfw=amss2c(as,rsigma)
      rs_nfw=rvir/c_nfw
      return
      end
      
      function amss2c(as,r)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /conc/ b0,c0
      common /defs/ pi,tps,ch0,rn,onethr
      common /omegas/ omem,omev

ccc Bullock et al model
!!! 2.775e11/1e14=2.775e-3
      amass=2.775e-3*4.0*pi*omem*r**3/3.0
      amss2c=c0*as*amass**b0

      return
      end

ccc delta_c and delta_vir 
ccc VALID only FLAT model
      subroutine deltas(as,deltac,deltavir)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev
      common /defs/ pi,tps,ch0,rn,onethr
      x=as*(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=-80)
      parameter (irmax=40)
      parameter (delr=0.25)
      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=-80)
      parameter (irmax=40)
      parameter (delr=0.25)
      parameter (irall=irmax-irmin+1)

      dimension sigmar0(irall)
      dimension r0(irall)

      do 10 i=1,irall
         ii=irmax-i+1
         r0(i)=delr*ii
         x=exp(r0(i))
         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=-80)
      parameter (irmax=40)
      parameter (delr=0.25)
      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 /pkshape/ 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=-20.0
      wmax=20.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 /pkshape/ 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)
      
      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,igrflag)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /omegas/ omem,omev
      save grnorm
      if (igrflag.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
         igrflag=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
      e=a*omem+a**4*omev
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc integration pakege
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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=20, JMAXP=JMAX+1, K=5, 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
      write(*,*)'too many steps in qromb'
      write(*,*)'STOP'
      stop
      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=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.)then
             write(*,*)'failure in polint'
             write(*,*)'STOP'
             stop
          endif
          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.

ccc 
      SUBROUTINE cisi(x,ci,si)
      INTEGER MAXIT
      REAL ci,si,x,EPS,EULER,PIBY2,FPMIN,TMIN
      PARAMETER (EPS=6.e-8,EULER=.57721566,MAXIT=100,PIBY2=1.5707963,
     *FPMIN=1.e-30,TMIN=2.)
      INTEGER i,k
      REAL a,err,fact,sign,sum,sumc,sums,t,term,absc
      COMPLEX h,b,c,d,del
      LOGICAL odd
      absc(h)=abs(real(h))+abs(aimag(h))
      t=abs(x)
      if(t.eq.0.)then
        si=0.
        ci=-1./FPMIN
        return
      endif
      if(t.gt.TMIN)then
        b=cmplx(1.,t)
        c=1./FPMIN
        d=1./b
        h=d
        do 11 i=2,MAXIT
          a=-(i-1)**2
          b=b+2.
          d=1./(a*d+b)
          c=b+a/c
          del=c*d
          h=h*del
          if(absc(del-1.).lt.EPS)goto 1
11      continue
        write(*,*)'cf failed in cisi'
        write(*,*)'STOP'
        stop
1       continue
        h=cmplx(cos(t),-sin(t))*h
        ci=-real(h)
        si=PIBY2+aimag(h)
      else
        if(t.lt.sqrt(FPMIN))then
          sumc=0.
          sums=t
        else
          sum=0.
          sums=0.
          sumc=0.
          sign=1.
          fact=1.
          odd=.true.
          do 12 k=1,MAXIT
            fact=fact*t/k
            term=fact/k
            sum=sum+sign*term
            err=term/abs(sum)
            if(odd)then
              sign=-sign
              sums=sum
              sum=sumc
            else
              sumc=sum
              sum=sums
            endif
            if(err.lt.EPS)goto 2
            odd=.not.odd
12        continue
          write(*,*)'maxits exceeded in cisi'
          write(*,*)'STOP'
          stop
        endif
2       si=sums
        ci=sumc+log(t)+EULER
      endif
      if(x.lt.0.)si=-si
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.