Coverage for src/galsbi/ucat/lensing_util.py: 97%

1# Copyright (c) 2017 ETH Zurich, Cosmology Research Group 

2""" 

3Created on Jul 2, 2021 

4@author: Tomasz Kacprzak 

5""" 

6 

7import numpy as np 

8 

9 

10def sigma_to_fwhm(sigma): 

11 return sigma * 2 * np.sqrt(2.0 * np.log(2.0)) 

12 

13 

14def fwhm_to_sigma(fwhm): 

15 return fwhm / (2 * np.sqrt(2.0 * np.log(2.0))) 

16 

17 

18def apply_reduced_shear_to_ellipticities( 

19 int_e1, int_e2, g1, g2, ellipticity_unit="distortion" 

20): 

21 """ 

22 Applies reduced shear to the intrinsic ellipticities, 

23 if unit=='distortion', then following Bartelmann & Schneider 2001, 

24 https://arxiv.org/pdf/astro-ph/9912508.pdf. We use the ellipticity defined 

25 in eq. (4.4), which is sheared according to eq. (4.6), "by interchanging χ and χ(s) 

26 and replacing g by −g". 

27 if unit=='shear', then following ... 

28 

29 :param int_e1: intrinsic ellipticity 1-component 

30 :param int_e2: intrinsic ellipticity 2-component 

31 :param kappa: kappa 

32 :param g1: reduced shear 1-component 

33 :param g2: reduced shear 2-component 

34 :return e1: sheared ellipticity 1-component 

35 :return e2: sheared ellipticity 2-component 

36 """ 

37 

38 # get complex reduced shear 

39 g = g1 + 1j * g2 

40 

41 # compute complex intrinsic ellipticity 

42 int_e = int_e1 + 1j * int_e2 

43 

44 if ellipticity_unit == "distortion": 

45 int_e_conj = np.conjugate(int_e) 

46 

47 # compute complex sheared ellipticity (eq. (4.6) with χ and χ(s) interchanged 

48 # and g --> -g)) 

49 e = (int_e + 2 * g + g**2 * int_e_conj) / ( 

50 1 + np.absolute(g) ** 2 + 2 * np.real(g * int_e_conj) 

51 ) 

52 

53 elif ellipticity_unit == "shear": 

54 e = (int_e + g) / (1 + np.conjugate(g) * int_e) 

55 

56 else: 

57 raise ValueError(f"shape unit {ellipticity_unit} not supported") 

58 

59 return np.real(e), np.imag(e) 

60 

61 

62def apply_shear_to_ellipticities( 

63 int_e1, int_e2, kappa, gamma1, gamma2, ellipticity_unit="distortion" 

64): 

65 """ 

66 Applies shear to the intrinsic ellipticities, following Bartelmann & Schneider 2001, 

67 https://arxiv.org/pdf/astro-ph/9912508.pdf. We use the ellipticity defined in eq. 

68 (4.4), which is sheared according to eq. (4.6), "by interchanging χ and χ(s) and 

69 replacing g by −g". 

70 

71 :param int_e1: intrinsic ellipticity 1-component 

72 :param int_e2: intrinsic ellipticity 2-component 

73 :param kappa: kappa 

74 :param gamma1: shear 1-component 

75 :param gamma2: shear 2-component 

76 :return e1: sheared ellipticity 1-component 

77 :return e2: sheared ellipticity 2-component 

78 """ 

79 

80 # compute complex reduced shear 

81 g1, g2 = shear_to_reduced_shear(gamma1, gamma2, kappa) 

82 

83 # add shear 

84 e1, e2 = apply_reduced_shear_to_ellipticities( 

85 int_e1, int_e2, g1, g2, ellipticity_unit=ellipticity_unit 

86 ) 

87 

88 return e1, e2 

89 

90 

91def shear_to_reduced_shear(gamma1, gamma2, kappa): 

92 """ 

93 Calculate reduced shear 

94 https://arxiv.org/pdf/astro-ph/9407032.pdf1 

95 Eq 2.14 

96 """ 

97 g = (gamma1 + 1j * gamma2) / (1 - kappa) 

98 

99 return g.real, g.imag 

100 

101 

102def distortion_to_shear(e1, e2): 

103 """ 

104 Convert shape in distortion units to shear units 

105 https://arxiv.org/abs/astro-ph/0107431 

106 eq 2-7, 2-8 

107 """ 

108 e = e1 + 1j * e2 

109 g_abs = np.tanh(np.arctanh(np.abs(e)) / 2.0) 

110 g = g_abs * np.exp(1j * np.angle(e)) 

111 return g.real, g.imag 

112 

113 

114def shear_to_distortion(g1, g2): 

115 """ 

116 Convert shape in shear units to distortion units 

117 # https://arxiv.org/abs/astro-ph/0107431 

118 # eq 2-7, 2-8 

119 """ 

120 g = g1 + 1j * g2 

121 e_abs = np.tanh(2 * np.arctanh(np.abs(g))) 

122 e = e_abs * np.exp(1j * np.angle(g)) 

123 return e.real, e.imag 

124 

125 

126def distortion_to_moments(fwhm, e1, e2, xx_out=None, yy_out=None, xy_out=None): 

127 if xx_out is None: 

128 xx_out = np.zeros(len(fwhm), dtype=fwhm.dtype) 

129 if yy_out is None: 

130 yy_out = np.zeros(len(fwhm), dtype=fwhm.dtype) 

131 if xy_out is None: 

132 xy_out = np.zeros(len(fwhm), dtype=fwhm.dtype) 

133 

134 r_sq = fwhm_to_sigma(fwhm) ** 2 

135 xx_out[:] = r_sq * (1 + e1) 

136 yy_out[:] = r_sq * (1 - e1) 

137 xy_out[:] = r_sq * e2 

138 

139 return xx_out, yy_out, xy_out 

140 

141 

142def moments_to_distortion(xx, yy, xy, e1_out=None, e2_out=None, fwhm_out=None): 

143 if e1_out is None: 

144 e1_out = np.zeros(len(xx), dtype=xx.dtype) 

145 if e2_out is None: 

146 e2_out = np.zeros(len(xx), dtype=xx.dtype) 

147 if fwhm_out is None: 

148 fwhm_out = np.zeros(len(xx), dtype=xx.dtype) 

149 

150 xx_yy_sum = xx + yy 

151 e1_out[:] = (xx - yy) / xx_yy_sum 

152 e2_out[:] = 2 * xy / xx_yy_sum 

153 fwhm_out[:] = sigma_to_fwhm(np.sqrt(xx_yy_sum / 2.0)) 

154 return e1_out, e2_out, fwhm_out 

155 

156 

157def moments_to_shear(xx, yy, xy, e1_out=None, e2_out=None, fwhm_out=None): 

158 if e1_out is None: 

159 e1_out = np.zeros(len(xx), dtype=xx.dtype) 

160 if e2_out is None: 

161 e2_out = np.zeros(len(xx), dtype=xx.dtype) 

162 if fwhm_out is None: 

163 fwhm_out = np.zeros(len(xx), dtype=xx.dtype) 

164 

165 zz = xx + yy + 2 * np.sqrt(xx * yy - xy**2) 

166 e1_out[:] = (xx - yy) / zz 

167 e2_out[:] = 2 * xy / zz 

168 fwhm_out[:] = sigma_to_fwhm(np.sqrt((xx + yy) / 2.0)) 

169 

170 return e1_out, e2_out, fwhm_out 

171 

172 

173def shear_to_moments(g1, g2, fwhm, xx_out=None, yy_out=None, xy_out=None): 

174 """ 

175 Convert shape in shear unit to moments 

176 PhD Thesis Tomasz Kacprzak Eqn 2.13 

177 """ 

178 

179 if xx_out is None: 

180 xx_out = np.zeros(len(fwhm), dtype=fwhm.dtype) 

181 if xy_out is None: 

182 xy_out = np.zeros(len(fwhm), dtype=fwhm.dtype) 

183 if yy_out is None: 

184 yy_out = np.zeros(len(fwhm), dtype=fwhm.dtype) 

185 

186 zz = 1 + g1**2 + g2**2 

187 # TODO: check if r2 should be a parameter 

188 r2 = 1 

189 xx_out[:] = r2 * (1 + g1**2 + g2**2 + 2 * g1) / zz 

190 yy_out[:] = r2 * (1 + g1**2 + g2**2 - 2 * g1) / zz 

191 xy_out[:] = r2 * (2 * g2) / zz 

192 return xx_out, yy_out, xy_out 

193 

194 

195def calculate_flux_magnification(kappa, gamma1, gamma2): 

196 # magnification in terms of flux, see Bartelmann & Schneider 2001, 

197 # https://arxiv.org/pdf/astro-ph/9912508.pdf, eq. (3.14) 

198 magnification_flux = 1 / ((1 - kappa) ** 2 - gamma1**2 - gamma2**2) 

199 return magnification_flux 

200 

201 

202def calculate_size_magnification(r, kappa): 

203 magnified_r = r / (1 - kappa) 

204 return magnified_r