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

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

2 

3""" 

4Created on Sept 2021 

5author: Tomasz Kacprzak 

6Code author: Joerg Herbel 

7""" 

8 

9from collections import OrderedDict 

10 

11import numpy as np 

12from cosmic_toolbox import logger 

13 

14from galsbi.ucat import sed_templates_util, spectrum_util 

15from galsbi.ucat.galaxy_population_models.galaxy_luminosity_function import ( 

16 find_closest_ind, 

17) 

18 

19LOGGER = logger.get_logger(__file__) 

20 

21speed_of_light_c_in_micrometer_per_sec = 3 * 1e8 * 1e6 

22 

23TEMPL_INTEGRALS_CACHE = OrderedDict() 

24 

25 

26def flux_to_AB_mag(flux): 

27 mag = -2.5 * np.log10(flux) - 48.6 

28 mag[np.isnan(mag)] = np.inf 

29 return mag 

30 

31 

32def AB_mag_to_flux(mag): 

33 return 10 ** ((mag + 48.6) / (-2.5)) 

34 

35 

36class MagCalculatorDirect: 

37 """ 

38 Computes galaxy magnitudes by integrating redshifted and reddened galaxy spectra 

39 which are given by a sum of five template spectra. See also 

40 docs/jupyter_notebooks/spectra_to_magnitudes.ipynb and 

41 docs/jupyter_notebooks/extinction.ipynb. 

42 

43 We use the following units: 

44 

45 - Wavelength: micrometer 

46 - SED: erg/s/m2/Å 

47 """ 

48 

49 def __init__(self, filters, sed_templates): 

50 self.filters = filters 

51 self.sed_templates = sed_templates 

52 self.extinction_spline = spectrum_util.spline_ext_coeff() 

53 self.n_templates = sed_templates["amp"].shape[0] 

54 

55 def __call__( 

56 self, redshifts, excess_b_v, coeffs, filter_names, return_fluxes=False 

57 ): 

58 # check if iterable 

59 from collections.abc import Iterable 

60 

61 if not isinstance(filter_names, Iterable): 

62 filter_names = [filter_names] 

63 

64 n_obj = redshifts.size 

65 fluxes = {f: np.empty(n_obj) for f in filter_names} 

66 

67 for i in LOGGER.progressbar( 

68 range(n_obj), desc="getting magnitudes", at_level="debug" 

69 ): 

70 lam_obs_in_micrometer = self.sed_templates["lam"] * (1 + redshifts[i]) 

71 

72 spec = spectrum_util.construct_reddened_spectrum( 

73 lam_obs=lam_obs_in_micrometer, 

74 templates_amp=self.sed_templates["amp"], # in erg/s/m2/Å 

75 coeff=coeffs[i], 

76 excess_b_v=excess_b_v[i], 

77 extinction_spline=self.extinction_spline, 

78 ) 

79 

80 for f in filter_names: 

81 filt_amp = np.interp( 

82 lam_obs_in_micrometer, 

83 self.filters[f]["lam"], 

84 self.filters[f]["amp"], 

85 ) 

86 fluxes[f][i] = np.trapz( 

87 lam_obs_in_micrometer * spec * filt_amp, x=lam_obs_in_micrometer 

88 ).item() / ( 

89 speed_of_light_c_in_micrometer_per_sec * self.filters[f]["integ"] 

90 ) 

91 

92 if return_fluxes: 

93 return fluxes 

94 

95 else: 

96 mags = {} 

97 for f in filter_names: 

98 mags[f] = flux_to_AB_mag(fluxes[f]) 

99 mags[f][np.isnan(mags[f])] = np.inf 

100 del fluxes[f] 

101 

102 return mags 

103 

104 

105class MagCalculatorTable: 

106 """ 

107 Computes galaxy magnitudes by looking up pre-computed values of the integrals of our 

108 template spectra as a function of redshift and extinction E(B-V). The integrals need 

109 to be pre-computed for every filter band separately, such that for every filter 

110 band, we have five (number of template spectra) 2-dim. tables of integrals with 

111 redshift on the one and E(B-V) on the other axis. See also 

112 docs/jupyter_notebooks/tabulate_template_integrals.ipynb and 

113 docs/jupyter_notebooks/extinction.ipynb. 

114 """ 

115 

116 def __init__( 

117 self, 

118 filter_names, 

119 filepath_sed_integrals, 

120 reload_cache=False, 

121 copy_to_cwd=False, 

122 ): 

123 # TODO: check if this cache helped at all and if it did, fix the memory leak 

124 """ 

125 if reload_cache: 

126 for key in TEMPL_INTEGRALS_CACHE: 

127 del TEMPL_INTEGRALS_CACHE[key] 

128 sed_templates_util.load_sed_integrals( 

129 filepath_sed_integrals, 

130 filter_names, 

131 crop_negative=True, 

132 sed_templates=TEMPL_INTEGRALS_CACHE, 

133 ) 

134 self.templates_int_table_dict = TEMPL_INTEGRALS_CACHE 

135 self.z_grid = TEMPL_INTEGRALS_CACHE.z_grid 

136 self.excess_b_v_grid = TEMPL_INTEGRALS_CACHE.excess_b_v_grid 

137 

138 """ 

139 self.templates_int_table_dict = sed_templates_util.load_sed_integrals( 

140 filepath_sed_integrals, 

141 filter_names, 

142 crop_negative=True, 

143 sed_templates=None, 

144 copy_to_cwd=copy_to_cwd, 

145 ) 

146 self.z_grid = self.templates_int_table_dict.z_grid 

147 self.excess_b_v_grid = self.templates_int_table_dict.excess_b_v_grid 

148 

149 self.n_templates = self.templates_int_table_dict.n_templates 

150 

151 def __call__( 

152 self, redshifts, excess_b_v, coeffs, filter_names, return_fluxes=False 

153 ): 

154 # check if iterable 

155 from collections.abc import Iterable 

156 

157 if not isinstance(filter_names, Iterable): 

158 filter_names = [filter_names] 

159 

160 # find the values of redshift and extinction in the grid for all objects 

161 z_ind = find_closest_ind(self.z_grid, redshifts) 

162 excess_b_v_ind = find_closest_ind(self.excess_b_v_grid, excess_b_v) 

163 

164 fluxes = {} 

165 

166 for f in filter_names: 

167 templates_int_tables = self.templates_int_table_dict[f] 

168 fluxes[f] = np.zeros(redshifts.size, dtype=np.float64) 

169 

170 # create fluxes from templates 

171 for i in range(coeffs.shape[1]): 

172 fluxes[f] += ( 

173 coeffs[:, i] * templates_int_tables[i][z_ind, excess_b_v_ind] 

174 ) 

175 

176 if return_fluxes: 

177 return fluxes 

178 else: 

179 # convert flux to AB mag 

180 magnitudes = {} 

181 for f in filter_names: 

182 magnitudes[f] = flux_to_AB_mag(fluxes[f]) 

183 del fluxes[f] 

184 return magnitudes