dustatt_calzleit.py 12.6 KB
 Yannick Roehlly committed Apr 03, 2013 1 ``````# -*- coding: utf-8 -*- `````` Yannick Roehlly committed Sep 16, 2013 2 3 4 ``````# 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 `````` Yannick Roehlly committed Apr 03, 2013 5 6 `````` import numpy as np `````` Yannick Roehlly committed Sep 17, 2013 7 ``````from collections import OrderedDict `````` Yannick Roehlly committed Apr 03, 2013 8 9 10 11 12 13 14 15 ``````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 16 17 18 `````` 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 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 `````` 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 34 `````` mask = (wavelength < 630) `````` Yannick Roehlly committed Apr 03, 2013 35 36 37 38 39 `````` 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 40 `````` mask = (wavelength >= 630) `````` Yannick Roehlly committed Apr 03, 2013 41 42 43 44 45 46 47 48 49 `````` 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 50 51 52 `````` 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 53 54 55 56 57 58 59 60 61 62 63 64 `````` 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 65 66 67 `````` result = (5.472 + 0.671e3 / wavelength - 9.218e3 / wavelength ** 2 + 2.620e6 / wavelength ** 3) `````` Yannick Roehlly committed Apr 03, 2013 68 69 70 71 `````` return result `````` Yannick Roehlly committed May 13, 2013 72 73 74 75 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 ``````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 162 163 164 165 166 167 168 169 170 ``````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. """ `````` Yannick Roehlly committed Sep 17, 2013 171 172 `````` parameter_list = OrderedDict([ ("E_BVs_young", ( `````` Yannick Roehlly committed Apr 03, 2013 173 174 175 176 `````` "float", "E(B-V)*, the colour excess of the stellar continuum light for " "the young population.", None `````` Yannick Roehlly committed Sep 17, 2013 177 178 `````` )), ("E_BVs_old_factor", ( `````` Yannick Roehlly committed Apr 03, 2013 179 180 181 182 `````` "float", "Reduction factor for the E(B-V)* of the old population compared " "to the young one (<1).", None `````` Yannick Roehlly committed Sep 17, 2013 183 184 `````` )), ("young_contribution_name", ( `````` Yannick Roehlly committed Apr 03, 2013 185 186 187 `````` "string", "Name of the contribution containing the spectrum of the " "young population.", `````` Yannick Roehlly committed Sep 17, 2013 188 `````` "ssp_young" `````` Yannick Roehlly committed Sep 17, 2013 189 190 `````` )), ("old_contribution_name", ( `````` Yannick Roehlly committed Apr 03, 2013 191 192 193 194 `````` "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 Sep 17, 2013 195 `````` "ssp_old" `````` Yannick Roehlly committed Sep 17, 2013 196 197 `````` )), ("uv_bump_wavelength", ( `````` Yannick Roehlly committed May 13, 2013 198 199 200 `````` "float", "Central wavelength of the UV bump in nm.", 217.5 `````` Yannick Roehlly committed Sep 17, 2013 201 202 `````` )), ("uv_bump_width", ( `````` Yannick Roehlly committed May 13, 2013 203 204 205 `````` "float", "Width (FWHM) of the UV bump in nm.", None `````` Yannick Roehlly committed Sep 17, 2013 206 207 `````` )), ("uv_bump_amplitude", ( `````` Yannick Roehlly committed May 13, 2013 208 209 210 `````` "float", "Amplitude of the UV bump in nm.", None `````` Yannick Roehlly committed Sep 17, 2013 211 212 `````` )), ("powerlaw_slope", ( `````` Yannick Roehlly committed May 13, 2013 213 214 215 `````` "float", "Slope delta of the power law modifying the attenuation curve.", None `````` Yannick Roehlly committed Sep 17, 2013 216 217 `````` )), ("filters", ( `````` Yannick Roehlly committed Sep 17, 2013 218 219 220 221 222 `````` "string", "Filters for which the attenuation will be computed and added to " "the SED information dictionary. You can give several filter " "names separated by a & (don't use commas).", "V_B90 & FUV" `````` Yannick Roehlly committed Sep 17, 2013 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 `````` )) ]) out_parameter_list = OrderedDict([ ("E_BVs_young", "E(B-V)*, the colour excess of the stellar continuum " "light for the young population."), ("E_BVs_old", "E(B-V)*, the colour excess of the stellar " "continuum light for the old population."), ("attenuation_young", "Amount of luminosity attenuated from the " "young population in W."), ("E_BVs_old_factor", "Ratio of the old population E(B-V)* to the " "young one."), ("attenuation_old", "Amount of luminosity attenuated from the " "old population in W."), ("attenuation", "Total amount of luminosity attenuated in W."), ("uv_bump_wavelength", "Central wavelength of UV bump in nm."), ("uv_bump_width", "Width of the UV bump in nm."), ("uv_bump_amplitude", "Amplitude of the UV bump in nm."), ("powerlaw_slope", "Slope of the power law."), `````` Yannick Roehlly committed Sep 17, 2013 242 `````` ("FILTER_attenuation", "Attenuation in the FILTER filter.") `````` Yannick Roehlly committed Sep 17, 2013 243 `````` ]) `````` Yannick Roehlly committed Apr 03, 2013 244 `````` `````` Yannick Roehlly committed Jun 28, 2013 245 246 `````` def _init_code(self): """Get the filters from the database""" `````` Yannick Roehlly committed Sep 17, 2013 247 248 `````` filter_list = [item.strip() for item in self.parameters["filters"].split("&")] `````` Yannick Roehlly committed Jun 28, 2013 249 250 `````` self.filters = {} base = Database() `````` Yannick Roehlly committed Sep 17, 2013 251 `````` for filter_name in filter_list: `````` Yannick Roehlly committed Jun 28, 2013 252 253 254 255 `````` 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 """ wavelength = sed.wavelength_grid `````` Yannick Roehlly committed Jun 28, 2013 265 266 267 268 269 270 271 272 273 `````` 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 274 `````` `````` Yannick Roehlly committed Jul 02, 2013 275 `````` # Fλ fluxes (only from continuum) in each filter before attenuation. `````` Yannick Roehlly committed Apr 03, 2013 276 `````` flux_noatt = {} `````` Yannick Roehlly committed Jun 28, 2013 277 `````` for filter_name, filter_ in filters.items(): `````` Yannick Roehlly committed Apr 03, 2013 278 `````` flux_noatt[filter_name] = sed.compute_fnu( `````` Yannick Roehlly committed Jun 28, 2013 279 `````` filter_.trans_table, `````` Yannick Roehlly committed Jul 02, 2013 280 281 `````` filter_.effective_wavelength, add_line_fluxes=False) `````` Yannick Roehlly committed Apr 03, 2013 282 `````` `````` Yannick Roehlly committed May 13, 2013 283 284 285 286 `````` # 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 287 288 289 `````` # Young population attenuation luminosity = sed.get_lumin_contribution(young_contrib) `````` Yannick Roehlly committed May 13, 2013 290 291 292 `````` attenuated_luminosity = luminosity * 10 ** (ebvs_young * sel_attenuation / -2.5) attenuated_luminosity[wavelength < 91.2] = 0 `````` Yannick Roehlly committed Apr 03, 2013 293 294 295 296 297 `````` 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 Sep 17, 2013 298 299 300 301 302 `````` sed.add_module(self.name, self.parameters) sed.add_info("E_BVs_young" + self.postfix, ebvs_young) sed.add_info("attenuation_young" + self.postfix, attenuation_young) sed.add_contribution("attenuation_young" + self.postfix, wavelength, attenuation_spectrum) `````` Yannick Roehlly committed Apr 03, 2013 303 304 305 306 `````` # Old population (if any) attenuation if old_contrib: luminosity = sed.get_lumin_contribution(old_contrib) `````` Yannick Roehlly committed May 13, 2013 307 308 309 310 `````` attenuated_luminosity = (luminosity * 10 ** (ebvs_old * sel_attenuation / -2.5)) attenuated_luminosity[wavelength < 91.2] = 0 `````` Yannick Roehlly committed Apr 03, 2013 311 312 313 `````` attenuation_spectrum = attenuated_luminosity - luminosity attenuation_old = -1 * np.trapz(attenuation_spectrum, wavelength) `````` Yannick Roehlly committed Sep 17, 2013 314 315 `````` sed.add_info("E_BVs_old" + self.postfix, ebvs_old) sed.add_info("E_BVs_old_factor" + self.postfix, `````` Yannick Roehlly committed Jun 28, 2013 316 `````` self.parameters["E_BVs_old_factor"]) `````` Yannick Roehlly committed Sep 17, 2013 317 318 `````` sed.add_info("attenuation_old" + self.postfix, attenuation_old) sed.add_contribution("attenuation_old" + self.postfix, `````` Yannick Roehlly committed Apr 03, 2013 319 320 321 322 `````` wavelength, attenuation_spectrum) else: attenuation_old = 0 `````` Yannick Roehlly committed Jul 02, 2013 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 `````` # 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 Sep 17, 2013 339 `````` sed.add_info("attenuation" + self.postfix, `````` Yannick Roehlly committed Apr 03, 2013 340 341 `````` attenuation_young + attenuation_old) `````` Yannick Roehlly committed Jul 02, 2013 342 `````` # Fλ fluxes (only from continuum) in each filter after attenuation. `````` Yannick Roehlly committed Apr 03, 2013 343 `````` flux_att = {} `````` Yannick Roehlly committed Jun 28, 2013 344 345 346 `````` for filter_name, filter_ in filters.items(): flux_att[filter_name] = sed.compute_fnu( filter_.trans_table, `````` Yannick Roehlly committed Jul 02, 2013 347 348 `````` filter_.effective_wavelength, add_line_fluxes=False) `````` Yannick Roehlly committed Apr 03, 2013 349 350 `````` # Attenuation in each filter `````` Yannick Roehlly committed Jun 28, 2013 351 `````` for filter_name in filters: `````` Yannick Roehlly committed Sep 17, 2013 352 `````` sed.add_info(filter_name + "_attenuation" + self.postfix, `````` Yannick Roehlly committed Apr 03, 2013 353 354 `````` -2.5 * np.log(flux_att[filter_name] / flux_noatt[filter_name]))``````