Coverage for src/ufig/plugins/add_lensing.py: 90%

1# Copyright (c) 2014 ETH Zurich, Institute of Astronomy, Claudio Bruderer 

2# <claudio.bruderer@phys.ethz.ch> 

3""" 

4Created on Nov 4, 2014 

5@author: Claudio Bruderer 

6 

7""" 

8 

9import healpy as hp 

10import numpy as np 

11from ivy.plugin.base_plugin import BasePlugin 

12 

13from ufig import io_util 

14 

15from .. import coordinate_util 

16 

17NAME = "add shear" 

18 

19 

20def shear_constant(numgalaxies, par): 

21 """ 

22 Generates a constant shear field at the location of each galaxy 

23 

24 :param numgalaxies: number of galaxies, i.e. number of samples 

25 :param par: ctx.parameters; part of ctx containing parameters 

26 :return: Shear 1-component 

27 :return: Shear 2-component 

28 """ 

29 

30 gamma1 = np.ones(numgalaxies) * par.g1_constant 

31 gamma2 = np.ones(numgalaxies) * par.g2_constant 

32 return gamma1, gamma2 

33 

34 

35def load_shear_skymaps(shear_maps, g1_prefactor, maps_remote_dir=""): 

36 """ 

37 Load the map containing a realization on the whole sphere of a given input angular 

38 power spectrum for the shear 

39 

40 :param shear_maps: File name of the shear maps 

41 :param g1_prefactor: Prefactor to account for a potential flipping of the shear 

42 g1-map 

43 :param maps_remote_dir: Root directory where maps are stored 

44 :return gamma1_map: Realization of a shear-cl, 1-component 

45 :return gamma2_map: Realization of a shear-cl, 2-component 

46 """ 

47 

48 if shear_maps == "zeros": 

49 npix = hp.nside2npix(1) 

50 g1_map = np.zeros(npix, dtype=np.float32) 

51 g2_map = np.zeros(npix, dtype=np.float32) 

52 

53 else: 

54 maps = io_util.load_hpmap(shear_maps, maps_remote_dir) 

55 g1_map = g1_prefactor * maps[1].astype( 

56 np.float32 

57 ) # Minus sign is HEALPix convention 

58 g2_map = maps[2].astype(np.float32) 

59 

60 return g1_map, g2_map 

61 

62 

63def shear_from_sky_map(gamma1_map, gamma2_map, w, x, y): 

64 """ 

65 Given an input gamma1 and gamma2 map, sample it at positions (x, y) to get input 

66 shear 

67 

68 :param gamma1_map: Map containing gamma1 information 

69 :param gamma2_map: Map containing gamma2 information 

70 :param w: wcs-object containing all the relevant wcs-transformation information 

71 :param x: pixel x-coordinate 

72 :param y: pixel y-coordinate 

73 :return gamma1: Array containing the sampled gamma1 values 

74 :return gamma2: Array containing the sampled gamma2 values 

75 """ 

76 

77 theta, phi = coordinate_util.xy2thetaphi(w, x, y) 

78 

79 # gamma1 = hp.get_interp_val(gamma1_map, theta=theta, phi=phi, nest=False) 

80 # gamma2 = hp.get_interp_val(gamma2_map, theta=theta, phi=phi, nest=False) 

81 nside = hp.get_nside(gamma1_map) 

82 

83 pix_indices = hp.ang2pix(nside=nside, theta=theta, phi=phi, nest=False) 

84 gamma1 = gamma1_map[pix_indices] 

85 gamma2 = gamma2_map[pix_indices] 

86 

87 return gamma1, gamma2 

88 

89 

90def add_shear(int_e1, int_e2, gamma1, gamma2): 

91 """ 

92 Function that adds shear and the intrinsic ellipticity in the weak shear regime 

93 (small values of gamma1 and gamma2) for the ellipticity definition 

94 (a**2 - b**2) / (a**2 + b**2) 

95 

96 :param int_e1: Intrinsic ellipticity 1-component 

97 :param int_e2: Intrinsic ellipticity 2-component 

98 :param gamma1: Shear 1-component 

99 :param gamma2: Shear 1-component 

100 :return e1: Effective ellipticity 1-component 

101 :return e2: Effective ellipticity 2-component 

102 """ 

103 

104 int_e1_e2 = int_e1 * int_e2 

105 e1 = int_e1 + 2 * gamma1 * (1 - int_e1**2) - 2 * gamma2 * int_e1_e2 

106 e2 = int_e2 + 2 * gamma2 * (1 - int_e2**2) - 2 * gamma1 * int_e1_e2 

107 return e1, e2 

108 

109 

110class Plugin(BasePlugin): 

111 """ 

112 Shear and magnification are applied to the input galaxy catalog and 

113 combined into an effective ellipticity, magnitude, and size 

114 

115 :param numgalaxies: number of galaxies 

116 :param shear_type: whether a constant or variable shear is added 

117 :param int_e1: Intrinsic ellipticity 1-component 

118 :param int_e2: Intrinsic ellipticity 2-component 

119 :param int_mag: Intrinsic magnitude 

120 :param int_r50: Intrinsic r50-radius 

121 

122 :return gamma1: Shear 1-component at every galaxy location 

123 :return gamma2: Shear 2-component at every galaxy location 

124 :return e1: Effective ellipticity 1-component at every galaxy location 

125 :return e2: Effective ellipticity 1-component at every galaxy location 

126 :return kappa: Kappa at every galaxy location 

127 :return mag: Effective magnitude of every galaxy 

128 :return r50: Effective size of every galaxy 

129 """ 

130 

131 def __str__(self): 

132 return NAME 

133 

134 def __call__(self): 

135 par = self.ctx.parameters 

136 

137 if par.shear_type == "constant": 

138 self.ctx.galaxies.gamma1, self.ctx.galaxies.gamma2 = shear_constant( 

139 self.ctx.numgalaxies, par 

140 ) 

141 elif par.shear_type == "grf_sky": 141 ↛ 174line 141 didn't jump to line 174 because the condition on line 141 was always true

142 gamma1_map, gamma2_map = load_shear_skymaps( 

143 par.shear_maps, par.g1_prefactor, par.maps_remote_dir 

144 ) 

145 try: 

146 w = coordinate_util.tile_in_skycoords( 

147 pixscale=par.pixscale, 

148 ra0=par.ra0, 

149 dec0=par.dec0, 

150 crpix_ra=par.crpix_ra, 

151 crpix_dec=par.crpix_dec, 

152 ) 

153 except ( 

154 AttributeError 

155 ): # TODO: This is only here as transition! Remove at some point 

156 par.crpix_ra = par.size_x / 2 + 0.5 

157 par.crpix_dec = par.size_y / 2 + 0.5 

158 

159 w = coordinate_util.tile_in_skycoords( 

160 pixscale=par.pixscale, 

161 ra0=par.ra0, 

162 dec0=par.dec0, 

163 crpix_ra=par.crpix_ra, 

164 crpix_dec=par.crpix_dec, 

165 ) 

166 

167 self.ctx.galaxies.gamma1, self.ctx.galaxies.gamma2 = shear_from_sky_map( 

168 gamma1_map=gamma1_map, 

169 gamma2_map=gamma2_map, 

170 w=w, 

171 x=self.ctx.galaxies.x, 

172 y=self.ctx.galaxies.y, 

173 ) 

174 elif par.shear_type == "zero": 

175 self.ctx.galaxies.gamma1, self.ctx.galaxies.gamma2 = ( 

176 np.zeros_like(self.ctx.galaxies.int_e1), 

177 np.zeros_like(self.ctx.galaxies.int_e1), 

178 ) 

179 

180 else: 

181 raise ValueError(f"unknown shear_type={par.shear_type}") 

182 

183 self.ctx.galaxies.e1, self.ctx.galaxies.e2 = add_shear( 

184 self.ctx.galaxies.int_e1, 

185 self.ctx.galaxies.int_e2, 

186 self.ctx.galaxies.gamma1, 

187 self.ctx.galaxies.gamma2, 

188 ) 

189 

190 # TODO: implement magnification; placeholder at the moment 

191 self.ctx.galaxies.kappa = np.zeros(self.ctx.numgalaxies, dtype=np.float32) 

192 self.ctx.galaxies.r50 = self.ctx.galaxies.int_r50 * 1.0 

193 self.ctx.galaxies.mag = self.ctx.galaxies.int_mag * 1.0