dustatt_calzleit.py 12.3 KB
 Yannick Roehlly committed Apr 03, 2013 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ``````# -*- coding: utf-8 -*- """ Copyright (C) 2013 Centre de données Astrophysiques de Marseille Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt @author: Yannick Roehlly """ import numpy as np from . import common from ...data import Database def k_calzetti2000(wavelength): """Compute the Calzetti et al. (2000) A(λ)/E(B-V)∗ Given a wavelength grid, this function computes the selective attenuation `````` Yannick Roehlly committed May 13, 2013 20 21 22 `````` A(λ)/E(B-V)∗ using the formula from Calzetti at al. (2000). This formula is given for wavelengths between 120 nm and 2200 nm, but this function makes the computation outside. `````` Yannick Roehlly committed Apr 03, 2013 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 `````` Parameters ---------- wavelength : array of floats Wavelength grid in nm. Returns ------- a numpy array of floats """ wavelength = np.array(wavelength) result = np.zeros(len(wavelength)) # Attenuation between 120 nm and 630 nm `````` Yannick Roehlly committed May 13, 2013 38 `````` mask = (wavelength < 630) `````` Yannick Roehlly committed Apr 03, 2013 39 40 41 42 43 `````` result[mask] = 2.659 * (-2.156 + 1.509e3 / wavelength[mask] - 0.198e6 / wavelength[mask] ** 2 + 0.011e9 / wavelength[mask] ** 3) + 4.05 # Attenuation between 630 nm and 2200 nm `````` Yannick Roehlly committed May 13, 2013 44 `````` mask = (wavelength >= 630) `````` Yannick Roehlly committed Apr 03, 2013 45 46 47 48 49 50 51 52 53 `````` result[mask] = 2.659 * (-1.857 + 1.040e3 / wavelength[mask]) + 4.05 return result def k_leitherer2002(wavelength): """Compute the Leitherer et al. (2002) A(λ)/E(B-V)∗ Given a wavelength grid, this function computes the selective attenuation `````` Yannick Roehlly committed May 13, 2013 54 55 56 `````` A(λ)/E(B-V)∗ using the formula from Leitherer at al. (2002). This formula is given for wavelengths between 91.2 nm and 180 nm, but this function makes the computation outside. `````` Yannick Roehlly committed Apr 03, 2013 57 58 59 60 61 62 63 64 65 66 67 68 `````` Parameters ---------- wavelength : array of floats Wavelength grid in nm. Returns ------- a numpy array of floats """ wavelength = np.array(wavelength) `````` Yannick Roehlly committed May 13, 2013 69 70 71 `````` result = (5.472 + 0.671e3 / wavelength - 9.218e3 / wavelength ** 2 + 2.620e6 / wavelength ** 3) `````` Yannick Roehlly committed Apr 03, 2013 72 73 74 75 `````` return result `````` Yannick Roehlly committed May 13, 2013 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 ``````def uv_bump(wavelength, central_wave, gamma, ebump): """Compute the Lorentzian-like Drude profile. Parameters ---------- wavelength : array of floats Wavelength grid in nm. central_wave : float Central wavelength of the bump in nm. gamma : float Width (FWHM) of the bump in nm. ebump : float Amplitude of the bump. Returns ------- a numpy array of floats """ return (ebump * wavelength ** 2 * gamma ** 2 / ((wavelength ** 2 - central_wave ** 2) ** 2 + wavelength ** 2 * gamma ** 2)) def power_law(wavelength, delta): """Power law 'centered' on 550 nm.. Parameters ---------- wavelength : array of floats The wavelength grid in nm. delta : float The slope of the power law. Returns ------- array of floats """ return (wavelength / 550) ** delta def a_vs_ebv(wavelength, bump_wave, bump_width, bump_ampl, power_slope): """Compute the complete attenuation curve A(λ)/E(B-V)* The Leitherer et al. (2002) formula is used between 91.2 and 150 nm and the Calzetti et al. (2000) formula is used after 150 (we do an extrapolation after 2200 nm). When the attenuation becomes negative, it is kept to 0. This continuum is multiplied by the power law and then the UV bump is added. Parameters ---------- wavelength : array of floats The wavelength grid (in nm) to compute the attenuation curve on. bump_wave : float Central wavelength (in nm) of the UV bump. bump_width : float Width (FWHM, in nm) of the UV bump. bump_ampl : float Amplitude of the UV bump. power_slope : float Slope of the power law. Returns ------- attenuation : array of floats The A(λ)/E(B-V)* attenuation at each wavelength of the grid. """ attenuation = np.zeros(len(wavelength)) # Leitherer et al. mask = (wavelength >= 91.2) & (wavelength < 150) attenuation[mask] = k_leitherer2002(wavelength[mask]) # Calzetti et al. mask = (wavelength >= 150) attenuation[mask] = k_calzetti2000(wavelength[mask]) # We set attenuation to 0 where it becomes negative mask = (attenuation < 0) attenuation[mask] = 0 # Power law attenuation = attenuation * power_law(wavelength, power_slope) # UV bump attenuation = attenuation + uv_bump(wavelength, bump_wave, bump_width, bump_ampl) return attenuation `````` Yannick Roehlly committed Apr 03, 2013 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 ``````class Module(common.SEDCreationModule): """Add CCM dust attenuation based on the Calzetti formula If a contribution name is given in the parameter list, the attenuation will be applied only to the flux of this contribution; else, it will be applied to the whole spectrum. """ parameter_list = { "E_BVs_young": ( "float", "E(B-V)*, the colour excess of the stellar continuum light for " "the young population.", None ), "E_BVs_old_factor": ( "float", "Reduction factor for the E(B-V)* of the old population compared " "to the young one (<1).", None ), "young_contribution_name": ( "string", "Name of the contribution containing the spectrum of the " "young population.", `````` Yannick Roehlly committed Apr 03, 2013 192 `````` "m2005_young" `````` Yannick Roehlly committed Apr 03, 2013 193 194 195 196 197 198 `````` ), "old_contribution_name": ( "string", "Name of the contribution containing the spectrum of the " "old population. If it is set to 'None', only one population " "is considered.", `````` Yannick Roehlly committed Apr 03, 2013 199 `````` "m2005_old" `````` Yannick Roehlly committed Apr 03, 2013 200 `````` ), `````` Yannick Roehlly committed May 13, 2013 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 `````` "uv_bump_wavelength": ( "float", "Central wavelength of the UV bump in nm.", 217.5 ), "uv_bump_width": ( "float", "Width (FWHM) of the UV bump in nm.", None ), "uv_bump_amplitude": ( "float", "Amplitude of the UV bump in nm.", None ), "powerlaw_slope": ( "float", "Slope delta of the power law modifying the attenuation curve.", None ), `````` Yannick Roehlly committed Apr 03, 2013 221 222 223 224 225 226 227 228 229 230 `````` "filters": ( "list of strings", "List of the filters for which the attenuation will be computed.", ['V_B90', 'FUV'] ) } out_parameter_list = { "NAME_E_BVs_young": "E(B-V)*, the colour excess of the stellar " "continuum light for the young population.", `````` Yannick Roehlly committed May 22, 2013 231 232 `````` "NAME_E_BVs_old": "E(B-V)*, the colour excess of the stellar " "continuum light for the old population.", `````` Yannick Roehlly committed Apr 03, 2013 233 234 235 236 237 238 239 `````` "NAME_attenuation_young": "Amount of luminosity attenuated from the " "young population in W.", "NAME_E_BVs_old_factor": "Ratio of the old population E(B-V)* to the " "young one.", "NAME_attenuation_old": "Amount of luminosity attenuated from the " "old population in W.", "NAME_attenuation": "Total amount of luminosity attenuated in W.", `````` Yannick Roehlly committed May 13, 2013 240 241 242 243 `````` "NAME_uv_bump_wavelength": "Central wavelength of UV bump in nm.", "NAME_uv_bump_width": "Width of the UV bump in nm.", "NAME_uv_bump_amplitude": "Amplitude of the UV bump in nm.", "NAME_powerlaw_slope": "Slope of the power law.", `````` Yannick Roehlly committed Apr 03, 2013 244 245 246 `````` "NAME_FILTER": "Attenuation in the FILTER filter.", } `````` Yannick Roehlly committed Jun 28, 2013 247 248 249 250 251 252 253 254 255 `````` def _init_code(self): """Get the filters from the database""" self.filters = {} base = Database() for filter_name in self.parameters["filters"]: self.filters[filter_name] = base.get_filter(filter_name) base.close() def process(self, sed): `````` Yannick Roehlly committed Apr 03, 2013 256 257 258 259 260 261 262 263 264 `````` """Add the CCM dust attenuation to the SED. Parameters ---------- sed : pcigale.sed.SED object """ # Base name for adding information to the SED. `````` Yannick Roehlly committed May 22, 2013 265 `````` name = self.name or 'dustatt_calzleit_' `````` Yannick Roehlly committed Apr 03, 2013 266 267 `````` wavelength = sed.wavelength_grid `````` Yannick Roehlly committed Jun 28, 2013 268 269 270 271 272 273 274 275 276 `````` ebvs_young = float(self.parameters["E_BVs_young"]) ebvs_old = float(self.parameters["E_BVs_old_factor"]) * ebvs_young young_contrib = self.parameters["young_contribution_name"] old_contrib = self.parameters["old_contribution_name"] uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"]) uv_bump_width = float(self.parameters["uv_bump_wavelength"]) uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"]) powerlaw_slope = float(self.parameters["powerlaw_slope"]) filters = self.filters `````` Yannick Roehlly committed Apr 03, 2013 277 `````` `````` Yannick Roehlly committed Jul 02, 2013 278 `````` # Fλ fluxes (only from continuum) in each filter before attenuation. `````` Yannick Roehlly committed Apr 03, 2013 279 `````` flux_noatt = {} `````` Yannick Roehlly committed Jun 28, 2013 280 `````` for filter_name, filter_ in filters.items(): `````` Yannick Roehlly committed Apr 03, 2013 281 `````` flux_noatt[filter_name] = sed.compute_fnu( `````` Yannick Roehlly committed Jun 28, 2013 282 `````` filter_.trans_table, `````` Yannick Roehlly committed Jul 02, 2013 283 284 `````` filter_.effective_wavelength, add_line_fluxes=False) `````` Yannick Roehlly committed Apr 03, 2013 285 `````` `````` Yannick Roehlly committed May 13, 2013 286 287 288 289 `````` # Compute attenuation curve sel_attenuation = a_vs_ebv(wavelength, uv_bump_wavelength, uv_bump_width, uv_bump_amplitude, powerlaw_slope) `````` Yannick Roehlly committed Apr 03, 2013 290 291 292 `````` # Young population attenuation luminosity = sed.get_lumin_contribution(young_contrib) `````` Yannick Roehlly committed May 13, 2013 293 294 295 `````` attenuated_luminosity = luminosity * 10 ** (ebvs_young * sel_attenuation / -2.5) attenuated_luminosity[wavelength < 91.2] = 0 `````` Yannick Roehlly committed Apr 03, 2013 296 297 298 299 300 `````` attenuation_spectrum = attenuated_luminosity - luminosity # We integrate the amount of luminosity attenuated (-1 because the # spectrum is negative). attenuation_young = -1 * np.trapz(attenuation_spectrum, wavelength) `````` Yannick Roehlly committed Jun 28, 2013 301 `````` sed.add_module(name, self.parameters) `````` Yannick Roehlly committed Apr 03, 2013 302 303 304 305 306 307 308 `````` sed.add_info(name + "_E_BVs_young", ebvs_young) sed.add_info(name + "_attenuation_young", attenuation_young) sed.add_contribution(name + "_young", wavelength, attenuation_spectrum) # Old population (if any) attenuation if old_contrib: luminosity = sed.get_lumin_contribution(old_contrib) `````` Yannick Roehlly committed May 13, 2013 309 310 311 312 `````` attenuated_luminosity = (luminosity * 10 ** (ebvs_old * sel_attenuation / -2.5)) attenuated_luminosity[wavelength < 91.2] = 0 `````` Yannick Roehlly committed Apr 03, 2013 313 314 315 316 317 `````` attenuation_spectrum = attenuated_luminosity - luminosity attenuation_old = -1 * np.trapz(attenuation_spectrum, wavelength) sed.add_info(name + "_E_BVs_old", ebvs_old) sed.add_info(name + "_E_BVs_old_factor", `````` Yannick Roehlly committed Jun 28, 2013 318 `````` self.parameters["E_BVs_old_factor"]) `````` Yannick Roehlly committed Apr 03, 2013 319 320 321 322 323 324 `````` sed.add_info(name + "_attenuation_old", attenuation_old) sed.add_contribution(name + "_old", wavelength, attenuation_spectrum) else: attenuation_old = 0 `````` Yannick Roehlly committed Jul 02, 2013 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 `````` # Attenuation of spectral lines if young_contrib in sed.lines: line_attenuation = a_vs_ebv(sed.lines[young_contrib][0], uv_bump_wavelength, uv_bump_width, uv_bump_amplitude, powerlaw_slope) sed.lines[young_contrib][1] *= 10 ** (ebvs_young * line_attenuation / -2.5) if old_contrib in sed.lines: line_attenuation = a_vs_ebv(sed.lines[old_contrib][0], uv_bump_wavelength, uv_bump_width, uv_bump_amplitude, powerlaw_slope) sed.lines[old_contrib][1] *= 10 ** (ebvs_old * line_attenuation / -2.5) # Total attenuation (we don't take into account the energy attenuated # in the spectral lines) `````` Yannick Roehlly committed Apr 03, 2013 341 342 343 `````` sed.add_info(name + "_attenuation", attenuation_young + attenuation_old) `````` Yannick Roehlly committed Jul 02, 2013 344 `````` # Fλ fluxes (only from continuum) in each filter after attenuation. `````` Yannick Roehlly committed Apr 03, 2013 345 `````` flux_att = {} `````` Yannick Roehlly committed Jun 28, 2013 346 347 348 `````` for filter_name, filter_ in filters.items(): flux_att[filter_name] = sed.compute_fnu( filter_.trans_table, `````` Yannick Roehlly committed Jul 02, 2013 349 350 `````` filter_.effective_wavelength, add_line_fluxes=False) `````` Yannick Roehlly committed Apr 03, 2013 351 352 `````` # Attenuation in each filter `````` Yannick Roehlly committed Jun 28, 2013 353 `````` for filter_name in filters: `````` Yannick Roehlly committed Apr 03, 2013 354 355 356 `````` sed.add_info(name + "_" + filter_name, -2.5 * np.log(flux_att[filter_name] / flux_noatt[filter_name]))``````