PROGRAM CosmicVariance USE cosmo USE growth USE linearpk USE angular_distance ! A code for computing the cosmic variance (sample variance) of ! the mean density for a given survey geometry. More specifically, ! it gives "sigma," which is a fractional error on the measured mean ! density of galaxies within a survey volume. It is calculated as: ! ! sigma^2 = int dkx dky dkz P(sqrt(kx^2+ky^2+kz^2)) ! *[sinc(kxLx)sinc(kuLy)sinc(kzLz)]^2 ! ! where P(k) includes the effect of the linear redshift space ! distortion. If you do not want this, set beta=0d0 in the code by hand. ! December 5, 2010: E.Komatsu IMPLICIT none double precision :: k_ov_h,k0_ov_h,kmax_ov_h ! h Mpc^-1 double precision :: linear_pk,sig2,mu2,z1,z2 double precision :: dk,factor,dummy,thetax,thetay,one_over_h double precision :: pk,pkd,z,zin,beta,bias,g,dgdlna,f,da double precision :: kx,ky,kz,Lx,Ly,Lz ! h^-1 Mpc integer :: n character(len=128) :: filename external linear_pk,dgdlna,g,da,one_over_h ! maximum wavenumber to which P(k) is integrated. kmax_ov_h=0.5d0 ! this works as long as L>10 h^-1 Mpc. ! Specify three cosmological parameters for computing the growth factor ! and the angular diameter distance. ! The data type has been defined in MODULE cosmo. ode0=0.723d0 om0=0.277d0 w=-1d0 CALL setup_growth ! tabulate the growth factor CALL setup_da ! tabulate the proper angular diameter distance ! read in linear P(k) filename='wmap5baosn_max_likelihood_matterpower_at_z=30.dat' n=896 ! # of lines in the file zin=30d0 ! redshift of the input power spectrum CALL open_linearpk(filename,n) ! read in redshift print*,'minimum redshift and maximum redshift? (eg: 1.9 3.5)' read*,z1,z2 z=0.5d0*(z1+z2) print*,'average linear bias? (eg: 2)' read*,bias f=(1d0+dgdlna(z)/g(z)) beta=f/bias ! survey size print*,'angular extensions in x and y directions (in degrees)? (eg: 10 10)' read*,thetax,thetay thetax=thetax*(3.14159d0/180d0) ! radian thetay=thetay*(3.14159d0/180d0) ! radian Lx=thetax*da(z)*(1d0+z) ! h^-1 Mpc, comoving Ly=thetay*da(z)*(1d0+z) ! h^-1 Mpc, comoving Lz=(z2-z1)*(2998d0*one_over_h(z)) print*,'-- survey geometry (comoving) --' print*,'Lx=',Lx,' h^-1 Mpc' print*,'Ly=',Ly,' h^-1 Mpc' print*,'Lz=',Lz,' h^-1 Mpc' print*,'--------------------------------' open(11,file=filename,status='old') read(11,*)k0_ov_h,dummy close(11) sig2=0d0 dk=5d-3 kz=k0_ov_h do while (kz<=kmax_ov_h) ky=k0_ov_h do while (ky<=kmax_ov_h) kx=k0_ov_h do while (kx<=kmax_ov_h) k_ov_h=dsqrt(kx**2d0+ky**2d0+kz**2d0) pkd=linear_pk(k_ov_h)*(g(z)/g(zin)*(1d0+zin)/(1d0+z))**2d0 factor=8d0*dk**3d0/(2d0*3.1415926535d0)**3d0 ! h^3 Mpc^-3 mu2=kz**2d0/k_ov_h**2d0 sig2=sig2+bias**2d0*pkd*(1d0+beta*mu2)**2d0*factor & *(dsin(kx*Lx/2d0)/(kx*Lx/2d0))**2d0 & *(dsin(ky*Ly/2d0)/(ky*Ly/2d0))**2d0 & *(dsin(kz*Lz/2d0)/(kz*Lz/2d0))**2d0 kx=kx+dk enddo ky=ky+dk enddo kz=kz+dk enddo CALL close_linearpk print*,'Cosmic Variance for (delta N)/N =',dsqrt(sig2) print*,'- This is the fractional accuracy with which you can determine' print*,'the mean number density of galaxies in the survey, and should be' print*,'compared to the shot noise given by 1/SQRT(N).' print*,'(Here, N is the number of objects in the survey.)' END PROGRAM CosmicVariance