ccc pskap_pd.f Time-stamp: <2006-08-14 09:02:07 hamanatk> implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv common /dpower/rn,gamfin,sig8,anrm external dw parameter (ismin=10) parameter (ismax=60) parameter (dels=0.1) dimension s(ismin:ismax),dps(ismin:ismax),pskap(ismin:ismax) dimension dpsnl(ismin:ismax),pskapnl(ismin:ismax) parameter (ia=100) dimension asf(0:ia),cd(0:ia),d(0:ia),dod(0:ia) ccc parameters for PD's nonlinear Pk fitting function parameter (klmin=-80,klmax=180) parameter (kl=klmax-klmin+1) parameter (delkl=0.1) parameter (delkpm=0.01) dimension wl(klmin:klmax) dimension wnl(klmin:klmax) dimension dell(klmin:klmax) dimension pnl(klmin:klmax) dimension dpnl(klmin:klmax) dimension a(klmin:klmax) dimension b(klmin:klmax) dimension alp(klmin:klmax) dimension bet(klmin:klmax) dimension v(klmin:klmax) pi=4.0*atan(1.0) twopi=2.0*pi tps=2.0*pi**2 onethr=1.0/3.0 cv=3000.0 igrfrag=-1 write(*,*)'### pskap_pd.f' write(*,*)'### version-1 [14/8/2006]' write(*,*)'###' write(*,*)'### Use this at your own risk' write(*,*)'### Bug report to hamanatk [at-mark] cc.nao.ac.jp' write(*,*)'###' write(*,*)'### This uses BBKS CDM P(k) with ' write(*,*)'### Sugiyama (1995) Gamma factor and' write(*,*)'### Peacock & Dodds (1996) nonlinear fitting function' write(*,*)'###' write(*,*)'### This will output a data file pskap_pd.dat' write(*,*)'###' write(*,*)'### Hubble parameter is defined by' write(*,*)'### H_0=100h[km/sec/Mpc]' write(*,*)'#########################################' 11 write(*,*)'Cosmology 1, 2 or 3 ?' write(*,*)'1: Lambda-flat (Omega_m=0.3, Lambda=0.7)' write(*,*)'2: Eds (Omega_m=1.0, Lambda=0.0)' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then om=0.3 ov=0.7 elseif (icosmo.eq.2) then om=1.0 ov=0.0 elseif (icosmo.eq.3) then write(*,*)'Omega_m ?' read(*,*)om write(*,*)'Lambda ?' read(*,*)ov else goto 11 endif crv=om+ov-1.0 scrv=sqrt(abs(crv)) kcv=0 if (crv.gt.0.001) kcv=1 if (crv.lt.-0.001) kcv=-1 12 write(*,*)'Hubble parameter 1, 2 or 3 ?' write(*,*)'1: h=0.7' write(*,*)'2: h=0.5' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then hubp=0.7 elseif (icosmo.eq.2) then hubp=0.5 elseif (icosmo.eq.3) then write(*,*)'h ?' read(*,*)hubp else goto 12 endif 13 write(*,*)'Baryon density 1, 2 or 3 ?' write(*,*)'1: Omega_b=0.04' write(*,*)'2: Omega_b=0.0' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then ob=0.04 elseif (icosmo.eq.2) then ob=0.0 elseif (icosmo.eq.3) then write(*,*)'Omega_b ?' read(*,*)ob else goto 13 endif ccc gamma factor by Sugiyama (1995) gamfin=exp(ob+sqrt(2.0*hubp)*ob/om)/(om*hubp*cv) 14 write(*,*)'sigma_8 1, 2 or 3 ?' write(*,*)'1: sigma_8=0.9' write(*,*)'2: sigma_8=0.5' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then sig8=0.9 elseif (icosmo.eq.2) then sig8=0.5 elseif (icosmo.eq.3) then write(*,*)'sigma_8 ?' read(*,*)sig8 else goto 14 endif ccccc normalization of power spectrum call normPk coef=9.0/4.0*om*om 15 write(*,*)'source redshift 1, 2 or 3 ?' write(*,*)'1: zs=1.0' write(*,*)'2: zs=2.0' write(*,*)'3: something else' read(*,*)icosmo if (icosmo.eq.1) then zmax=1.0 elseif (icosmo.eq.2) then zmax=2.0 elseif (icosmo.eq.3) then write(*,*)'zs ?' read(*,*)zmax else goto 15 endif ccc source redshift amin=1.0/(1.0+zmax) dela=(1.0-amin)/ia dela2=0.5*dela ccc complute a-w a-r relation asf(0)=1.0 cd(0)=0.0 d(0)=0.0 do j=1,ia asf(j)=1.0-j*dela call qromb(dw,asf(j),1.0,ss) cd(j)=ss d(j)=fk(ss,kcv) enddo wi=cd(ia) fks=1.0/d(ia) do j=0,ia wij=wi-cd(j) dod1=fk(wij,kcv)*fks dod(j)=dod1*dod1 enddo ccc compute parameters of PD's fitting function dwp=0.5*exp(delkpm) dwm=0.5*exp(-delkpm) do ik=klmin,klmax wn=exp(ik*delkl) wl(ik)=wn q=wn*gamfin pslin=anrm*trnf2(q)*wn dell(ik)=pslin*wn**3/tps w2=wn*0.5 w1=wn*dwm w3=wn*dwp delk31=w3-w1 q1=w1*gamfin q2=w2*gamfin q3=w3*gamfin ps1=trnf2(q1)*w1 ps2=trnf2(q2)*w2 ps3=trnf2(q3)*w3 dlplk=w2*(ps3-ps1)/(delk31*ps2) en=1.0/(1.0+dlplk/3.0) a(ik)=0.482*en**0.947 b(ik)=0.226*en**1.778 alp(ik)=3.310*en**0.244 bet(ik)=0.862*en**0.287 v(ik)=11.55*en**0.423 enddo ccc set initial do i=ismin,ismax s(i)=10.0**(i*dels) pskap(i)=0.0 dps(i)=0.0 pskapnl(i)=0.0 dpsnl(i)=0.0 enddo ccc integration over z do i=1,ia z=1.0/asf(i)-1.0 gr2=gr(asf(i),igrfrag)**2 ga2=gr2*asf(i)**2 g3=grpd(asf(i))**3 weight=dod(i)/sqrt(e(asf(i))) weightnl=weight/asf(i)**2 ccc PD's fitting function of the nonlinear Pk do j=klmin,klmax delnl2=fnl(a(j),b(j),alp(j),bet(j),v(j),dell(j)*ga2,g3) wknl=wl(j)*(1.0+delnl2)**onethr pnl(j)=delnl2*tps/wknl**3 wnl(j)=wknl enddo dpdkmin=(pnl(klmin+1)-pnl(klmin))/ $ (wnl(klmin+1)-wnl(klmin)) dpdkmax=(pnl(klmax)-pnl(klmax-1))/ $ (wnl(klmax)-wnl(klmax-1)) call spline(wnl,pnl,kl,dpdkmin,dpdkmax,dpnl) do ks=ismin,ismax ak=s(ks)/d(i) !!! linear dps2=weight*gr2*plin(ak) pskap(ks)=pskap(ks)+ $ (dps(ks)+dps2)*dela2 dps(ks)=dps2 !!! non-linear if (ak.le.wnl(klmin)) then y=pnl(klmin)*(ak/wnl(klmin)) elseif (ak.ge.wnl(klmax)) then y=pnl(klmax)*(ak/wnl(klmax))**(-2.5) else call splint(wnl,pnl,dpnl,kl,ak,y) endif dps2=weightnl*y pskapnl(ks)=pskapnl(ks)+ $ (dpsnl(ks)+dps2)*dela2 dpsnl(ks)=dps2 enddo enddo open(10,file='pskap_pd.dat',status='unknown') write(10,'(a41)')'# s; Plin; s^2*Plin/2pi; Pnl; s^2*Pnl/2pi' do ks=ismin,ismax ps2d=coef*pskap(ks) delps2d=ps2d*s(ks)**2/twopi ps2dnl=coef*pskapnl(ks) delps2dnl=ps2dnl*s(ks)**2/twopi write(10,'(5(1pe13.5))')s(ks),ps2d,delps2d,ps2dnl,delps2dnl enddo close(10) end !!! function used in the PD fitting function function fnl(a,b,alp,bet,v,x,g3) implicit real (a-h,o-z) implicit integer (i-n) fnl=x*((1.0+b*bet*x+(a*x)**(alp*bet))/ $ (1.0+((a*x)**alp*g3/(v*sqrt(x)))**bet))**(1.0/bet) return end c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ccc CDM power spectrum function plin(w) implicit real (a-h,o-z) implicit integer (i-n) common /dpower/rn,gamfin,sig8,anrm q=w*gamfin plin=anrm*trnf2(q)*w return end ccc CDM transfer function by BBKS 1986 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 ccc linear growth factor function gr(a,igrfrag) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv save grnorm if (igrfrag.eq.-1) then et=1.0/e(1.0) ot=om*et vt=ov*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=om*a*et vt=ov*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 linear grwoth rate function grpd(a) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv et=1.0/e(a) ot=om*a*et vt=ov*et*a**4 graw=2.5*ot/(ot**0.571428571-vt+ $ (1.0+0.5*ot)*(1.0+vt*0.014285714)) grpd=graw return end ccc function dw(a) implicit real (a-h,o-z) implicit integer (i-n) dw=1.0/sqrt(e(a)) return end ccc function e(a) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv a2=a*a e=a*om-a2*crv+a2*a2*ov return end ccc function fk(w,kcv) implicit real (a-h,o-z) implicit integer (i-n) common /cosmo/ om,ov,crv,scrv if (kcv.eq.-1) then fk=sinh(scrv*w)/scrv elseif (kcv.eq.0) then fk=w else fk=sin(scrv*w)/scrv endif return end ccccccc normalization of power spectrum subroutine normPk implicit real (a-h,o-z) implicit integer (i-n) common /dpower/rn,gamfin,sig8,anrm external delta pi=4.0*atan(1.0) coe=4.5/pi/pi rn=8.0/3000.0 wmin=-10.0 wmax=10.0 call qromb(delta,wmin,wmax,seki) y=coe*seki anrm=sig8*sig8/y return end ccccccccc k k^2 t^2 W ccccccccccccccccccccccccccccccccccc function delta(x) implicit real (a-h,o-z) implicit integer (i-n) common /dpower/rn,gamfin,sig8,anrm w=exp(x) q=w*gamfin r=w*rn filt=(sin(r)-r*cos(r))/r**3 tranf=trnf2(q) delta=w**4*filt*filt*tranf return end cccccc integration pakege cccccccccccccccccccccccccccccccc 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 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 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. cccccc spline fitting ccccccccccccccccccccccccccccccccc SUBROUTINE splint(xa,ya,y2a,n,x,y) INTEGER n REAL x,y,xa(n),y2a(n),ya(n) INTEGER k,khi,klo REAL a,b,h klo=1 khi=n 1 if (khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x)then khi=k else klo=k endif goto 1 endif h=xa(khi)-xa(klo) if (h.eq.0.) pause 'bad xa input in splint' a=(xa(khi)-x)/h b=(x-xa(klo))/h y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** *2)/6. return END C (C) Copr. 1986-92 Numerical Recipes Software 41m. SUBROUTINE spline(x,y,n,yp1,ypn,y2) INTEGER n,NMAX REAL yp1,ypn,x(n),y(n),y2(n) PARAMETER (NMAX=5000) INTEGER i,k REAL p,qn,sig,un,u(NMAX) if (yp1.gt..99e30) then y2(1)=0. u(1)=0. else y2(1)=-0.5 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do 11 i=2,n-1 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) p=sig*y2(i-1)+2. y2(i)=(sig-1.)/p u(i)=(6.*((y(i+1)-y(i))/(x(i+ *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig* *u(i-1))/p 11 continue if (ypn.gt..99e30) then qn=0. un=0. else qn=0.5 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) do 12 k=n-1,1,-1 y2(k)=y2(k)*y2(k+1)+u(k) 12 continue return END C (C) Copr. 1986-92 Numerical Recipes Software 41m.