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.