Expand the dust attenuation modules to attenuate emission lines.

pcigale/sed_modules/dustatt_2powerlaws.py

@@ -141,6 +141,11 @@ class TwoPowerLawAtt(SedModule):
         self.slope_ISM = float(self.parameters['slope_ISM'])
         self.filter_list = [item.strip() for item in
                             self.parameters["filters"].split("&")]
+        self.Av_ISM = self.Av_BC / self.BC_to_ISM_factor
+        # We cannot compute the attenuation until we know the wavelengths. Yet,
+        # we reserve the object.
+        self.contatt = {}
+        self.lineatt = {}
 
     def process(self, sed):
         """Add the dust attenuation to the SED.
@@ -151,51 +156,51 @@ class TwoPowerLawAtt(SedModule):
 
         """
 
-        def ism_lumin_attenuation(wave):
-            """Compute the luminosity attenuation factor for the ISM"""
-            return 10 ** (self.Av_BC * self.BC_to_ISM_factor *
-                          alambda_av(wave, self.slope_ISM) / -2.5)
-
-        def bc_ism_lumin_attenuation(wave):
-            """Compute the luminosity attenuation factor for ISM + BC"""
-            return 10 ** (self.Av_BC *
-                          alambda_av(wave, self.slope_BC, self.slope_ISM,
-                                     self.BC_to_ISM_factor) / -2.5)
-
-        wavelength = sed.wavelength_grid
+        wl = sed.wavelength_grid
+
+        # Compute the attenuation curves on the continuum wavelength grid
+        if len(self.contatt) == 0:
+            self.contatt['old'] = 10. ** (-.4 * alambda_av(wl, self.slope_ISM) *
+                                          self.Av_ISM)
+            # Emission from the young population is attenuated by both
+            # components
+            self.contatt['young'] = 10. ** (-.4 * alambda_av(wl, self.slope_BC) *
+                                            self.Av_BC) * self.contatt['old']
+
+        # Compute the attenuation curves on the line wavelength grid
+        if len(self.lineatt) == 0:
+            names = [k for k in sed.lines]
+            linewl = np.array([sed.lines[k][0] for k in names])
+            old_curve =  10. ** (-.4 * alambda_av(linewl, self.slope_ISM) *
+                                 self.Av_ISM)
+            young_curve = 10. ** (-.4 * alambda_av(linewl, self.slope_BC) *
+                                  self.Av_BC) * old_curve
+
+            for name, old, young in zip(names, old_curve, young_curve):
+                self.lineatt[name] = (old, young)
 
         # Fλ fluxes in each filter before attenuation.
         flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
 
-        attenuation_total = 0.
+        dust_lumin = 0.
         contribs = [contrib for contrib in sed.contribution_names if
                     'absorption' not in contrib]
 
         for contrib in contribs:
-
             age = contrib.split('.')[-1].split('_')[-1]
-
             luminosity = sed.get_lumin_contribution(contrib)
 
-            # Use the ISM attenuation unless we are dealing with young
-            # populations.
-            if age == 'young':
-                extinction_factor = bc_ism_lumin_attenuation(wavelength)
-            else:
-                extinction_factor = ism_lumin_attenuation(wavelength)
-
-            attenuated_luminosity = luminosity * extinction_factor
-            attenuation_spectrum = attenuated_luminosity - luminosity
-            # We integrate the amount of luminosity attenuated (-1 because the
-            # spectrum is negative).
-            attenuation = -1 * np.trapz(attenuation_spectrum, wavelength)
-            attenuation_total += attenuation
+            attenuation_spectrum = luminosity * (self.contatt[age] - 1.)
+            dust_lumin -= np.trapz(attenuation_spectrum, wl)
 
             sed.add_module(self.name, self.parameters)
-            sed.add_info("attenuation." + contrib, attenuation, True)
-            sed.add_contribution("attenuation." + contrib, wavelength,
+            sed.add_contribution("attenuation." + contrib, wl,
                                  attenuation_spectrum)
 
+        for name, (linewl, old, young) in sed.lines.items():
+            sed.lines[name] = (linewl, old * self.lineatt[name][0],
+                               young * self.lineatt[name][1])
+
         sed.add_info('attenuation.Av_BC', self.Av_BC)
         sed.add_info('attenuation.slope_BC', self.slope_BC)
         sed.add_info('attenuation.BC_to_ISM_factor', self.BC_to_ISM_factor)
@@ -204,10 +209,9 @@ class TwoPowerLawAtt(SedModule):
         # Total attenuation
         if 'dust.luminosity' in sed.info:
             sed.add_info("dust.luminosity",
-                         sed.info["dust.luminosity"]+attenuation_total, True,
-                         True)
+                         sed.info["dust.luminosity"] + dust_lumin, True, True)
         else:
-            sed.add_info("dust.luminosity", attenuation_total, True)
+            sed.add_info("dust.luminosity", dust_lumin, True)
 
         # Fλ fluxes (only in continuum) in each filter after attenuation.
         flux_att = {filt: sed.compute_fnu(filt) for filt in self.filter_list}

pcigale/sed_modules/dustatt_calzleit.py

@@ -252,7 +252,8 @@ class CalzLeit(SedModule):
                             self.parameters["filters"].split("&")]
         # We cannot compute the attenuation until we know the wavelengths. Yet,
         # we reserve the object.
-        self.curve = {}
+        self.contatt = {}
+        self.lineatt = {}
 
     def process(self, sed):
         """Add the CCM dust attenuation to the SED.
@@ -268,12 +269,25 @@ class CalzLeit(SedModule):
         flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
 
         # Compute attenuation curve
-        if len(self.curve) == 0:
+        if len(self.contatt) == 0:
             sel_att = a_vs_ebv(wavelength, self.uv_bump_wavelength,
                                self.uv_bump_width, self.uv_bump_amplitude,
                                self.powerlaw_slope)
             for age in ['young', 'old']:
-                self.curve[age] = 10 ** (-.4 * self.ebvs[age] * sel_att)
+                self.contatt[age] = 10 ** (-.4 * self.ebvs[age] * sel_att)
+
+        # Compute the attenuation curves on the line wavelength grid
+        if len(self.lineatt) == 0:
+            names = [k for k in sed.lines]
+            linewl = np.array([sed.lines[k][0] for k in names])
+            sel_att = a_vs_ebv(linewl, self.uv_bump_wavelength,
+                               self.uv_bump_width, self.uv_bump_amplitude,
+                               self.powerlaw_slope)
+            old_curve = 10 ** (-.4 * self.ebvs['old'] * sel_att)
+            young_curve = 10 ** (-.4 * self.ebvs['young'] * sel_att)
+
+            for name, old, young in zip(names, old_curve, young_curve):
+                self.lineatt[name] = (old, young)
 
         attenuation_total = 0.
         contribs = [contrib for contrib in sed.contribution_names if
@@ -281,7 +295,7 @@ class CalzLeit(SedModule):
         for contrib in contribs:
             age = contrib.split('.')[-1].split('_')[-1]
             luminosity = sed.get_lumin_contribution(contrib)
-            attenuation_spectrum = luminosity * (self.curve[age] - 1.)
+            attenuation_spectrum = luminosity * (self.contatt[age] - 1.)
             # We integrate the amount of luminosity attenuated (-1 because the
             # spectrum is negative).
             attenuation = -.5 * np.dot(np.diff(wavelength),
@@ -295,6 +309,10 @@ class CalzLeit(SedModule):
             sed.add_contribution("attenuation." + contrib, wavelength,
                                  attenuation_spectrum)
 
+        for name, (linewl, old, young) in sed.lines.items():
+            sed.lines[name] = (linewl, old * self.lineatt[name][0],
+                               young * self.lineatt[name][1])
+
         # Total attenuation
         if 'dust.luminosity' in sed.info:
             sed.add_info("dust.luminosity",

pcigale/sed_modules/dustatt_powerlaw.py

@@ -165,7 +165,8 @@ class PowerLawAtt(SedModule):
                             self.parameters["filters"].split("&")]
         # We cannot compute the attenuation until we know the wavelengths. Yet,
         # we reserve the object.
-        self.sel_attenuation = None
+        self.contatt = {}
+        self.lineatt = {}
 
     def process(self, sed):
         """Add the dust attenuation to the SED.
@@ -181,32 +182,46 @@ class PowerLawAtt(SedModule):
         flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
 
         # Compute attenuation curve
-        if self.sel_attenuation is None:
-            self.sel_attenuation = alambda_av(wavelength, self.powerlaw_slope,
-                                              self.uv_bump_wavelength,
-                                              self.uv_bump_width,
-                                              self.uv_bump_amplitude)
-
-        attenuation_total = 0.
+        if len(self.contatt) == 0:
+            att = alambda_av(wavelength, self.powerlaw_slope,
+                             self.uv_bump_wavelength, self.uv_bump_width,
+                             self.uv_bump_amplitude)
+            for age in ['old', 'young']:
+                self.contatt[age] = 10. ** (-.4 * self.av[age] * att)
+
+        # Compute the attenuation curves on the line wavelength grid
+        if len(self.lineatt) == 0:
+            names = [k for k in sed.lines]
+            linewl = np.array([sed.lines[k][0] for k in names])
+            att = alambda_av(linewl, self.powerlaw_slope,
+                             self.uv_bump_wavelength, self.uv_bump_width,
+                             self.uv_bump_amplitude)
+            old_curve = 10 ** (-.4 * self.av['old'] * att)
+            young_curve = 10 ** (-.4 * self.av['young'] * att)
+
+            for name, old, young in zip(names, old_curve, young_curve):
+                print(name, old_curve, young_curve)
+                self.lineatt[name] = (old, young)
+
+        dust_lumin = 0.
         contribs = [contrib for contrib in sed.contribution_names if
                     'absorption' not in contrib]
         for contrib in contribs:
             age = contrib.split('.')[-1].split('_')[-1]
             luminosity = sed.get_lumin_contribution(contrib)
-            attenuated_luminosity = (luminosity * 10. **
-                                     (self.av[age]*self.sel_attenuation/-2.5))
-            attenuation_spectrum = attenuated_luminosity - luminosity
-            # We integrate the amount of luminosity attenuated (-1 because the
-            # spectrum is negative).
-            attenuation = -1 * np.trapz(attenuation_spectrum, wavelength)
-            attenuation_total += attenuation
+
+            attenuation_spectrum = luminosity * (self.contatt[age] - 1.)
+            dust_lumin -= np.trapz(attenuation_spectrum, wavelength)
 
             sed.add_module(self.name, self.parameters)
             sed.add_info("attenuation.Av." + contrib, self.av[age])
-            sed.add_info("attenuation." + contrib, attenuation, True)
             sed.add_contribution("attenuation." + contrib, wavelength,
                                  attenuation_spectrum)
 
+        for name, (linewl, old, young) in sed.lines.items():
+            sed.lines[name] = (linewl, old * self.lineatt[name][0],
+                               young * self.lineatt[name][1])
+
         sed.add_info("attenuation.av_old_factor", self.av_old_factor)
         sed.add_info('attenuation.uv_bump_wavelength', self.uv_bump_wavelength)
         sed.add_info('attenuation.uv_bump_width', self.uv_bump_width)
@@ -216,10 +231,9 @@ class PowerLawAtt(SedModule):
         # Total attenuation
         if 'dust.luminosity' in sed.info:
             sed.add_info("dust.luminosity",
-                         sed.info["dust.luminosity"]+attenuation_total, True,
-                         True)
+                         sed.info["dust.luminosity"] + dust_lumin, True, True)
         else:
-            sed.add_info("dust.luminosity", attenuation_total, True)
+            sed.add_info("dust.luminosity", dust_lumin, True)
 
         # Fλ fluxes (only in continuum) in each filter after attenuation.
         flux_att = {filt: sed.compute_fnu(filt) for filt in self.filter_list}