!----------------------------------------------------------------------------- ! ! Read in, tabulate, and interpolate the linear matter power spectrum. ! The input file is in the following format (same as CAMB format): ! ! [1st column] k/h (i.e., wavenumber in units of h Mpc^-1) ! [2nd column] h^3 P(k) (i.e., P(k) in units of h^-3 Mpc^3) ! ! << sample code >> ! ! USE linearpk ! double precision :: linear_pk, k_ov_h ! character(len=128) :: filename ! integer :: n ! external linear_pk ! filename='wmap5baosn_max_likelihood_matterpower.dat' ! n=896 ! # of lines in the file ! CALL open_linearpk(filename,n) ! k_ov_h=1d0 ! print*,'P(k) at k=',k_ov_h,' h Mpc^-1 is',linear_pk(k_ov_h),' h^3 Mpc^-3' ! end ! ! August 23, 2008 ! E. Komatsu ! !----------------------------------------------------------------------------- MODULE LinearPk IMPLICIT none INTEGER :: jlo,ndata DOUBLE PRECISION, dimension(:), allocatable :: xa,ya,y2a contains SUBROUTINE Open_LinearPk(filename,n) IMPLICIT none DOUBLE PRECISION, dimension(:), allocatable :: aklin,pk_lin CHARACTER(len=128) :: filename INTEGER, intent(IN) :: n INTEGER :: i ndata=n ALLOCATE(xa(ndata),ya(ndata),y2a(ndata)) ALLOCATE(aklin(ndata),pk_lin(ndata)) open(1,file=filename,status='old',form='formatted') print*,'read in '//trim(filename) do i=1,ndata read(1,*)aklin(i),pk_lin(i) ! k/h & h^3 P(k) in CAMB format! enddo close(1) xa=dlog(aklin) ya=dlog(pk_lin) DEALLOCATE(aklin,pk_lin) CALL spline(xa,ya,ndata,1.d30,1.d30,y2a) return END SUBROUTINE Open_LinearPk SUBROUTINE Close_LinearPk DEALLOCATE(xa,ya,y2a) return END SUBROUTINE Close_LinearPk END MODULE LinearPk DOUBLE PRECISION FUNCTION Linear_Pk(ak) Use LinearPk IMPLICIT none DOUBLE PRECISION :: a,b,h,x,y,ak x = dlog(ak) CALL hunt(xa,ndata,x,jlo) h=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/h b=(x-xa(jlo))/h y=a*ya(jlo)+b*ya(jlo+1)+((a**3-a)*y2a(jlo)+(b**3-b)*y2a(jlo+1))*(h**2)/6. Linear_Pk = dexp(y) return END FUNCTION Linear_Pk DOUBLE PRECISION FUNCTION dlnPkdlnk(ak) Use LinearPk IMPLICIT none DOUBLE PRECISION :: a,b,h,x,y,ak x = dlog(ak) CALL hunt(xa,ndata,x,jlo) h=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/h b=(x-xa(jlo))/h y=(ya(jlo+1)-ya(jlo))/h+(-(3.*a**2-1.)*y2a(jlo)+(3.*b**2-1.)*y2a(jlo+1))*h/6. dlnPkdlnk = y return END FUNCTION DlnPkdlnk