When using different line widths during a single run, ensure that the fluxes...

When using different line widths during a single run, ensure that the fluxes and other quantities are always computed correctly

CHANGELOG.md

@@ -34,6 +34,7 @@
 - The formula from Sawicki et al. (2012) used to compute the χ² in the presence of upper limits was not correct. This led the χ² to depend directly on the absolute value of the upper limit. The formula has been rederived and corrected. (Médéric Boquien & Denis Burgarella)
 - For some reason the wavelengths of the SCUBA 450 μm filter were a factor 10 too small. (Médéric Boquien)
 - Ensure that the computation of the continuum level is correct when determining the equivalent width, in particular when the line width is very narrow. (Médéric Boquien)
+- When using different line widths during a single run, ensure that the fluxes and other quantities are always computed correctly. (Médéric Boquien, special thanks to Genoveva Micheva)
 
 ### Optimised
 - The cache architecture has been simplified, making it somewhat faster. It speeds up the model generation by ~1%. (Médéric Boquien)

pcigale/sed/__init__.py

@@ -270,11 +270,20 @@ class SED(object):
         # (spectrum+filter), we store the resampled filter. That way we only
         # have to resample to spectrum.
         if 'universe.redshift' in self.info:
-            key = (wavelength.size, filter_name,
-                   self.info['universe.redshift'])
+            if 'nebular.lines_width' in self.info:
+                key = (wavelength.size, filter_name,
+                       self.info['nebular.lines_width'],
+                       self.info['universe.redshift'])
+            else:
+                key = (wavelength.size, filter_name,
+                       self.info['universe.redshift'])
             dist = self.info['universe.luminosity_distance']
         else:
-            key = (wavelength.size, filter_name, 0.)
+            if 'nebular.lines_width' in self.info:
+                key = (wavelength.size, filter_name,
+                       self.info['nebular.lines_width'], 0.)
+            else:
+                key = (wavelength.size, filter_name, 0.)
             dist = 10. * parsec
 
         if key in self.cache_filters:

pcigale/sed_modules/restframe_parameters.py

@@ -82,8 +82,13 @@ class RestframeParam(SedModule):
         # Attenuated (observed) UV slopes beta as defined in Calzetti et al.
         # (1994, ApJ 429, 582, Tab. 2) that excludes the 217.5 nm bump
         # wavelength range and other spectral features
-        if wl.size in self.w_calz94:
-            w_calz94 = self.w_calz94[wl.size]
+        if 'nebular.lines_width' in sed.info:
+            key = (wl.size, sed.info['nebular.lines_width'])
+        else:
+            key = (wl.size, )
+
+        if key in self.w_calz94:
+            w_calz94 = self.w_calz94[key]
         else:
             calz_wl = [(126.8, 128.4), (130.9, 131.6), (134.2, 137.1),
                        (140.7, 151.5), (156.2, 158.3), (167.7, 174.0),
@@ -92,7 +97,7 @@ class RestframeParam(SedModule):
             w_calz94 = np.where(np.any([(wl >= wlseg[0]) & (wl <= wlseg[1])
                                         for wlseg in calz_wl], axis=0))
 
-            self.w_calz94[wl.size] = w_calz94
+            self.w_calz94[key] = w_calz94
 
         # We compute the regression directly from the covariance matrix as the
         # numpy/scipy regression routines are quite slow.
@@ -110,27 +115,33 @@ class RestframeParam(SedModule):
         # i.e., ratio of the flux in the red continuum to that in the blue
         # continuum: Blue continuum: 385.0-395.0 nm & red continuum:
         # 400.0-410.0 nm.
-        if wl.size in self.w_D4000blue:
-            w_D4000blue = self.w_D4000blue[wl.size]
-            w_D4000red = self.w_D4000red[wl.size]
+        if 'nebular.lines_width' in sed.info:
+            key = (wl.size, sed.info['nebular.lines_width'])
+        else:
+            key = (wl.size, )
+
+        if key in self.w_D4000blue:
+            w_D4000blue = self.w_D4000blue[key]
+            w_D4000red = self.w_D4000red[key]
         else:
             w_D4000blue = np.where((wl >= 385.0) & (wl <= 395.0))
             w_D4000red = np.where((wl >= 400.0) & (wl <= 410.0))
-            self.w_D4000blue[wl.size] = w_D4000blue
-            self.w_D4000red[wl.size] = w_D4000red
+            self.w_D4000blue[key] = w_D4000blue
+            self.w_D4000red[key] = w_D4000red
 
         return np.mean(lumin[w_D4000red]) / np.mean(lumin[w_D4000blue])
 
     def EW(self, sed):
         wl = sed.wavelength_grid
 
-        if wl.size in self.w_lines:
-            w_lines = self.w_lines[wl.size]
+        key = (wl.size, sed.info['nebular.lines_width'])
+        if key in self.w_lines:
+            w_lines = self.w_lines[key]
         else:
             w_lines = {line: np.where((wl >= line[0]-line[1]) &
                                       (wl <= line[0]+line[1]))
                        for line in self.lines}
-            self.w_lines[wl.size] = w_lines
+            self.w_lines[key] = w_lines
 
         lumin_line = np.sum([sed.get_lumin_contribution(name)
                              for name in sed.contribution_names
@@ -141,7 +152,7 @@ class RestframeParam(SedModule):
         for line in self.lines:
             w_line = w_lines[line]
             wl_line = wl[w_line]
-            key = (wl_line.size, line[0], 0.)
+            key = (wl_line.size, sed.info['nebular.lines_width'], line[0], 0.)
             EW[line] = (flux_trapz(lumin_line[w_line], wl_line, key) /
                         flux_trapz(lumin_cont[w_line], wl_line, key) *
                         (wl_line[-1]-wl_line[0]))