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.