Description |
This page describes the installation of the modules WeakLen (file 'WeakLen.f90') and Angular (file 'AngularCorr.f90'), written by Julien Lesgourgues (LAPTH, Annecy), to be used with the public MCMC code CosmoMC by Sarah Bridle and Antony Lewis.
The modules WeakLen and Angular extend CosmoMC to incorporate weak lensing (cosmic shear) data into the analysis. The modules were written for the paper "A combined analysis of Lyman-α forest, 3D Weak Lensing and WMAP year three data" by Lesgourgues et al., JCAP 11 (2007) 8. They are designed and optimized for the HST COSMOS data released in "COSMOS: 3D weak lensing and the growth of structure" by Massey et al., ApJS 172 (2007) 239. The COSMOS data is included. However, the modules could also be easily adapted to other weak lensing data.
If it is a while since the papers were published, it might be worth checking with the authors for any subsequent updates to code or lensing data.
Download cosmos.tar.gz
. If
you expand it inside your cosmomc/
directory via
"tar -zxvf cosmos.tar.gz
" you will get the files
source/WeakLen.f90 source/AngularCorr.f90 data/cosmos_ctheta.txt data/cosmos_eta.txt data/cosmos_invcov.txt
You may be asked whether you
want to overwrite the existing file
source/WeakLen.f90
. In the current distribution of
CosmoMC (October 2006), this existing file is essentially empty: you can safely
overwrite it. We have not included documentation for the fortan
code, but the files contain embedded explanatory comments.
If you have CosmoMC installed (we refer here to the version of October 2006), you have to make a few simple modifications
in order to make it compatible with the two new source files source/WeakLen.f90
and
source/AngularCorr.f90
:
num_norm
in source/settings.f90
by 3. In the public code, the default value is 3 and should be increased to 6: integer, parameter :: num_norm = 6
norm_A, norm_B, norm_C
to the index of the three new parameters. In source/cmbtypes.f90
, just after the line: integer, parameter :: norm_As=1, norm_amp_ratio=2
integer, parameter :: norm_A=4, norm_B=5, norm_C=6
params.ini
: at the very end of this file, add:
# nuisance param A (COSMOS systematics)
param14 = 1. 0. 2. 0.06 0.06
# nuisance param B (COSMOS systematics)
param15 = 1. 1. 2. 0.05 0.05
# nuisance param C (COSMOS systematics)
param16 = 1. 1. 2. 0.10 0.10
distparams.ini
the number of columns and parameters for the covariance matrix should be increase.
Weaklen
. The few additions remaining to be done are:
source/CMB_Cls_simple.f90
, just after the line: logical :: use_nonlinear = .false.
logical :: Use_WeakLen = .false.
source/Makefile
, in the definition of OBJFILES
, add AngularCorr.o
betwwen cmbdata.o
and WeakLen.o
params.ini
, together with the lines: use_CMB = T
use_min_zre = 0
use_WeakLen = T
use_WeakLen = T
don't forget to include non-linear corrections by setting also nonlinear_pk = T
, otherwise your results will be incorrect.
source/cmbtypes.f90
and change the value of three parameters at the beginning. the recommended values are:
integer, parameter :: num_matter_power = 98 ! default was 74; 98 is fine for COSMOS
real, parameter :: matter_power_minkh = 0.999e-4 ! unchanged
real, parameter :: matter_power_dlnkh = 0.143911568 ! unchanged
real, parameter :: matter_power_maxz = 3. ! default was 0; 3 is fine for COSMOS
integer, parameter :: matter_power_lnzsteps = 10 ! default was 1; 10 is fine for COSMOS
source/CMB_Cls_simple.f90
and change the value of one parameter:P%Transfer%kmax = 100. ! default value was 1.2; 100 is OK for COSMOS
z_and_r(:,:)
where these values will be stored. In source/cmbtypes.f90
, inside the definition of the type Type CosmoTheory
, insert new lines:Type CosmoTheory
...
real z_and_r(3,matter_power_lnzsteps)
end Type CosmoTheory
source/CMB_Cls_simple.f90
, in the line starting by use CAMB, only :
, add another function DeltaTime
to the list:
use CAMB, only : CAMB_GetResults, CAMB_GetAge, CAMBParams, CAMB_SetDefParams,Transfer_GetMatterPower, &
AccuracyBoost, Cl_scalar, Cl_tensor, Cl_lensed, outNone, w_lam, &
CAMBParams_Set,Re_OpticalDepthAtZ, MT, CAMBdata, NonLinear_Pk,&
CAMB_GetTransfers,CAMB_FreeCAMBdata,CAMB_InitCAMBdata, CAMB_TransfersToPowers, &
initial_adiabatic,initial_vector,initial_iso_baryon,initial_iso_neutrino, initial_iso_neutrino_vel, &
DeltaTime
source/CMB_Cls_simple.f90
, in the subroutine subroutine GetCls(CMB,Info, Cls, error)
, at the end of the declaration of variables, add:
real*8 aaa
real*8 dtauda
external dtauda
if (Feedback > 1) write (*,*) 'CAMB done'
if (Use_WeakLen) then
do zix = 1,matter_power_lnzsteps
Info%Theory%z_and_r(1,zix)=P%Transfer%redshifts(zix)
aaa=1.d0/(1.+P%Transfer%redshifts(zix))
Info%Theory%z_and_r(2,zix)=DeltaTime(aaa,1.d0)
Info%Theory%z_and_r(3,zix)=1./(aaa**2*dtauda(aaa))
end do
end if
You are done! In principle you should be able to compile as usual ( cd camb, make clean, make all, cd ../source,
make clean, make cosmomc, cd ..
) and run with ./cosmomc params.ini
. If you start in chatty mode
(Feedback = 2
) the new module will output some useful information if use_WeakLen = T
.
When using weak lensing data (use_WeakLen = T
) it is
crucial to include non-linear corrections by setting also
nonlinear_pk = T
. However, some other data sets require
nonlinear_pk
to be set to F. This is the case for current
Lyman-alpha data modules, which provide constraints on the
linear power spectrum (non linear corrections are accounted for
earlier during the analysis pipeline of Lyman-alpha raw data).
This raises a contradiction which is not difficult to solve.
Our solution is to duplicate various subroutines in camb/module.f90
and source/cmbtypes.f90
in order to define two independent functions:
function MatterPowerAt(T,kh) ! non-linear matter power spectrum
function MatterPowerLinearAt(T,kh) ! linear matter power spectrum
Once this is done, the trick is simply to call MatterPowerLinearAt(T,kh)
in the Lyman-alpha module(s), and MatterPowerAt(T,kh)
in the other modules (weak lensing and galaxy power spectrum). In order to see the list of necessary changes, you can:
camb/module.f90
and
source/cmbtypes.f90
files. Modifications are always indicated by the comment ! LAPTH NEW
. In addition:
source/CMB_Cls_simple.f90
, you should adapt accordingly the following two lines:
call Transfer_GetMatterPower(Info%Transfers%MTrans,&
Info%Theory%matter_power(:,zix),&
Info%Theory%matter_power_linear(:,zix),& !LAPTH NEW: note the new argument
matter_power_lnzsteps-zix+1,&
1,matter_power_minkh, matter_power_dlnkh,num_matter_power)
call Transfer_GetMatterPower(MT,&
Theory%matter_power(:,zix),&
Theory%matter_power_linear(:,zix),& !LAPTH NEW: note the new argument
matter_power_lnzsteps-zix+1,1,matter_power_minkh,&
matter_power_dlnkh,num_matter_power)
If you use the code and publish results, please refer to having used this code and cite the paper JCAP 11 (2007) 8.
Last updated by Richard Massey on 2^{nd} May 2008. |