Coverage for UFalcon/lensing_weights.py : 72%
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
4import numpy as np
5import scipy.integrate as integrate
6from scipy.interpolate import interp1d
8from UFalcon import utils, constants
13class Continuous:
14 """
15 Computes the lensing weights for a continuous, user-defined n(z) distribution.
16 """
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 """
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.")
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)
45 # get the upper bound if necessary
46 if z_lim_up is None:
47 z_lim_up = nz[-1, 0]
49 # get the callable function
50 self.nz_intpt = interp1d(nz[:, 0] - shift_nz, nz[:, 1], bounds_error=False, fill_value=0.0)
52 # points for integration
53 self.lightcone_points = nz[np.logical_and(z_lim_low < nz[:, 0], nz[:, 0] < z_lim_up), 0]
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)
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]
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)
87 return numerator / norm
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)
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]
120 return quad_y(x)
124class Dirac:
125 """
126 Computes the lensing weights for a single-source redshift.
127 """
129 def __init__(self, z_source):
130 """
131 Constructor
132 :param z_source: source redshift
133 """
134 self.z_source = z_source
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 """
145 # source is below the shell --> zero weight
146 if self.z_source <= z_low:
147 w = 0
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')
154 # source is above the shell --> usual weight
155 else:
157 numerator = integrate.quad(self._integrand,
158 z_low,
159 z_up,
160 args=(cosmo,))[0]
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.)
166 w = numerator / norm
168 return w
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)
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
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
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]
210 # Calculate the growth factor today
211 g_norm = 5.0 * OmegaM / 2.0 * integrate.quad(growth, 0, 1)[0]
213 # divide out a
214 g = g / g_norm
216 # critical density today = 3*H^2/(8piG)
217 rho_c = cosmo.critical_density0.to("Msun Mpc^-3").value
219 # Proportionality constant Msun^-1 Mpc^3
220 C1 = 5e-14 / (H0/100.0) ** 2
222 return -IA * rho_c * C1 * OmegaM / g
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
237 """
239 def f(x):
240 return (F_NIA_model(x, IA, cosmo) * nz_intpt(x))
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]
247 return dbl