Implement the nebular continuum emission. This has required some refactoring...

Implement the nebular continuum emission. This has required some refactoring to have some homogeneous naming.

database_builder/__init__.py

@@ -20,7 +20,7 @@ import numpy as np
 from scipy import interpolate
 import scipy.constants as cst
 from pcigale.data import (Database, Filter, SspM2005, SspBC03, AgnFritz2006,
-                          Dale2014, DL2007, Lines)
+                          Dale2014, DL2007, NebularLines, NebularContinuum)
 
 
 def read_bc03_ssp(filename):
@@ -483,37 +483,66 @@ def build_fritz2006(base):
         base.add_fritz2006_agn(agn)
 
 
-def build_lines(base):
-    lines_dir = os.path.join(os.path.dirname(__file__), 'lines/')
-    
+def build_nebular(base):
+    lines_dir = os.path.join(os.path.dirname(__file__), 'nebular/')
+
+    # Number of Lyman continuum photon to normalize the nebular continuum
+    # templates
+    nlyc_continuum = {'0.0001': 2.68786E+53, '0.0004': 2.00964E+53,
+                      '0.004': 1.79593E+53, '0.008': 1.58843E+53,
+                      '0.02': 1.24713E+53, '0.05': 8.46718E+52}
+
     for Z in ['0.0001', '0.0004', '0.004', '0.008', '0.02', '0.05']:
         filename = "{}lines_{}.dat".format(lines_dir, Z)
         print("Importing {}...".format(filename))
         wave, ratio1, ratio2, ratio3 = np.genfromtxt(filename, unpack=True,
                                                      usecols=(0, 3, 7, 11))
-        
+
         # Convert wavelength from Å to nm
         wave *= 0.1
-        
+
         # Convert log(flux) into flux (arbitrary units)
         ratio1 = 10**(ratio1-38.)
         ratio2 = 10**(ratio2-38.)
         ratio3 = 10**(ratio3-38.)
-        
+
         # Normalize all lines to Hβ
-        w = np.where(wave==486.1)
+        w = np.where(wave == 486.1)
         ratio1 = ratio1/ratio1[w]
         ratio2 = ratio2/ratio2[w]
         ratio3 = ratio3/ratio3[w]
-        
-        lines = Lines(np.float(Z), -3., wave, ratio1)
-        base.add_lines(lines)
-        
-        lines = Lines(np.float(Z), -2., wave, ratio2)
-        base.add_lines(lines)
 
-        lines = Lines(np.float(Z), -1., wave, ratio3)
-        base.add_lines(lines)
+        lines = NebularLines(np.float(Z), -3., wave, ratio1)
+        base.add_nebular_lines(lines)
+
+        lines = NebularLines(np.float(Z), -2., wave, ratio2)
+        base.add_nebular_lines(lines)
+
+        lines = NebularLines(np.float(Z), -1., wave, ratio3)
+        base.add_nebular_lines(lines)
+
+        filename = "{}continuum_{}.dat".format(lines_dir, Z)
+        print("Importing {}...".format(filename))
+        wave, cont1, cont2, cont3 = np.genfromtxt(filename, unpack=True,
+                                                  usecols=(0, 3, 7, 11))
+
+        # Convert wavelength from Å to nm
+        wave *= 0.1
+
+        # Normalize flux from erg s¯¹ Hz¯¹ (Msun/yr)¯¹ to W nm¯¹ photon¯¹ s¯¹
+        conv = 1e-7 * cst.c * 1e9 / (wave * wave) / nlyc_continuum[Z]
+        cont1 *= conv
+        cont2 *= conv
+        cont3 *= conv
+
+        cont = NebularContinuum(np.float(Z), -3., wave, cont1)
+        base.add_nebular_continuum(cont)
+
+        cont = NebularContinuum(np.float(Z), -2., wave, cont2)
+        base.add_nebular_continuum(cont)
+
+        cont = NebularContinuum(np.float(Z), -1., wave, cont3)
+        base.add_nebular_continuum(cont)
 
 
 def build_base():
@@ -556,8 +585,8 @@ def build_base():
     print("\nDONE\n")
     print('#' * 78)
     
-    print("8- Build nebular lines\n")
-    build_lines(base)
+    print("8- Importing nebular lines and continuum\n")
+    build_nebular(base)
     print("\nDONE\n")
     print('#' * 78)
 

pcigale/creation_modules/lines.py → pcigale/creation_modules/nebular.py

@@ -12,15 +12,16 @@ from . import CreationModule
 
 class Module(CreationModule):
     """
-    Module computing the nebular emission lines from the ultraviolet
-    to the near-infrared
+    Module computing the nebular emission from the ultraviolet to the
+    near-infrared. It includes both the nebular lines and the nubular
+    continuum. It takes into account the escape fraction and the absorption by
+    dust.
 
     Given the number of Lyman continuum photons, we compute the Hβ line
     luminosity. We then compute the other lines using the
     metallicity-dependent templates that provide the ratio between individual
-    lines and Hβ.
-
-    Information added to the SED: NAME_escape_fraction.
+    lines and Hβ. The nebular continuum is scaled directly from the number of
+    ionizing photons.
 
     """
 
@@ -30,7 +31,7 @@ class Module(CreationModule):
             "Radiation field intensity",
             -2.
         )),
-        ('escape_fraction', (
+        ('f_esc', (
             'float',
             "Fraction of Lyman continuum photons escaping the galaxy",
             0.
@@ -40,7 +41,7 @@ class Module(CreationModule):
             "Fraction of Lyman continuum photons absorbed by dust",
             0.
         )),
-        ('lines_width', (
+        ('nebular_lines_width', (
             'float',
             "Line width in km s¯¹",
             300.
@@ -49,24 +50,45 @@ class Module(CreationModule):
 
     out_parameter_list = OrderedDict([
         ('logU', "Radiation field intensity"),
-        ('escape_fraction', "Fraction of Lyman continuum photons escaping "
+        ('f_esc', "Fraction of Lyman continuum photons escaping "
          "the galaxy"),
         ('f_dust', "Fraction of Lyman continuum photons absorbed by dust"),
-        ('lines_width', "Width of the lines")
+        ('nebular_lines_width', "Width of the nebular lines")
     ])
 
     def _init_code(self):
         """Get the nebular emission lines out of the database and resample
-           them to see the line profile.
+           them to see the line profile. Compute scaling coefficients.
         """
 
-        with Database() as database:
-            self.lines_template = {m:database.
-                                   get_lines(m,self.parameters['logU'])
-                                   for m in database.get_lines_metallicities()
+        fesc = self.parameters['f_esc']
+        fdust = self.parameters['f_dust']
+
+        if fesc < 0. or fesc > 1:
+            raise Exception("Escape fraction must be between 0 and 1")
+
+        if fdust < 0 or fdust > 1:
+            raise Exception("Fraction of lyman photons absorbed by dust must "
+                            "be between 0 and 1")
+
+        if fesc + fdust > 1:
+            raise Exception("Escape fraction+f_dust>1")
+
+        with Database() as db:
+            self.lines_template = {m:
+                                   db.
+                                   get_nebular_lines(m,
+                                                     self.parameters['logU'])
+                                   for m in db.get_nebular_metallicities()
+                                   }
+            self.cont_template = {m:
+                                  db.
+                                  get_nebular_continuum(m,
+                                                        self.parameters['logU'])
+                                  for m in db.get_nebular_metallicities()
                                   }
-        
-        lines_width = self.parameters['lines_width'] * 1e3
+
+        lines_width = self.parameters['nebular_lines_width'] * 1e3
         for lines in self.lines_template.values():
             new_wave = np.array([])
             for line_wave in lines.wave:
@@ -74,32 +96,32 @@ class Module(CreationModule):
                 new_wave = np.concatenate((new_wave, np.linspace(line_wave -
                                                                  3. * width,
                                                                  line_wave
-                                                                 + 3. * 
+                                                                 + 3. *
                                                                  width,
-                                                       19)))
+                                                                 19)))
             new_wave = np.sort(new_wave)
             new_flux = np.zeros_like(new_wave)
             for line_flux, line_wave in zip(lines.ratio, lines.wave):
                 width = line_wave * lines_width / cst.c
                 new_flux += (line_flux * np.exp(- 4. * np.log(2.) *
-                            (new_wave - line_wave) ** 2. /  (width * width)) /
+                            (new_wave - line_wave) ** 2. / (width * width)) /
                             (width * np.sqrt(np.pi / np.log(2.)) / 2.))
             lines.wave = new_wave
             lines.ratio = new_flux
 
         # We compute the conversion coefficient to determine the fluxes using
-        # the formula of Krüger+95: h×ν(Hβ)*α(Hβ,10000K)/αB(10000K)
-        # Osterbrock 1989 gives:
-        # α(Hβ,10000K)=3.03×10¯¹⁴
-        # αB(10000K)2.59×10¯¹³
-        self.conv = (13.598 * cst.electron_volt * (1. / 4. - 1. / 16.) *
-                    3.03e-14 / 2.59e-13 *
-                    (1. - self.parameters['escape_fraction'] -
-                     self.parameters['f_dust'])
-                    )
-        if self.conv < 0.:
-            raise Exception("Escape fraction+f_dust>1")
+        # the formula of Inoue 2011: LHβ=Q(H)*γHβ(10000K)/αβ(10000K)
+        gamma_Hbeta = 1.23e-38  # Inoue 2011, W m³
+        alpha_B = 2.58e-19  # Ferland 1980, m³ s¯¹
 
+        # To take into acount the escape fraction and the fraction of Lyman
+        # continuum photons absorbed by dust we correct by a factor
+        # k=(1-fesc-fdust)/(1+(α1/αβ)*(fesc+fdust))
+        alpha_1 = 1.81e-19  # Ferland 1980, m³ s¯¹
+        k = (1. - fesc - fdust) / (1. + alpha_1 / alpha_B * (fesc + fdust))
+
+        self.conv_line = gamma_Hbeta / alpha_B * k
+        self.conv_cont = k
 
     def process(self, sed):
         """Add the nebular emission lines
@@ -112,19 +134,25 @@ class Module(CreationModule):
         """
 
         # Base name for adding information to the SED.
-        name = self.name or 'lines'
+        name = self.name or 'nebular'
 
-        escape_fraction = self.parameters['escape_fraction']
+        f_esc = self.parameters['f_esc']
         f_dust = self.parameters['f_dust']
         NLy_old = sed.info['ssp_n_ly_old']
         NLy_young = sed.info['ssp_n_ly_young']
         lines = self.lines_template[sed.info['ssp_metallicity']]
+        cont = self.cont_template[sed.info['ssp_metallicity']]
 
         sed.add_module(name, self.parameters)
-        sed.add_info(name + '_escape_fraction', escape_fraction)
+        sed.add_info(name + '_f_esc', f_esc)
         sed.add_info(name + '_f_dust', f_dust)
 
-        sed.add_contribution('lines_old', lines.wave, NLy_old * lines.ratio *
-                             self.conv)
-        sed.add_contribution('lines_young', lines.wave, NLy_young * lines.ratio *
-                             self.conv)
+        sed.add_contribution('nebular_lines_old', lines.wave, lines.ratio *
+                             NLy_old * self.conv_line)
+        sed.add_contribution('nebular_lines_young', lines.wave, lines.ratio *
+                             NLy_young * self.conv_line)
+
+        sed.add_contribution('nebular_continuum_old', cont.wave, cont.lumin *
+                             NLy_old * self.conv_cont)
+        sed.add_contribution('nebular_continuum_young', cont.wave, cont.lumin *
+                             NLy_young * self.conv_cont)

pcigale/data/__init__.py

@@ -27,7 +27,8 @@ from .ir_templates_dh2002 import IrTemplatesDH2002
 from .ir_agn_templates_dale2014 import Dale2014
 from .ir_models_dl2007 import DL2007
 from .agn_fritz2006 import AgnFritz2006
-from .lines import Lines
+from .nebular_cont import NebularContinuum
+from .nebular_lines import NebularLines
 
 
 DATABASE_FILE = pkg_resources.resource_filename(__name__, 'data.db')
@@ -51,7 +52,6 @@ class DatabaseInsertError(Exception):
     """
 
 
-
 class _Filter(BASE):
     """ Storage for filters
     """
@@ -201,21 +201,38 @@ class _Fritz2006AGN(BASE):
         self.luminosity = agn.luminosity
 
 
-class _Lines(BASE):
+class _NebularLines(BASE):
     """Storage for line templates
     """
-    
-    __tablename__ = 'lines'
+
+    __tablename__ = 'nebular_lines'
     metallicity = Column(Float, primary_key=True)
     logU = Column(Float, primary_key=True)
     wave = Column(PickleType)
     ratio = Column(PickleType)
-    
-    def __init__(self, lines):
-        self.metallicity = lines.metallicity
-        self.logU = lines.logU
-        self.wave = lines.wave
-        self.ratio = lines.ratio
+
+    def __init__(self, nebular_lines):
+        self.metallicity = nebular_lines.metallicity
+        self.logU = nebular_lines.logU
+        self.wave = nebular_lines.wave
+        self.ratio = nebular_lines.ratio
+
+
+class _NebularContinuum(BASE):
+    """Storage for nebular continuum templates
+    """
+
+    __tablename__ = 'nebular_continuum'
+    metallicity = Column(Float, primary_key=True)
+    logU = Column(Float, primary_key=True)
+    wave = Column(PickleType)
+    lumin = Column(PickleType)
+
+    def __init__(self, nebular_continuum):
+        self.metallicity = nebular_continuum.metallicity
+        self.logU = nebular_continuum.logU
+        self.wave = nebular_continuum.wave
+        self.lumin = nebular_continuum.lumin
 
 
 class Database(object):
@@ -409,19 +426,34 @@ class Database(object):
         else:
             raise Exception('The database is not writable.')
 
-    def add_lines(self, lines):
+    def add_nebular_lines(self, nebular_lines):
         """
         Add ultraviolet and optical line templates to the database.
         """
         if self.is_writable:
-            self.session.add(_Lines(lines))
+            self.session.add(_NebularLines(nebular_lines))
+            try:
+                self.session.commit()
+            except exc.IntegrityError:
+                self.session.rollback()
+                raise Exception('The line is already in the base')
+        else:
+            raise Exception('The database is not writable')
+
+    def add_nebular_continuum(self, nebular_continuum):
+        """
+        Add nebular continuum templates to the database.
+        """
+        if self.is_writable:
+            self.session.add(_NebularContinuum(nebular_continuum))
             try:
                 self.session.commit()
             except exc.IntegrityError:
                 self.session.rollback()
-                raise StandardError('The line is already in the base')
+                raise Exception('The continuum template is already in the '
+                                'base')
         else:
-            raise StandardError('The database is not writable')
+            raise Exception('The database is not writable')
 
     def get_filter(self, name):
         """
@@ -649,37 +681,57 @@ class Database(object):
                 "The DL2007 model for qpah <{0}>, umin <{1}>, and umax <{2}> "
                 "is not in the database.".format(qpah, umin, umax))
 
-    def get_lines(self, metallicity, logU):
+    def get_nebular_lines(self, metallicity, logU):
         """
         Get the line ratios corresponding to the given set of parameters.
-        
+
         Parameters
         ----------
         metallicity: float
             Gas phase metallicity
-        
         logU: float
             Radiation field intensity
         """
-        result = (self.session.query(_Lines).
-                 filter(_Lines.metallicity == metallicity).
-                 filter(_Lines.logU == logU).
-                 first())
+        result = (self.session.query(_NebularLines).
+                  filter(_NebularLines.metallicity == metallicity).
+                  filter(_NebularLines.logU == logU).
+                  first())
         if result:
-            return Lines(result.metallicity, result.logU, result.wave,
-                         result.ratio)
+            return NebularLines(result.metallicity, result.logU, result.wave,
+                                result.ratio)
         else:
             return None
 
-    def get_lines_metallicities(self):
-        """Get the list of metallicities in the lines database
-        
+    def get_nebular_continuum(self, metallicity, logU):
+        """
+        Get the nebular continuum corresponding to the given set of parameters.
+
+        Parameters
+        ----------
+        metallicity: float
+            Gas phase metallicity
+        logU: float
+            Radiation field intensity
+        """
+        result = (self.session.query(_NebularContinuum).
+                  filter(_NebularContinuum.metallicity == metallicity).
+                  filter(_NebularContinuum.logU == logU).
+                  first())
+        if result:
+            return NebularContinuum(result.metallicity, result.logU,
+                                    result.wave, result.lumin)
+        else:
+            return None
+
+    def get_nebular_metallicities(self):
+        """Get the list of metallicities for the nebular emission.
+
         Returns
         -------
         metallicities: list
-            list of the line templates metallicities
+            list of the nebular emission metallicities
         """
-        result = self.session.query(_Lines.metallicity).all()
+        result = self.session.query(_NebularLines.metallicity).all()
         return [_[0] for _ in result]
 
     def get_filter_list(self):

pcigale/data/nebular_cont.py

+## -*- coding: utf-8 -*-
+# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
+# Copyright (C) 2014 Institute of Astronomy, University of Cambridge
+# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
+# Author: Médéric Boquien
+
+
+class NebularContinuum(object):
+    """Nebular lines templates.
+
+    This class holds the data associated with the line templates
+
+    """
+
+    def __init__(self, metallicity, logU, wave, lumin):
+        """Create a new nebular lines template
+
+        Parameters
+        ----------
+        metallicity: float
+            Gas phase metallicity
+        logU: float
+            Radiation field intensity
+        wave: array
+            Vector of the λ grid used in the templates [nm]
+        lumin: array
+            Luminosity density of the nebular continuum in Fλ
+
+        """
+
+        self.metallicity = metallicity
+        self.logU = logU
+        self.wave = wave
+        self.lumin = lumin

pcigale/data/lines.py → pcigale/data/nebular_lines.py

@@ -5,7 +5,7 @@
 # Author: Médéric Boquien
 
 
-class Lines(object):
+class NebularLines(object):
     """Nebular lines templates.
 
     This class holds the data associated with the line templates