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 :param int_e1: intrinsic ellipticity 1-component 

29 :param int_e2: intrinsic ellipticity 2-component 

30 :param kappa: kappa 

31 :param g1: reduced shear 1-component 

32 :param g2: reduced shear 2-component 

33 :return e1: sheared ellipticity 1-component 

34 :return e2: sheared ellipticity 2-component 

35 """ 

36 

37 # get complex reduced shear 

38 g = g1 + 1j * g2 

39 

40 # compute complex intrinsic ellipticity 

41 int_e = int_e1 + 1j * int_e2 

42 

43 if ellipticity_unit == "distortion": 

44 int_e_conj = np.conjugate(int_e) 

45 

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

47 # and g --> -g)) 

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

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

50 ) 

51 

52 elif ellipticity_unit == "shear": 

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

54 

55 else: 

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

57 

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

59 

60 

61def apply_shear_to_ellipticities( 

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

63): 

64 """ 

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

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

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

68 replacing g by −g". 

69 :param int_e1: intrinsic ellipticity 1-component 

70 :param int_e2: intrinsic ellipticity 2-component 

71 :param kappa: kappa 

72 :param gamma1: shear 1-component 

73 :param gamma2: shear 2-component 

74 :return e1: sheared ellipticity 1-component 

75 :return e2: sheared ellipticity 2-component 

76 """ 

77 

78 # compute complex reduced shear 

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

80 

81 # add shear 

82 e1, e2 = apply_reduced_shear_to_ellipticities( 

83 int_e1, int_e2, g1, g2, ellipticity_unit=ellipticity_unit 

84 ) 

85 

86 return e1, e2 

87 

88 

89def shear_to_reduced_shear(gamma1, gamma2, kappa): 

90 """ 

91 Calculate reduced shear 

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

93 Eq 2.14 

94 """ 

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

96 

97 return g.real, g.imag 

98 

99 

100def distortion_to_shear(e1, e2): 

101 """ 

102 Convert shape in distortion units to shear units 

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

104 eq 2-7, 2-8 

105 """ 

106 e = e1 + 1j * e2 

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

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

109 return g.real, g.imag 

110 

111 

112def shear_to_distortion(g1, g2): 

113 """ 

114 Convert shape in shear units to distortion units 

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

116 # eq 2-7, 2-8 

117 """ 

118 g = g1 + 1j * g2 

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

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

121 return e.real, e.imag 

122 

123 

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

125 if xx_out is None: 

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

127 if yy_out is None: 

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

129 if xy_out is None: 

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

131 

132 r_sq = fwhm_to_sigma(fwhm) ** 2 

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

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

135 xy_out[:] = r_sq * e2 

136 

137 return xx_out, yy_out, xy_out 

138 

139 

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

141 if e1_out is None: 

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

143 if e2_out is None: 

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

145 if fwhm_out is None: 

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

147 

148 xx_yy_sum = xx + yy 

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

150 e2_out[:] = 2 * xy / xx_yy_sum 

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

152 return e1_out, e2_out, fwhm_out 

153 

154 

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

156 if e1_out is None: 

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

158 if e2_out is None: 

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

160 if fwhm_out is None: 

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

162 

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

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

165 e2_out[:] = 2 * xy / zz 

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

167 

168 return e1_out, e2_out, fwhm_out 

169 

170 

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

172 """ 

173 Convert shape in shear unit to moments 

174 PhD Thesis Tomasz Kacprzak Eqn 2.13 

175 """ 

176 

177 if xx_out is None: 

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

179 if xy_out is None: 

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

181 if yy_out is None: 

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

183 

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

185 # TODO: check if r2 should be a parameter 

186 r2 = 1 

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

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

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

190 return xx_out, yy_out, xy_out 

191 

192 

193def calculate_flux_magnification(kappa, gamma1, gamma2): 

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

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

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

197 return magnification_flux 

198 

199 

200def calculate_size_magnification(r, kappa): 

201 magnified_r = r / (1 - kappa) 

202 return magnified_r