Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# Copyright (C) 2020 ETH Zurich, Institute for Particle Physics and Astrophysics 

2# Author: Raphael Sgier and Jörg Herbel 

3 

4import numpy as np 

5import scipy.integrate as integrate 

6from scipy.interpolate import interp1d 

7 

8from UFalcon import utils, constants 

9 

10 

11 

12 

13class Continuous: 

14 """ 

15 Computes the lensing weights for a continuous, user-defined n(z) distribution. 

16 """ 

17 

18 def __init__(self, n_of_z, z_lim_low=0, z_lim_up=None, shift_nz=0.0, IA=0.0): 

19 """ 

20 Constructor. 

21 :param n_of_z: either path to file containing n(z), assumed to be a text file readable with numpy.genfromtext 

22 with the first column containing z and the second column containing n(z), or a callable that 

23 is directly a redshift distribution 

24 :param z_lim_low: lower integration limit to use for n(z) normalization, default: 0 

25 :param z_lim_up: upper integration limit to use for n(z) normalization, default: last z-coordinate in n(z) file 

26 :param shift_nz: Can shift the n(z) function by some redshift (intended for easier implementation of photo z bias) 

27 :param IA: Intrinsic Alignment. If not None computes the lensing weights for IA component 

28 (needs to be added to the weights without IA afterwards) 

29 """ 

30 

31 # we handle the redshift dist depending on its type 

32 if callable(n_of_z): 

33 if z_lim_up is None: 

34 raise ValueError("An upper bound of the redshift normalization has to be defined if n_of_z is not a " 

35 "tabulated function.") 

36 

37 self.nz_intpt = n_of_z 

38 # set the integration limit and integration points 

39 self.lightcone_points = None 

40 self.limit = 1000 

41 else: 

42 # read from file 

43 nz = np.genfromtxt(n_of_z) 

44 

45 # get the upper bound if necessary 

46 if z_lim_up is None: 

47 z_lim_up = nz[-1, 0] 

48 

49 # get the callable function 

50 self.nz_intpt = interp1d(nz[:, 0] - shift_nz, nz[:, 1], bounds_error=False, fill_value=0.0) 

51 

52 # points for integration 

53 self.lightcone_points = nz[np.logical_and(z_lim_low < nz[:, 0], nz[:, 0] < z_lim_up), 0] 

54 

55 # check if there are any points remaining for the integration 

56 if len(self.lightcone_points) == 0: 

57 self.lightcone_points = None 

58 self.limit = 1000 

59 else: 

60 self.limit = 10 * len(self.lightcone_points) 

61 

62 self.z_lim_up = z_lim_up 

63 self.z_lim_low = z_lim_low 

64 self.IA = IA 

65 # Normalization 

66 self.nz_norm = integrate.quad(lambda x: self.nz_intpt(x), z_lim_low, self.z_lim_up, 

67 points=self.lightcone_points, limit=self.limit)[0] 

68 

69 def __call__(self, z_low, z_up, cosmo): 

70 """ 

71 Computes the lensing weights for the redshift interval [z_low, z_up]. 

72 :param z_low: lower end of the redshift interval 

73 :param z_up: upper end of the redshift interval 

74 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

75 :return: lensing weight 

76 """ 

77 norm = utils.dimensionless_comoving_distance(z_low, z_up, cosmo) * self.nz_norm 

78 norm *= (utils.dimensionless_comoving_distance(0., (z_low + z_up)/2., cosmo) ** 2.) 

79 if abs(self.IA - 0.0) < 1e-10: 

80 # lensing weights without IA 

81 numerator = integrate.quad(self._integrand_1d, z_low, z_up, args=(cosmo,))[0] 

82 else: 

83 # lensing weights for IA 

84 numerator = (2.0/(3.0*cosmo.Om0)) * \ 

85 w_IA(self.IA, z_low, z_up, cosmo, self.nz_intpt, points=self.lightcone_points) 

86 

87 return numerator / norm 

88 

89 def _integrand_2d(self, y, x, cosmo): 

90 """ 

91 The 2d integrant of the continous lensing weights 

92 :param y: redhsift that goes into the n(z) 

93 :param x: redshift for the Dirac part 

94 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

95 :return: the 2d integrand function 

96 """ 

97 return self.nz_intpt(y) * \ 

98 utils.dimensionless_comoving_distance(0, x, cosmo) * \ 

99 utils.dimensionless_comoving_distance(x, y, cosmo) * \ 

100 (1 + x) * \ 

101 cosmo.inv_efunc(x) / \ 

102 utils.dimensionless_comoving_distance(0, y, cosmo) 

103 

104 def _integrand_1d(self, x, cosmo): 

105 """ 

106 Function that integrates out y from the 2d integrand 

107 :param x: at which x (redshfit to eval) 

108 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

109 :return: the 1d integrant at x 

110 """ 

111 if self.lightcone_points is not None: 

112 points = self.lightcone_points[np.logical_and(self.z_lim_low < self.lightcone_points, 

113 self.lightcone_points < self.z_lim_up)] 

114 quad_y = lambda x: integrate.quad(lambda y: self._integrand_2d(y, x, cosmo), x, self.z_lim_up, 

115 limit=self.limit, points=points)[0] 

116 else: 

117 quad_y = lambda x: integrate.quad(lambda y: self._integrand_2d(y, x, cosmo), x, self.z_lim_up, 

118 limit=self.limit)[0] 

119 

120 return quad_y(x) 

121 

122 

123 

124class Dirac: 

125 """ 

126 Computes the lensing weights for a single-source redshift. 

127 """ 

128 

129 def __init__(self, z_source): 

130 """ 

131 Constructor 

132 :param z_source: source redshift 

133 """ 

134 self.z_source = z_source 

135 

136 def __call__(self, z_low, z_up, cosmo): 

137 """ 

138 Computes the lensing weights for the redshift interval [z_low, z_up]. 

139 :param z_low: lower end of the redshift interval 

140 :param z_up: upper end of the redshift interval 

141 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

142 :return: lensing weight 

143 """ 

144 

145 # source is below the shell --> zero weight 

146 if self.z_source <= z_low: 

147 w = 0 

148 

149 # source is inside the shell --> error 

150 elif self.z_source < z_up: 

151 raise NotImplementedError('Attempting to call UFalcon.lensing_weights.Dirac with z_low < z_source < z_up, ' 

152 'this is not implemented') 

153 

154 # source is above the shell --> usual weight 

155 else: 

156 

157 numerator = integrate.quad(self._integrand, 

158 z_low, 

159 z_up, 

160 args=(cosmo,))[0] 

161 

162 norm = utils.dimensionless_comoving_distance(z_low, z_up, cosmo) * \ 

163 utils.dimensionless_comoving_distance(0, self.z_source, cosmo) 

164 norm *= (utils.dimensionless_comoving_distance(0., (z_low + z_up)/2., cosmo) ** 2.) 

165 

166 w = numerator / norm 

167 

168 return w 

169 

170 def _integrand(self, x, cosmo): 

171 return utils.dimensionless_comoving_distance(0, x, cosmo) * \ 

172 utils.dimensionless_comoving_distance(x, self.z_source, cosmo) * \ 

173 (1 + x) * \ 

174 cosmo.inv_efunc(x) 

175 

176 

177def kappa_prefactor(n_pix, n_particles, boxsize, cosmo): 

178 """ 

179 Computes the prefactor to transform from number of particles to convergence, see https://arxiv.org/abs/0807.3651, 

180 eq. (A.1). 

181 :param n_pix: number of healpix pixels used 

182 :param n_particles: number of particles 

183 :param boxsize: size of the box in Gigaparsec 

184 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

185 :return: convergence prefactor 

186 """ 

187 convergence_factor = (3.0 * cosmo.Om0 / 2.0) * \ 

188 (n_pix / (4.0 * np.pi)) * \ 

189 (cosmo.H0.value / constants.c) ** 3 * \ 

190 (boxsize * 1000.0) ** 3 / n_particles 

191 return convergence_factor 

192 

193 

194def F_NIA_model(z, IA, cosmo): 

195 """ 

196 Calculates the NIA kernel used to calculate the IA shell weight 

197 :param z: Redshift where to evaluate 

198 :param IA: Galaxy intrinsic alignments amplitude 

199 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

200 :return: NIA kernel at redshift z 

201 """ 

202 OmegaM = cosmo.Om0 

203 H0 = cosmo.H0.value 

204 

205 # growth factor calculation 

206 growth = lambda a: 1.0 / (a * cosmo.H(1.0 / a - 1.0).value)**3.0 

207 a = 1.0 / (1.0 + z) 

208 g = 5.0 * OmegaM / 2.0 * cosmo.efunc(z) * integrate.quad(growth, 0, a)[0] 

209 

210 # Calculate the growth factor today 

211 g_norm = 5.0 * OmegaM / 2.0 * integrate.quad(growth, 0, 1)[0] 

212 

213 # divide out a 

214 g = g / g_norm 

215 

216 # critical density today = 3*H^2/(8piG) 

217 rho_c = cosmo.critical_density0.to("Msun Mpc^-3").value 

218 

219 # Proportionality constant Msun^-1 Mpc^3 

220 C1 = 5e-14 / (H0/100.0) ** 2 

221 

222 return -IA * rho_c * C1 * OmegaM / g 

223 

224 

225def w_IA(IA, z_low, z_up, cosmo, nz_intpt, points=None): 

226 """ 

227 Calculates the weight per slice for the NIA model given a 

228 distribution of source redshifts n(z). 

229 :param IA: Galaxy Intrinsic alignment amplitude 

230 :param z_low: Lower redshift limit of the shell 

231 :param z_up: Upper redshift limit of the shell 

232 :param cosmo: Astropy.Cosmo instance, controls the cosmology used 

233 :param nz_intpt: nz function 

234 :param points: Points in redshift where integrand is evaluated (used for better numerical integration), can be None 

235 :return: Shell weight for NIA model 

236 

237 """ 

238 

239 def f(x): 

240 return (F_NIA_model(x, IA, cosmo) * nz_intpt(x)) 

241 

242 if points is not None: 

243 dbl = integrate.quad(f, z_low, z_up, points=points[np.logical_and(z_low < points, points < z_up)])[0] 

244 else: 

245 dbl = integrate.quad(f, z_low, z_up)[0] 

246 

247 return dbl