!Define the data types and read/writes them to disk. Also change l_max here.
! March 2007: modified by Julien Lesgourgues in order to get both
! linear and non-linear power spectra. All changes can be identified
! by searching for "LAPTH" comments.
module cmbtypes
use settings
implicit none
!Number of Cls, 1 for just temperature, 3 (4) for polarization (with B)
integer, parameter :: num_cls = 3
!l_max. Tensors are not computed unless compute_tensors = T in input file
!Make these multiples of 50, should be 50 more than you need accurately
integer, parameter :: lmax = 2100, lmax_tensor = 400
!Parameters for calculating/storing the matter power spectrum
!Note that by default everything is linear
! integer, parameter :: num_matter_power = 49 !number of points computed in matter power spectrum
! real, parameter :: matter_power_minkh = 3.16227766E-03 !minimum value of k/h to store
integer, parameter :: num_matter_power = 101 !number of points computed in matter power spectrum
real, parameter :: matter_power_minkh = 0.66e-4 !1e-4 !minimum value of k/h to store
real, parameter :: matter_power_dlnkh = 0.143911568 !log spacing in k/h
real, parameter :: matter_power_maxz = 3. ! 0
integer, parameter :: matter_power_lnzsteps = 10 ! 1
!Only used in params_CMB
real :: pivot_k = 0.05 !Point for defining primordial power spectra
logical :: inflation_consistency = .false. !fix n_T or not
!Note these are the interpolated/extrapolated values. The k at which matter power is computed up to
!by CAMB is set in CMB_Cls_xxx with, e.g. P%Transfer%kmax = 0.6 (which is enough for 2dF)
!Number of scalar-only cls
!if num_cls=4 and CMB_lensing then increased to 4
integer :: num_clsS=min(num_cls,3)
integer, parameter :: norm_As=1, norm_amp_ratio=2
Type CMBParams
real norm(1:num_norm)
!These are fast parameters controling amplitudes, calibrations, etc.
real InitPower(1:num_initpower)
!These are fast paramters for the initial power spectrum
!Now remaining (non-independent) parameters
real omb, omc, omv, omnu, omk, omdm
real ombh2, omch2, omnuh2, omdmh2
real zre, nufrac
real h, H0
real w
real reserved(5)
! added by JL for Viel-SDSS-Lya
real LyaSDSSParams(26)
real A, B, C
end Type CMBParams
Type CosmoTheory
real Age, r10
real SN_loglike, reserved(3)
real cl(lmax,num_cls), cl_tensor(lmax_tensor,num_cls) !TT, TE, EE and BB in that order
real sigma_8
real matter_power(num_matter_power,matter_power_lnzsteps)
!second index is redshifts from 0 to matter_power_maxz, with equal spacing in
!log(1+z) and matter_power_lnzsteps points
! LAPTH NEW:
! linear matter power spectrum today (for lya)
real matter_power_linear(num_matter_power,matter_power_lnzsteps)
! zone for storing array of z, r, and dz/dr
real z_and_r(3,matter_power_lnzsteps)
! END LAPTH NEW
end Type CosmoTheory
logical, parameter :: Old_format = .false.
logical, parameter :: write_all_Cls = .false.
!if false use CAMB's flat interpolation scheme (lossless if models are flat except near lmax when lensed)
contains
subroutine WriteModel(i,CMB, T, like, mult)
integer i
real, intent(in), optional :: mult
Type(CosmoTheory) T
real like, amult
Type(CMBParams) CMB
integer j
if (present(mult)) then
amult = mult
else
amult = 1
end if
if (Old_format) then
stop 'old not supported'
else
j = 0 !format ID
if (write_all_cls) j=1
write(i) j
write(i) amult, num_matter_power, lmax, lmax_tensor, num_cls
write(i) T%SN_loglike, T%reserved
write(i) like
write(i) CMB
write(i) T%Age, T%r10, T%sigma_8, T%matter_power
if (write_all_cls) then
write(i) T%cl(2:lmax,1:num_cls)
write(i) T%cl_tensor(2:lmax_tensor,1:num_cls)
else
!Use interpolation scheme CAMB uses for flat models
!If using significantly non-flat, or increasing interpolation accuracy, save all th cls instead
write(i) T%cl(2:20,1:num_cls)
do j=30,90,10
write(i) T%cl(j,1:num_cls)
end do
do j=110,130, 20
write(i) T%cl(j,1:num_cls)
end do
do j=150,lmax, 50
write(i) T%cl(j,1:num_cls)
end do
if (lmax_tensor /= 0) then
if (lmax_tensor<150) stop 'lmax_tensor too small'
write(i) T%cl_tensor(2:20,1:num_cls)
do j=30,90,10
write(i) T%cl_tensor(j,1:num_cls)
end do
do j=110,130,20
write(i) T%cl_tensor(j,1:num_cls)
end do
do j=150,lmax_tensor, 50
write(i) T%cl_tensor(j,1:num_cls)
end do
end if
end if
end if
if (flush_write) call FlushFile(i)
end subroutine WriteModel
subroutine ReadModel(i,CMB, T, mult, like, error)
integer, intent(in) :: i
integer, intent(out) :: error
real, intent(out) :: mult
Type(CosmoTheory) T
real like
Type(CMBParams) CMB
real icl(lmax,1:num_cls)
integer allcl,j,ind, ix(lmax)
integer almax,almaxtensor, anumpowers, anumcls
error = 0
if (old_format) then
stop 'old not supported'
else
read(i,end=100,err=100) allcl
if (allcl/=0 .and. allcl/=1) stop 'wrong file format'
read(i,end=100,err=100) mult,anumpowers,almax, almaxtensor, anumcls
if (almax > lmax) stop 'reading file with larger lmax'
if (anumcls > num_cls) stop 'reading file with more Cls (polarization)'
read(i) T%SN_loglike, T%reserved
read(i,end = 100, err=100) like
read(i) CMB
read(i) T%Age, T%r10, T%sigma_8, T%matter_power(1:anumpowers,1:matter_power_lnzsteps)
T%cl = 0
T%cl_tensor = 0
if(allcl==1) then
read(i) T%cl(2:almax,1:anumcls)
read(i) T%cl_tensor(2:almaxtensor,1:anumcls)
else
read(i) icl(1:19,1:num_cls)
ind =1
do j =2,20
ix(ind)=j
ind=ind+1
end do
do j=30,90,10
read(i) icl(ind,1:num_cls)
ix(ind) = j
ind = ind + 1
end do
do j=110,130,20
read(i) icl(ind,1:num_cls)
ix(ind) = j
ind = ind + 1
end do
do j=150,almax, 50
read(i) icl(ind,1:num_cls)
ix(ind) = j
ind = ind+1
end do
ind = ind-1
call InterpCls(ix,icl, T%cl, ind, almax, num_Cls)
if (almaxtensor /= 0) then
read(i) icl(1:19,1:num_cls)
ind =1
do j =2,20
ix(ind)=j
ind=ind+1
end do
do j=30,90,10
read(i) icl(ind,1:num_cls)
ix(ind) = j
ind = ind + 1
end do
do j=110,130,20
read(i) icl(ind,1:num_cls)
ix(ind) = j
ind = ind + 1
end do
do j=150,almaxtensor, 50
read(i) icl(ind,1:num_cls)
ix(ind) = j
ind = ind+1
end do
ind = ind-1
call InterpCls(ix,icl, T%cl_tensor, ind, almaxtensor,num_cls)
end if
end if
return
100 error = 1
end if
end subroutine ReadModel
subroutine InterpCls(l,iCl, all_Cl, n, almax, ncls)
integer, intent(in) :: n, almax,ncls
real, intent(in) :: iCl(lmax,1:num_cls)
integer l(n),p
real all_Cl(:,:)
integer il,llo,lhi,xi
real xl(n), ddCl(n)
real a0,b0,ho
real inCl(n)
do p =1, ncls
do il=1,n
inCl(il) = iCl(il,p)*l(il)**2
end do
xl = l
call spline_real(xl,inCl,n,ddCl)
llo=1
do il=2,l(n)
xi=il
if ((xi > l(llo+1)).and.(llo < n)) then
llo=llo+1
end if
lhi=llo+1
ho=l(lhi)-l(llo)
a0=(xl(lhi)-xi)/ho
b0=(xi-xl(llo))/ho
all_Cl(il,p) = (a0*inCl(llo)+ b0*inCl(lhi)+((a0**3-a0)* ddCl(llo) &
+(b0**3-b0)*ddCl(lhi))*ho**2/6)/il**2
end do
end do
all_Cl(l(n)+1:almax,:) = 0
end subroutine InterpCls
subroutine ClsFromTheoryData(T, CMB, Cls)
Type(CosmoTheory) T
Type(CMBParams) CMB
real Cls(lmax,num_cls)
Cls(2:lmax,1:num_clsS) =T%cl(2:lmax,1:num_clsS) !CMB%norm(norm_As)*T%cl(2:lmax,1:num_clsS)
if (num_cls>3 .and. num_ClsS==3) Cls(2:lmax,num_cls)=0
if (CMB%norm(norm_amp_ratio) /= 0) then
Cls(2:lmax_tensor,:) = Cls(2:lmax_tensor,:)+ T%cl_tensor(2:lmax_tensor,:) !CMB%norm(norm_As)*CMB%norm(norm_amp_ratio)*T%cl_tensor(2:lmax_tensor,:)
end if
end subroutine ClsFromTheoryData
subroutine WriteTextCls(aname,T, CMB)
Type(CosmoTheory) T
Type(CMBParams) CMB
character (LEN=*), intent(in) :: aname
integer l
real Cls(lmax,num_cls)
call ClsFromTheoryData(T,CMB,Cls)
open(unit = tmp_file_unit, file = aname, form='formatted', status = 'replace')
do l=2, lmax
write (tmp_file_unit,*) l, Cls(l,:)*l*(l+1)/(2*pi)
end do
close(tmp_file_unit)
end subroutine WriteTextCls
function MatterPowerAt(T,kh)
!get matter power spectrum today at kh = k/h by interpolation from stored values
real, intent(in) :: kh
Type(CosmoTheory) T
real MatterPowerAt
real x, d
integer i
x = log(kh/matter_power_minkh) / matter_power_dlnkh
if (x < 0 .or. x >= num_matter_power-1) then
write (*,*) ' k/h out of bounds in MatterPowerAt (',kh,')'
stop
end if
i = int(x)
d = x - i
MatterPowerAt = exp(log(T%matter_power(i+1,1))*(1-d) + log(T%matter_power(i+2,1))*d)
!Just do linear interpolation in logs for now..
!(since we already cublic-spline interpolated to get the stored values)
end function
! LAPTH NEW : new function (replica of MatterPowerAT for the linear spectrum)
function MatterPowerLinearAt(T,kh)
!get matter power spectrum today at kh = k/h by interpolation from stored values
real, intent(in) :: kh
Type(CosmoTheory) T
real MatterPowerLinearAt
real x, d
integer i
x = log(kh/matter_power_minkh) / matter_power_dlnkh
if (x < 0 .or. x >= num_matter_power-1) then
write (*,*) ' k/h out of bounds in MatterPowerAt (',kh,')'
stop
end if
i = int(x)
d = x - i
MatterPowerLinearAt = exp(log(T%matter_power_linear(i+1,1))*(1-d) + log(T%matter_power_linear(i+2,1))*d)
!Just do linear interpolation in logs for now..
!(since we already cublic-spline interpolated to get the stored values)
end function
! END LAPTH NEW
function MatterPowerAt_Z(T,kh,z)
!get matter power spectrum at z at kh = k/h by interpolation from stored values
real, intent(in) :: kh
Type(CosmoTheory) T
real MatterPowerAt_Z
real x, d, z, y, dz, mup, mdn
real matter_power_dlnz
integer i, iz
matter_power_dlnz = log(matter_power_maxz+1) / (matter_power_lnzsteps -1 + 1e-13)
y = log(1.+ z) / matter_power_dlnz
if (z > matter_power_maxz ) then
write (*,*) ' z out of bounds in MatterPowerAt_Z (',z,')'
stop
end if
x = log(kh/matter_power_minkh) / matter_power_dlnkh
if (x < 0 .or. x >= num_matter_power-1) then
write (*,*) ' k/h out of bounds in MatterPowerAt_Z (',kh,')'
stop
end if
iz = int(y)
dz = y - iz
i = int(x)
d = x - i
mup = log(T%matter_power(i+1,iz+2))*(1-d) + log(T%matter_power(i+2,iz+2))*d
mdn = log(T%matter_power(i+1,iz+1))*(1-d) + log(T%matter_power(i+2,iz+1))*d
MatterPowerAt_Z = exp(mdn*(1-dz) + mup*dz)
end function MatterPowerAt_Z
end module cmbtypes