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.