c------------------------------------------------------------------------------ c Return the sigma_8 appropriate for COBE normalization given c Gamma, Omega_0 and Lambda_0 c c c Uses New fitting formulae from Bunn & White 96 in preparation c 1-sigma error bar now quoted at 7.5% c subroutine cobe_sigma8(Gamma,Omega_0,Lambda_0,nspec,itrans,igwave,sigma_8) **************************************variables******************************** implicit none integer N,i,igwave,itrans real rf,Gamma,Omega_0,Lambda_0,delta_H,EPS,COBE_4YR, x lnkmin,lnkmax,dlnk,sum,sigma_8,pw2k3,lnk,nspec parameter (EPS=1.0e-04,N=100000,COBE_4YR=3.1768) ******************************************************************************* rf=8.0 ! Top-Hat Filter Radius in Mpc/h c c Liddle et al 1996 fitting formulae to the Bunn & White 1997 COBE c 4yr data analysis c The first formula is from Liddle etal 1996 MNRAS 282,281 c equations 10 but the leading constant has a different value c due to the choice of units for wavenumber k and the Fourier c transform convention. if (abs(Lambda_0+Omega_0-1.0).le.EPS) then !flat models delta_H = COBE_4YR * Omega_0**(-0.785-0.05*log(Omega_0)) if (abs(nspec-1.0).ge.EPS) then c The following formulae are Liddle etal 1996 MNRAS 282,281 c equations 11 and 12 if (igwave.eq.0) then delta_H=delta_H*exp(-0.95*(nspec-1.0)-0.169*(nspec-1)**2) else if (igwave.eq.1) then delta_H=delta_H*exp( 1.00*(nspec-1.0)+1.97*(nspec-1)**2) d write(0,*)'assuming gravitational waves for power-law inflation' if (nspec.gt.1.0) stop 'which is not applicable for n>1' else d write(0,*) 'igwave must be 1 or 0' end if end if else if (abs(Lambda_0).le.EPS) then !open models if (abs(nspec-1.0).ge.EPS) then d write(0,*) 'Assuming nspec=1.0 see Liddle,Lyth Roberts and Viana.' nspec=1.0 end if delta_H = COBE_4YR * Omega_0**(-0.35-0.19*log(Omega_0)) else STOP 'Sorry cannot cope with this cosmology' end if c c Integrate k^2 P(k) W(k) and scale by delta_H. lnkmax=5.0-log(rf) lnkmin=-9.0-log(rf) dlnk=(lnkmax-lnkmin)/real(N-1) sum=0.0 do i=1,N lnk=lnkmin+dlnk*real(i-1) sum=sum+pw2k3(exp(lnk),rf,Gamma,nspec,itrans) end do lnkmin=lnkmin-dlnk lnkmax=lnkmax+dlnk sigma_8=(sum+0.5*pw2k3(exp(lnkmin),rf,Gamma,nspec,itrans)+ x 0.5*pw2k3(exp(lnkmax),rf,Gamma,nspec,itrans) )*dlnk sigma_8=delta_H*sqrt(sigma_8) return end c------------------------------------------------------------------------------ c CDM-like power spectrum multiplied by k^3 W^2(k) c c W(k) is the transform of the top-hat window function c c q is the scaled wavenumber k/Gamma and the shape parameter c Gamma = Omega_0 h exp( - Omega_b - Omega_b/Omega_0 ) c real function pw2k3(k,rf,Gamma,nspec,itrans) *************************************variables********************************* implicit none integer itrans real u,q,k,rf,Gamma,trans,win,nspec,KHORIZON parameter(KHORIZON=1.0/3000.0) !defined as c.H_0 in h Mpc^-1 ******************************************************************************* u=k*rf q=k/Gamma if (itrans.eq.0) then c BBKS transfer function trans = (log(1.0+2.34*q)/(2.34*q**2)) x / ( (1.0/q)**4+3.89/q**3 + (16.1/q)**2 + 5.46**3/q + 6.71**4 )**0.25 else c Bond & Efstathiou transfer function trans=(1.0+(((6.4*q)+((3.0*q)**1.5)+((1.7*q)**2.0))**1.13))**(-1.0/1.13) endif win=3.0*(sin(u)/u-cos(u))/u**2 ! top hat window function pw2k3= k**3 * trans**2 * win**2 * (k/KHORIZON)**nspec return end