ccc xi_halomodel.f Time-stamp: <2005-02-07 11:42:18 hamana> c c WAHT IS THIS ? c This program computes the halo model prediction for the dark matter c two-point correlation functions. c See Hamana, Yoshida & Suto, MNRAS, 568, 462 (2002) and references therein c for the description of the halo model. c c If this does not work, try a complile option that re-defines REAL as c DOUBLE PRECISION (REAL*8) such like -r8 or -dbl or -Ad. c c HISTORY: 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=-1000) parameter (ikmax=1000) parameter (delk=0.0025) parameter (ixmin=-20) parameter (ixmax=20) parameter (delx=0.1) dimension delhh(ikmin:ikmax) dimension delp(ikmin:ikmax) character*88 filename 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 11 write(*,*)'Cosmology 1 or 2 ?' write(*,*)'1: Lambda-flat (Omega_m=0.3, Lambda=0.7)' write(*,*)'2: something else' read(*,*)icosmo if (icosmo.eq.1) then omem=0.3 omev=0.7 elseif (icosmo.eq.2) then write(*,*)'Omega_m ?' read(*,*)omem write(*,*)'Lambda ?' read(*,*)omev else goto 11 endif ok=abs(omem+omev-1.0) if (ok.ge.1.0e-3) then write(*,*)'Flat model only' goto 11 endif 12 write(*,*)'Hubble parameter (H_0=100h) 1 or 2 ?' write(*,*)'1: h=0.7' write(*,*)'2: something else' read(*,*)icosmo if (icosmo.eq.1) then hubp=0.7 elseif (icosmo.eq.2) then write(*,*)'h ?' read(*,*)hubp else goto 12 endif 13 write(*,*)'Baryon density 1, 2 or 3 ?' write(*,*)'1: Omega_b=0.0' write(*,*)'2: Omega_b=0.04' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then omeb=0.0 elseif (icosmo.eq.2) then omeb=0.04 elseif (icosmo.eq.3) then write(*,*)'Omega_b ?' read(*,*)omeb else goto 13 endif omebh2=omeb*hubp**2 14 write(*,*)'sigma_8 1 or 2 ?' write(*,*)'1: sigma_8=0.9' write(*,*)'2: something else' read(*,*)icosmo if (icosmo.eq.1) then sig8=0.9 elseif (icosmo.eq.2) then write(*,*)'sigma_8 ?' read(*,*)sig8 else goto 14 endif write(*,*)'Redshift ?' read(*,*)zs write(*,*)'Output file name' read(*,'(a88)')filename ccc flag igrfrag=-1 ipf=0 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,igrfrag) 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 anumin=-30.0 anumax=10.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) delhh(i)=pk_lin*ss1**2 call qromb(fnunur3y2,anumin,anumax,ss2) delp(i)=(ss2+ss3)*ak**3/(1.5*pi) 10 continue ccc xi open(10,file=filename,STATUS= 'UNKNOWN') ipf=0 do 20 ix=ixmin,ixmax x=10.0**(ix*delx) xihh=0.0 xip=0.0 do 21 ik=ikmin,ikmax-1 akr=x*10.0**(ik*delk) wind=sin(akr)/akr xihh=xihh+delhh(ik)*wind xip=xip+delp(ik)*wind 21 continue xihh=xihh*delk*log(10.0) xip=xip*delk*log(10.0) if (xip.le.0.0) ipf=1 if (ipf.eq.1) xip=0.0 xihalo=xihh+xip write(10,100)x,xihh,xip,xihalo 20 continue write(10,99)'# col-1: r [Mpc/h] ' write(10,99)'# col-2: xi_2halo ' write(10,99)'# col-3: xi_1halo ' write(10,99)'# col-4: xi_halomodel' 99 format(a21) close(10) 100 format(4(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 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 c write(*,*)sig,ga2 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=-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 /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,igrfrag) implicit real (a-h,o-z) implicit integer (i-n) common /omegas/ omem,omev save grnorm 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 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 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=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 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. 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 pause 'cf failed in cisi' 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 pause 'maxits exceeded in cisi' 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.