diff --git a/README b/README new file mode 100644 index 0000000..ec492a8 --- /dev/null +++ b/README @@ -0,0 +1,70 @@ +********************************************************************** +Add-on switches to CAMB and cosmomc: +- compute all multipoles +- remove interpolation +- new fine tuned accuracy parameters +- use external files for the primordial spectra +- second argument in mpirun sets the file suffix +- mean likelihoods plotted in 3D as translucent surfaces +********************************************************************** + +-These patches modify the cosmomc or CAMB codes and therefore may have +bugs that are not present in the originals. The authors of the +original versions should not be bothered with problems resulting in +the use of these patches. Feel free to insult me, or at least contact +me, in case of bugs! + +-Apply the patch in the cosmomc/ directory: +patch -p1 < purpose_MMYY.patch + +*Modifications are traced by the "addon" in the source files. + +*Please note that the genuine Makefiles have been heavily changed and +may require edition to work on your system. They assume the existence +of a directory "cosmomc/WMAP" which points to the WMAP team likelihood +code directory. The WMAP data are also assumed to reside in the +"cosmomc/data" directory. + + +*Modified files and settings + +- camb/modules.f90 +New switches "DoClsInterpol" and "ComputeEveryl". + +- camb/power_tilt.f90 +Can read from external files the scalar and tensor power spectra. Uses +a spline interpolation in ln(k). + +- camb/inidriver.F90 +New keys for the above material. + +- camb/params.ini +New keys for the above material. + +-camb/subroutines.f90 +Added subroutine splinterpol + +- cosmomc/CMB_Cls_simple.f90 +Added screen output to check camb accuracy default values. New +switches lSampleLevel, lAccuracyLevel, CambComputeEveryl to control +the camb accuracy from cosmomc. + +- cosmomc/GetDist.f90 +3d transparent plots of the mean likelihoods. Try colormap('copper'). + +- cosmomc/driver.F90 +Second instance in mpirun sets the file name suffix. + +-cosmomc/MCMC.f90 +Add MpiRank output in the logs + +- cosmomc/settings.f90 +Define lAccuracyLevel,lSampleLevel, CambComputeEveryl, instance_shift. + +- cosmomc/params.ini +New keys for the above material. + +- cosmomc/distparams.ini +New keys for the above material. + + diff --git a/camb/Makefile b/camb/Makefile index 5c09a67..1aa151f 100644 --- a/camb/Makefile +++ b/camb/Makefile @@ -1,83 +1,115 @@ -#CAMB Makefile +# >>> DESIGNED FOR GMAKE <<< +#CAMB System unified Makefile + +ext=$(shell uname | cut -c1-3) +ifeq ($(ext),IRI) +F90C=f90 +FFLAGS = -Ofast -mp -n32 -LANG:recursive=ON +LFLAGS= +endif + +ifeq ($(ext),Lin) +F90C=gfortran +FFLAGS= -O -fopenmp #-DMPI +LFLAGS= +endif -#Set FISHER=Y to compile bispectrum fisher matrix code -FISHER= +ifeq ($(ext),OSF) +F90C=f90 +FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1 +LFLAGS= +endif -#Edit for your compiler -#Note there are many old ifc versions, some of which behave oddly +ifeq ($(ext),Sun) +F90C=f90 +FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none +LFLAGS= +endif -#Intel , -openmp toggles mutli-processor: -#note version 10.0 gives wrong result for lensed when compiled with -openmp [fixed in 10.1] -F90C = ifort -FFLAGS = -openmp -fast -W0 -WB -fpp2 -vec_report0 -ifneq ($(FISHER),) -FFLAGS += -mkl +ifeq ($(ext),AIX) +F90C=xlf90_r +FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90 +LFLAGS= endif -#Gfortran compiler: -#The options here work in v4.5, delete from RHS in earlier versions (15% slower) -#if pre v4.3 add -D__GFORTRAN__ -#With v4.6+ try -Ofast -march=native -fopenmp -#On my machine v4.5 is about 20% slower than ifort -#F90C = gfortran -#FFLAGS = -O3 -fopenmp -ffast-math -march=native -funroll-loops - - -#Old Intel ifc, add -openmp for multi-processor (some have bugs): -#F90C = ifc -#FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2 -#some systems can can also add e.g. -tpp7 -xW - -#G95 compiler -#F90C = g95 -#FFLAGS = -O2 - -#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems. -#F90C = f90 -#FFLAGS = -Ofast -mp - -#Digital/Compaq fortran, -omp toggles multi-processor -#F90C = f90 -#FFLAGS = -omp -O4 -arch host -math_library fast -tune host -fpe1 - -#Absoft ProFortran, single processor: -#F90C = f95 -#FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" - -#NAGF95, single processor: -#F90C = f95 -#FFLAGS = -DNAGF95 -O3 - -#PGF90 -#F90C = pgf90 -#FFLAGS = -O2 -DESCAPEBACKSLASH -Mpreprocess - -#Sun V880 -#F90C = mpf90 -#FFLAGS = -O4 -openmp -ftrap=%none -dalign -DMPI - -#Sun parallel enterprise: -#F90C = f95 -#FFLAGS = -O2 -xarch=native64 -openmp -ftrap=%none -#try removing -openmp if get bus errors. -03, -04 etc are dodgy. - -#IBM XL Fortran, multi-processor (run gmake) -#F90C = xlf90_r -#FFLAGS = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3 - -#Settings for building camb_fits -#Location of FITSIO and name of library -FITSDIR = /home/cpac/cpac-tools/lib -FITSLIB = cfitsio -#Location of HEALPIX for building camb_fits -HEALPIXDIR = /home/cpac/cpac-tools/healpix - -ifneq ($(FISHER),) + +#Files containing evolution equations initial power spectrum module +EQUATIONS = equations +POWERSPECTRUM = power_tilt +REIONIZATION = reionization +RECOMBINATION = recfast +BISPECTRUM = SeparableBispectrum +DENABLE_FISHER = + +#Module doing non-linear scaling +NONLINEAR = halofit + +#Driver program +DRIVER = inidriver.F90 +#DRIVER = sigma8.f90 +#DRIVER = tester.f90 + +CAMBLIB = libcamb.a + + +ifneq ($(DENABLE_FISHER),) FFLAGS += -DFISHER +LFLAGS += -llapack -lblas EXTCAMBFILES = Matrix_utils.o else EXTCAMBFILES = endif -include ./Makefile_main +HEALPIXDIR= + +#Shouldn't need to change anything else... + +F90FLAGS = $(FFLAGS) +HEALPIXLD = -L$(HEALPIXDIR)/lib -lhealpix + +CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o $(REIONIZATION).o modules.o \ + bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o + +default: camb.$(ext) + +all: camb.$(ext) $(CAMBLIB) + +subroutines.o: constants.o utils.o +$(POWERSPECTRUM).o: subroutines.o inifile.o +$(RECOMBINATION).o: subroutines.o inifile.o +$(REIONIZATION).o: constants.o inifile.o +modules.o: $(REIONIZATION).o $(POWERSPECTRUM).o $(RECOMBINATION).o +bessels.o: modules.o +$(EQUATIONS).o: bessels.o +$(NONLINEAR).o: modules.o +lensing.o: bessels.o +$(BISPECTRUM).o: lensing.o modules.o +cmbmain.o: lensing.o $(NONLINEAR).o $(EQUATIONS).o $(BISPECTRUM).o +camb.o: cmbmain.o +Matrix_utils.o: utils.o + + +camb.$(ext): $(CAMBOBJ) $(DRIVER) + $(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) $(LFLAGS) -o $@ + +$(CAMBLIB): $(CAMBOBJ) + ar -r $@ $? + +camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER) + $(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include tils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) $(LFLAGS) -DWRITE_FITS -o $@ + +%.o: %.f90 + $(F90C) $(F90FLAGS) -c $*.f90 + +utils.o: + $(F90C) $(F90FLAGS) -c utils.F90 + +$(BISPECTRUM).o: + $(F90C) $(F90FLAGS) -c $(BISPECTRUM).F90 + +Matrix_utils.o: + $(F90C) $(F90FLAGS) -c Matrix_utils.F90 + +clean: + -rm -f *.o *.a *.d core *.mod diff --git a/camb/inidriver.F90 b/camb/inidriver.F90 index b355761..81eb4ca 100644 --- a/camb/inidriver.F90 +++ b/camb/inidriver.F90 @@ -274,6 +274,16 @@ else lSampleBoost = Ini_Read_Double('l_sample_boost',lSampleBoost) end if +!addon + DoClsInterpol = Ini_Read_Logical('do_cls_interpol',.false.) + ComputeEveryl = Ini_Read_Logical('compute_every_l',.false.) + + if (FeedbackLevel > 0) then + write (*,*)'do_cls_interpol is',DoClsInterpol + write (*,*)'compute_every_l is',ComputeEveryl + endif +!end addon + if (outroot /= '') then if (InputFile /= trim(outroot) //'params.ini') then call Ini_SaveReadValues(trim(outroot) //'params.ini',1) @@ -332,6 +342,10 @@ call CAMB_cleanup + if (P%InitPower%WantExternalScalar.or.P%InitPower%WantExternalTensor) & + call FreeExternalData(P%InitPower) + + end program driver diff --git a/camb/modules.f90 b/camb/modules.f90 index 6053b4c..b24304f 100644 --- a/camb/modules.f90 +++ b/camb/modules.f90 @@ -212,6 +212,10 @@ ! lmax is max possible number of l's evaluated integer, parameter :: lmax_arr = l0max +!addon + logical :: DoClsInterpol=.true. + logical :: ComputeEveryl=.true. +!end addon character(LEN=1024) :: highL_unlensed_cl_template = 'HighLExtrapTemplate_lenspotentialCls.dat' !fiducial high-accuracy high-L C_L used for making small cosmology-independent numerical corrections @@ -581,6 +585,24 @@ integer lind, lvar, step,top,bot,ls(lmax_arr) real(dl) AScale +!addon + if (ComputeEveryl) then + write(*,*) 'Warning: computing every l' + if (max_l.gt.lmax_arr) then + write(*,*)'lmax_arr maxl',lmax_arr,max_l + stop + endif + lind=0 + do lvar=2, max_l + lind=lind+1 + ls(lind)=lvar + end do + + lSet%l0=lind + lSet%l = ls(1:lind) + + else +!end addon Ascale=scale/lSampleBoost if (lSampleBoost >=50) then @@ -721,6 +743,10 @@ lSet%l0=lind lSet%l(1:lind) = ls(1:lind) +!addon + endif +!end addon + end subroutine initlval subroutine InterpolateClArr(lSet,iCl, all_Cl, max_ind) @@ -738,7 +764,11 @@ if (max_ind > lSet%l0) stop 'Wrong max_ind in InterpolateClArr' xl = real(lSet%l(1:lSet%l0),dl) - call spline(xl,iCL(1),max_ind,cllo,clhi,ddCl(1)) + + call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1)) + +!addon + if (DoClsInterpol) then llo=1 do il=lmin,lSet%l(max_ind) @@ -755,7 +785,23 @@ +(b0**3-b0)*ddCl(lhi))*ho**2/6 end do - + + else + + llo=1 + do il=2,lSet%l(max_ind) + xi=il + if ((xi > lSet%l(llo+1)).and.(llo < max_ind)) then + llo=llo+1 + end if + + all_Cl(il) = iCl(llo+1) + + end do + all_Cl(2)=iCl(1) + endif +!end addon + end subroutine InterpolateClArr subroutine InterpolateClArrTemplated(lSet,iCl, all_Cl, max_ind, template_index) diff --git a/camb/params.ini b/camb/params.ini index cf5b821..426b751 100644 --- a/camb/params.ini +++ b/camb/params.ini @@ -67,6 +67,13 @@ nu_mass_fractions = 1 #Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}. initial_power_num = 1 + +#read the power spectra form ascii files in the form (Mpc^-1): k P(k) +use_external_scalar = F +use_external_tensor = F +scalar_input_file(1) = scalar_power.dat +tensor_input_file(1) = tensor_power.dat + pivot_scalar = 0.05 pivot_tensor = 0.05 scalar_amp(1) = 2.1e-9 @@ -218,3 +225,9 @@ l_accuracy_boost = 1 #interpolation errors may be up to 3% #Decrease to speed up non-flat models a bit l_sample_boost = 1 + +#Use spline interpolation for the Cls +do_cls_interpol = F + +# Remove undersampling of the Cls: if true lSampleBoost has no effect +compute_every_l = T diff --git a/camb/power_tilt.f90 b/camb/power_tilt.f90 index c23155f..53066d5 100644 --- a/camb/power_tilt.f90 +++ b/camb/power_tilt.f90 @@ -50,9 +50,28 @@ real(dl) rat(nnmax) !ratio of scalar to tensor initial power spectrum amplitudes real(dl) k_0_scalar, k_0_tensor real(dl) ScalarPowerAmp(nnmax) +!addon + logical :: WantExternalScalar + logical :: WantExternalTensor + character(80) :: ScalarPowerFileName(nnmax),TensorPowerFileName(nnmax) +!end addon end Type InitialPowerParams +!addon + type ExternalData +!only pointers may be defered shape in derived type + real(dl) :: ScalarMomentumMin, ScalarMomentumMax + real(dl), dimension(:), pointer :: PtrScalarPower, PtrScalarLnMomentum + real(dl), dimension(:), pointer :: PtrScalarPower2 + real(dl) :: TensorMomentumMin, TensorMomentumMax + real(dl), dimension(:), pointer :: PtrTensorPower, PtrTensorLnMomentum + real(dl), dimension(:), pointer :: PtrTensorPower2 + end type ExternalData + + type(ExternalData), dimension(:), pointer :: PtrExternalData => null() +!end addon + real(dl) curv !Curvature contant, set in InitializePowers Type(InitialPowerParams) :: P @@ -61,6 +80,11 @@ public InitialPowerParams, InitialPower_ReadParams, InitializePowers, ScalarPower, TensorPower public nnmax,Power_Descript, Power_Name, SetDefPowerParams + +!addon + public FreeExternalData +!end addon + ! public contains @@ -76,6 +100,11 @@ AP%k_0_scalar = 0.05 AP%k_0_tensor = 0.05 AP%ScalarPowerAmp = 1 +!addon +! nullify(PtrExternalData) + AP%WantExternalScalar = .false. + AP%WantExternalTensor = .false. +!end addon end subroutine SetDefPowerParams @@ -86,6 +115,14 @@ real(dl) acurv +!addon + integer :: i,ioflag,j + integer :: nPowerUnit=111 + integer :: nData + real(dl), save :: sbuffer1,sbuffer2 + real(dl), parameter :: spl_large=1.e40_dl +!end addon + if (AParamSet%nn > nnmax) then write (*,*) 'To use ',AParamSet%nn,'power spectra you need to increase' write (*,*) 'nnmax in power_tilt.f90, currently ',nnmax @@ -96,9 +133,163 @@ !Write implementation specific code here... +!addon + +!create and point to P%nn ExternalData objects (pointers and reals) + if ((P%WantExternalScalar).or.(P%WantExternalTensor)) then + if (.not.associated(PtrExternalData)) then + allocate(PtrExternalData(P%nn)) + else + write(*,*)'InitializePowers: PtrExternalData already Initialized!' + return + endif + endif + + if (P%WantExternalScalar) then + + do i=1,P%nn + + write(*,*)'using primordial scalar power spectrum: ',P%ScalarPowerFileName(i) + + open(unit=nPowerUnit,file=P%ScalarPowerFileName(i),action='read' & + ,form='formatted',status='old') +!count the number of records + nData=0 + do + read(nPowerUnit,*,iostat=ioflag) + if (ioflag.lt.0) exit + nData=nData + 1 + enddo +!create and point to + allocate(PtrExternalData(i)%PtrScalarPower(nData) & + ,PtrExternalData(i)%PtrScalarLnMomentum(nData)) + + rewind(nPowerUnit) + nData=1 + do + read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2 + if (ioflag.lt.0) exit + PtrExternalData(i)%PtrScalarLnMomentum(nData) = log(sbuffer1) + PtrExternalData(i)%PtrScalarPower(nData) = sbuffer2 + nData=nData+1 + enddo + nData=nData - 1 + close(nPowerUnit,iostat=ioflag) + +!keep trace of the accessible range + PtrExternalData(i)%ScalarMomentumMin = exp(PtrExternalData(i)%PtrScalarLnMomentum(1)) + PtrExternalData(i)%ScalarMomentumMax = exp(PtrExternalData(i)%PtrScalarLnMomentum(nData)) + + +!2nd order derivatives allocation + allocate(PtrExternalData(i)%PtrScalarPower2(nData)) +!splining in log k + call spline(PtrExternalData(i)%PtrScalarLnMomentum(:) & + ,PtrExternalData(i)%PtrScalarPower(:),nData, spl_large & + , spl_large,PtrExternalData(i)%PtrScalarPower2(:)) + + enddo + + endif + + +!the same for tensors + + if (P%WantExternalTensor) then + + do i=1,P%nn + + write(*,*)'using primordial tensor power spectrum: ',P%TensorPowerFileName(i) + + open(unit=nPowerUnit,file=P%TensorPowerFileName(i) & + ,action='read',form='formatted',status='old') +!count the number of records + nData=0 + do + read(nPowerUnit,*,iostat=ioflag) + if (ioflag.lt.0) exit + nData=nData + 1 + enddo +!store the data + allocate(PtrExternalData(i)%PtrTensorPower(nData) & + ,PtrExternalData(i)%PtrTensorLnMomentum(nData)) + + rewind(nPowerUnit) + nData=1 + do + read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2 + if (ioflag.lt.0) exit + PtrExternalData(i)%PtrTensorLnMomentum(nData) = log(sbuffer1) + PtrExternalData(i)%PtrTensorPower(nData) = sbuffer2 + nData=nData+1 + enddo + nData=nData - 1 + close(nPowerUnit) + +!keep trace of the accessible range + PtrExternalData(i)%TensorMomentumMin = exp(PtrExternalData(i)%PtrTensorLnMomentum(1)) + PtrExternalData(i)%TensorMomentumMax = exp(PtrExternalData(i)%PtrTensorLnMomentum(nData)) + + +!2nd order derivatives + allocate(PtrExternalData(i)%PtrTensorPower2(nData)) +!splining in log k + call spline(PtrExternalData(i)%PtrTensorLnMomentum(:) & + ,PtrExternalData(i)%PtrTensorPower(:),nData, spl_large & + , spl_large,PtrExternalData(i)%PtrTensorPower2(:)) + + enddo + + + endif + end subroutine InitializePowers + + subroutine FreeExternalData(AP) + type (InitialPowerParams), intent(in) :: AP + integer :: i + + if (AP%nn.ne.size(PtrExternalData,1)) then + write(*,*)'InitialPower.FreeExternalData : wrong initial power params' + stop + endif + + if (AP%WantExternalScalar) then + do i=1, AP%nn +!free heap arrays and disassociate ExternalData pointers + deallocate(PtrExternalData(i)%PtrScalarPower) + deallocate(PtrExternalData(i)%PtrScalarPower2) + deallocate(PtrExternalData(i)%PtrScalarLnMomentum) +!should still point to ExternalData objects + if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops' + enddo + endif + if (AP%WantExternalTensor) then + do i=1, AP%nn + deallocate(PtrExternalData(i)%PtrTensorPower) + deallocate(PtrExternalData(i)%PtrTensorPower2) + deallocate(PtrExternalData(i)%PtrTensorLnMomentum) + if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops' + enddo + endif +!Free remaining objects and disassociate pointer to ExternalData + if (associated(PtrExternalData)) then + deallocate(PtrExternalData) + PtrExternalData => null() + else + if (AP%WantExternalScalar.or.AP%WantExternalTensor) then +!should never happen + stop 'InitialPower.FreeExternalData : ooops' + endif + endif + + end subroutine FreeExternalData + +!end addon + + function ScalarPower(k,in) !"in" gives the index of the power to return for this k @@ -119,10 +310,30 @@ real(dl) ScalarPower,k, lnrat integer in +!addon + integer :: nData + + if (P%WantExternalScalar) then + nData=size(PtrExternalData(in)%PtrScalarLnMomentum,1) + + if ((k.ge.PtrExternalData(in)%ScalarMomentumMin) & + .and.(k.le.PtrExternalData(in)%ScalarMomentumMax)) then + call splinterpol(PtrExternalData(in)%PtrScalarLnMomentum & + ,PtrExternalData(in)%PtrScalarPower & + ,PtrExternalData(in)%PtrScalarPower2,nData,log(k),ScalarPower) +!confusing ScalarPower = P%ScalarPowerAmp(in)*ScalarPower + else + write(*,*)'InitPower.ScalarPower() : out of range: k= ',k + ScalarPower=0_dl + endif + + else + + lnrat = log(k/P%k_0_scalar) ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2) - -! ScalarPower = ScalarPower * (1 + 0.1*cos( lnrat*30 ) ) + endif +!end addon end function ScalarPower @@ -144,7 +355,30 @@ integer in +!addon + integer :: nData + + + if (P%WantExternalTensor) then + nData=size(PtrExternalData(in)%PtrTensorLnMomentum,1) + if ((k.ge.PtrExternalData(in)%TensorMomentumMin) & + .and.(k.le.PtrExternalData(in)%TensorMomentumMax)) then + call splinterpol(PtrExternalData(in)%PtrTensorLnMomentum & + ,PtrExternalData(in)%PtrTensorPower & + ,PtrExternalData(in)%PtrTensorPower2,nData,log(k),TensorPower) + +!confusing TensorPower = P%rat(in)*P%ScalarPowerAmp(in)*TensorPower + else + write(*,*)'InitPower.TensorPower() : out of range: k= ',k + TensorPower=0_dl + endif + else + TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor)) + + endif +!end addon + if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) @@ -195,17 +429,36 @@ InitPower%k_0_scalar = Ini_Read_Double_File(Ini,'pivot_scalar',InitPower%k_0_scalar) InitPower%k_0_tensor = Ini_Read_Double_File(Ini,'pivot_tensor',InitPower%k_0_tensor) + +!addon + InitPower%WantExternalScalar = Ini_Read_Logical('use_external_scalar',.false.) + InitPower%WantExternalTensor = Ini_Read_Logical('use_external_tensor',.false.) +!end addon InitPower%nn = Ini_Read_Int('initial_power_num') if (InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower' InitPower%rat(:) = 1 do i=1, InitPower%nn - +!addon + if (InitPower%WantExternalScalar) then + InitPower%ScalarPowerFileName(i) & + = Ini_Read_String_Array_File(Ini,'scalar_input_file',i) + else +!end addon InitPower%an(i) = Ini_Read_Double_Array_File(Ini,'scalar_spectral_index', i) InitPower%n_run(i) = Ini_Read_Double_Array_File(Ini,'scalar_nrun',i,0._dl) + endif if (WantTensors) then +!addon + if (InitPower%WantExternalTensor) then + InitPower%TensorPowerFileName(i) & + = Ini_Read_String_Array_File(Ini,'tensor_input_file',i) + else +!end addon InitPower%ant(i) = Ini_Read_Double_Array_File(Ini,'tensor_spectral_index',i) InitPower%rat(i) = Ini_Read_Double_Array_File(Ini,'initial_ratio',i) + endif + end if InitPower%ScalarPowerAmp(i) = Ini_Read_Double_Array_File(Ini,'scalar_amp',i,1._dl) diff --git a/camb/readme.html b/camb/readme.html index f5ada26..cff517b 100644 --- a/camb/readme.html +++ b/camb/readme.html @@ -112,7 +112,7 @@ where all quantities are in the synchronous gauge and evaluated at the requested
halofit.f90
-Implements the NonLinear module, to calculate non linear scalings of the matter power spectrum as a function of redshift. Uses HALOFIT (astro-ph/0207664, code thanks to Robert Smith. Note this is only reliable at the several percent level for standard ΛCDM models with power law initial power spectra. This module can be replaced to use a different non-linear fitting method. +Implements the NonLinear module, to calculate non linear scalings of the matter power spectrum as a function of redshift. Uses HALOFIT (astro-ph/0207664, code thanks to Robert Smith, with tweaks from arXiv:1109.4416. Note this is only reliable at the several percent level for standard ΛCDM models with power law initial power spectra. This module can be replaced to use a different non-linear fitting method.
@@ -563,7 +563,9 @@ Implements the NonLinear module, to calculate non linear scalings of the matter
-Scalar errors should rarely exceed 0.3% for min(2500, L well into the damping tail) at default accuracy setting, and 0.1% for 500<L<2000 with high_accuracy_default=T. Matter power spectrum errors are usually dominated by interpolation in the acoustic oscillations, with about 0.2% accuracy with high_accuracy_default (but much better rms accuracy). +Scalar numerical errors should rarely exceed 0.3% for min(2500, L well into the damping tail) at default accuracy setting, and 0.1% for 500<L<2000 with high_accuracy_default=T. Matter power spectrum errors are usually dominated by interpolation in the acoustic oscillations, with about 0.2% accuracy with high_accuracy_default (but much better rms accuracy). For a detailed study of numerical accuracy as of January 2012 see arXiv:1201.3654. + + See also comparison with CMBFAST. Accuracy of course assumes the model is correct, and is dependent on RECFAST being the correct ionization history. Lensed C_l TT, TE and EE are accurate at the same level (to within the approximation that the lensing potential is linear, or the accuracy of the the HALOFIT non-linear model).Extreme models (e.g. scale > 4, h>1) may give errors of 5% or more. @@ -594,10 +596,15 @@ See also comparison with CMBFAST. Accuracy of course assu
REFERENCES
-Some notes and relevant Maple derivations are given here (see also the Appendix of astro-ph/0406096). The CAMB notes outline the equations and approximations used, and relation to standard synchronous-gauge and Newtonian-gauge variables. +Some notes and relevant Maple derivations are given here (see also the Appendix of astro-ph/0406096). The CAMB notes outline the equations and approximations used, and relation to standard synchronous-gauge and Newtonian-gauge variables; see also arXiv:1201.3654. There is a BibTex file of references (including CosmoMC).
+CMB power spectrum parameter degeneracies in the era of precision cosmology
+Cullan Howlett, Antony Lewis, Alex Hall, Anthony Challinor arXiv:1201.3654. +Efficient computation of CMB anisotropies in closed FRW Models
Antony Lewis, Anthony Challinor and Anthony Lasenby astro-ph/9911177 Ap. J. 538:473-476, 2000. @@ -647,7 +654,8 @@ Antony Lewis, astro-ph/0403583 HALOFITStable clustering, the halo model and nonlinear cosmological power spectra
-Smith, R. E. and others, astro-ph/0207664 +Smith, R. E. and others, astro-ph/0207664. +RECOMBINATION
@@ -696,13 +704,25 @@ The Cosmic Linear Anisotropy Solving System (CLASS) II: Blas, Diego and Lesgourgues, Julien and Tram, Thomas. arXiv:1104.2933
+Massive Neutrinos ++CMB power spectrum parameter degeneracies in the era of precision cosmology +
+Cullan Howlett, Antony Lewis, Alex Hall, Anthony Challinor. +arXiv:1201.3654 ++
Evolution of cosmological dark matter perturbations
+Antony Lewis and Anthony Challinor astro-ph/0203507 +Phys. Rev. D66, 023531 (2002) + +Synchronous gauge theory and non-flat models
Complete treatment of CMB anisotropies in a FRW universe
Wayne Hu, Uros Seljak and Matias Zaldarriaga. Phys. Rev. D57:6, 3290-3301, 1998. astro-ph/9709066. -
WKB approx to hyperspherical Bessel functions @@ -719,10 +739,16 @@ Blas, Diego and Lesgourgues, Julien and Tram, Thomas. astro-ph/9603033 Ap.J. 469:2 437-444, 1996
+ +Integral solution for the microwave background + anisotropies in nonflat universes
+ Matias Zaldarriaga, Uros Seljak, Edmund Bertschinger. + ApJ. 494:491-501, 1998. astro-ph/9704265. + +CMBFAST for spatially closed universes
Uros Seljak and Matias Zaldariaga, astro-ph/9911219 --See also the references on the CMBFAST home page. +