Add the Draine and Li (2007) SED module

Add the necessary structures and methods to add and retrieve them from the
database.

database_builder/__init__.py

@@ -15,10 +15,13 @@ import sys
 import os
 sys.path.append(os.path.join(os.path.dirname(__file__), '../'))
 import glob
+import io
 import itertools
 import numpy as np
 from scipy import interpolate
-from pcigale.data import Database, Filter, SspM2005, SspBC03, AgnFritz2006
+import scipy.constants as cst
+from pcigale.data import (Database, Filter, SspM2005, SspBC03, AgnFritz2006,
+                          DL2007)
 
 
 def read_bc03_ssp(filename):
@@ -312,6 +315,78 @@ def build_dh2002(base):
     base.add_dh2002_infrared_templates(data)
 
 
+def build_dl2007(base):
+    dl2007_dir = os.path.join(os.path.dirname(__file__), 'dl2007/')
+
+    qpah = {
+        "00": 0.47,
+        "10": 1.12,
+        "20": 1.77,
+        "30": 2.50,
+        "40": 3.19,
+        "50": 3.90,
+        "60": 4.58
+    }
+
+    umaximum = ["1e3", "1e4", "1e5", "1e6"]
+    uminimum = ["0.10", "0.15", "0.20", "0.30", "0.40", "0.50", "0.70",
+                "0.80", "1.00", "1.20", "1.50", "2.00", "2.50", "3.00",
+                "4.00", "5.00", "7.00", "8.00", "10.0", "12.0", "15.0",
+                "20.0", "25.0"]
+
+    # Here we obtain the wavelength beforehand to avoid reading it each time.
+    datafile = open(dl2007_dir + "U{}/U{}_{}_MW3.1_{}.txt".format(umaximum[0],
+                                                                  umaximum[0],
+                                                                  umaximum[0],
+                                                                  "00"))
+    data = "".join(datafile.readlines()[-1001:])
+    datafile.close()
+
+    wave = np.genfromtxt(io.BytesIO(data.encode()), usecols=(0))
+    # For some reason wavelengths are decreasing in the model files
+    wave = wave[::-1]
+    # We convert wavelengths from μm to nm
+    wave *= 1000.
+
+    # The models are in Jy cm² sr¯¹ H¯¹. We convert this to W/nm.
+    conv = 4. * np.pi * 1e-30 / cst.m_p * cst.c / (wave * wave) * 1e9
+
+    for model in sorted(qpah.keys()):
+        for umin in uminimum:
+            filename = dl2007_dir + "U{}/U{}_{}_MW3.1_{}.txt".format(umin,
+                                                                     umin,
+                                                                     umin,
+                                                                     model)
+            print("Importing {}...".format(filename))
+            datafile = open(filename)
+            data = "".join(datafile.readlines()[-1001:])
+            datafile.close()
+            lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
+            # For some reason fluxes are decreasing in the model files
+            lumin = lumin[::-1]
+            # Conversion from Jy cm² sr¯¹ H¯¹ to W/nm
+            lumin *= conv
+
+            base.add_dl2007(DL2007(qpah[model], umin, umin, wave, lumin))
+            for umax in umaximum:
+                filename = dl2007_dir + "U{}/U{}_{}_MW3.1_{}.txt".format(umin,
+                                                                         umin,
+                                                                         umax,
+                                                                         model)
+                print("Importing {}...".format(filename))
+                datafile = open(filename)
+                data = "".join(datafile.readlines()[-1001:])
+                datafile.close()
+                lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
+                # For some reason fluxes are decreasing in the model files
+                lumin = lumin[::-1]
+
+                # Conversion from Jy cm² sr¯¹ H¯¹ to W/nm
+                lumin *= conv
+
+                base.add_dl2007(DL2007(qpah[model], umin, umax, wave, lumin))
+
+
 def build_fritz2006(base):
     fritz2006_dir = os.path.join(os.path.dirname(__file__), 'fritz2006/')
 
@@ -367,7 +442,12 @@ def build_base():
     print("\nDONE\n")
     print('#' * 78)
 
-    print("5- Importing Fritz et al. (2006) models\n")
+    print("5- Importing Draine and Li (2007) templates\n")
+    build_dl2007(base)
+    print("\nDONE\n")
+    print('#' * 78)
+
+    print("6- Importing Fritz et al. (2006) models\n")
     build_fritz2006(base)
     print("\nDONE\n")
     print('#' * 78)

pcigale/data/__init__.py

@@ -25,6 +25,7 @@ from .filters import Filter
 from .ssp_m2005 import SspM2005
 from .ssp_bc03 import SspBC03
 from .ir_templates_dh2002 import IrTemplatesDH2002
+from .ir_models_dl2007 import DL2007
 from .agn_fritz2006 import AgnFritz2006
 
 
@@ -119,6 +120,25 @@ class _DH2002InfraredTemplates(BASE):
         self.data = data
 
 
+class _DL2007(BASE):
+    """Storage for Draine and Li (2007) IR models
+    """
+
+    __tablename__ = 'DL2007_models'
+    qpah = Column(Float, primary_key=True)
+    umin = Column(Float, primary_key=True)
+    umax = Column(Float, primary_key=True)
+    wave = Column(PickleType)
+    lumin = Column(PickleType)
+
+    def __init__(self, model):
+        self.qpah = model.qpah
+        self.umin = model.umin
+        self.umax = model.umax
+        self.wave = model.wave
+        self.lumin = model.lumin
+
+
 class _Fritz2006AGN(BASE):
     """Storage for Fritz et al. (2006) models
     """
@@ -270,6 +290,25 @@ class Database(object):
         else:
             raise StandardError('The database is not writable.')
 
+    def add_dl2007(self, model):
+        """
+        Add a Draine and Li (2007) model to the databse.
+
+        Parameters
+        ----------
+        model: pcigale.data.DL2007
+
+        """
+        if self.is_writable:
+            self.session.add(_DL2007(model))
+            try:
+                self.session.commit()
+            except exc.IntegrityError:
+                self.session.rollback()
+                raise StandardError('The DL07 model is already in the base.')
+        else:
+            raise StandardError('The database is not writable.')
+
     def add_fritz2006_agn(self, agn):
         """
         Add a Fritz et al. (2006) AGN model to the database.
@@ -419,6 +458,34 @@ class Database(object):
         else:
             return None
 
+    def get_dl2007(self, qpah, umin, umax):
+        """
+        Get the Draine and Li (2007) model corresponding to the given set of
+        paramters.
+
+        Parameters
+        ----------
+        qpah: float
+            Mass fraction of PAH
+        umin: float
+            Minimum radiation field
+        umax: float
+            Maximum radiation field
+        gamma: float
+            Fraction of the dust exposed from Umin to Umax
+
+        """
+        result = (self.session.query(_DL2007).
+                  filter(_DL2007.qpah == qpah).
+                  filter(_DL2007.umin == umin).
+                  filter(_DL2007.umax == umax).
+                  first())
+        if result:
+            return DL2007(result.qpah, result.umin, result.umax, result.wave,
+                          result.lumin)
+        else:
+            return None
+
     def get_filter_list(self):
         """Get the list of the filters in the database.
 

pcigale/data/ir_models_dl2007.py

+# -*- 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: Médéric Boquien <mederic.boquien@oamp.fr>
+
+"""
+
+
+class DL2007(object):
+    """Draine and Li (2007) dust models
+
+    This class holds the data associated with the Draine and Li (2007)
+    dust models.
+
+    """
+
+    def __init__(self, qpah, umin, umax, wave, lumin):
+        """Create a new IR model
+
+        Parameters
+        ----------
+        qpah: float
+            Mass fraction of PAH
+        umin: float
+            Minimum radiation field illuminating the dust
+        umax: float
+            Maximum radiation field illuminating the dust
+        wave: array
+            Vector of the λ grid used in the templates [nm]
+        lumin: array
+            Model data in an array containing the luminosity density
+            versus the wavelength λ
+
+        """
+
+        self.qpah = qpah
+        self.umin = umin
+        self.umax = umax
+        self.wave = wave
+        self.lumin = lumin

pcigale/sed/modules/dl2007.py

+# -*- 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: Médéric Boquien <mederic.boquien@oamp.fr>
+
+"""
+
+
+from . import common
+import numpy as np
+from pcigale.data import Database
+
+
+class Module(common.SEDCreationModule):
+    """
+    Module computing the infra-red re-emission corresponding to an amount of
+    attenuation using the Draine and Li (2007) models.
+
+    Given an amount of attenuation (e.g. resulting from the action of a dust
+    attenuation module) this module normalises the Draine and Li (2007)
+    template corresponding to a given α to this amount of energy and add it
+    to the SED.
+
+    Information added to the SED: NAME_alpha.
+
+    """
+
+    parameter_list = {
+        'qpah': (
+            'float',
+            "Mass fraction of PAH",
+            None
+        ),
+        'umin': (
+            'float',
+            "Minimum radiation field",
+            None
+        ),
+        'umax': (
+            'float',
+            "Maximum radiation field",
+            None
+        ),
+        'gamma': (
+            'float',
+            "Fraction illuminated from Umin to Umax",
+            None
+        ),
+        'attenuation_value_names': (
+            'list of strings',
+            "List of attenuation value names (in the SED's info dictionary). "
+            "A new re-emission contribution will be added for each one.",
+            None
+        )
+    }
+
+    out_parameter_list = {
+                         'qpah': 'Mass fraction of PAH',
+                         'umin': 'Minimum radiation field',
+                         'umax': 'Maximum radiation field',
+                         'gamma': 'Fraction illuminated from Umin to Umax'
+                         }
+
+
+    def _init_code(self):
+        """Get the model out of the database"""
+
+        qpah = self.parameters["qpah"]
+        umin = self.parameters["umin"]
+        umax = self.parameters["umax"]
+        gamma = self.parameters["gamma"]
+
+        database = Database()
+        self.model_minmin = database.get_dl2007(qpah, umin, umin)
+        self.model_minmax = database.get_dl2007(qpah, umin, umax)
+        database.session.close_all()
+
+        # The models in memory are in W/nm for 1 kg of dust. At the same time
+        # we need to normalize them to 1 W here to easily scale them from the
+        # power absorbed in the UV-optical. If we want to retrieve the dust
+        # mass at a later point, we have to save their "emissivity" per unit
+        # mass in W kg¯¹, The gamma parameter does not affect the fact that it
+        # is for 1 kg because it represents a mass fraction of each component.
+        self.emissivity = np.trapz((1. - gamma) * self.model_minmin.lumin +
+                          gamma * self.model_minmax.lumin,
+                          x=self.model_minmin.wave)
+
+        # We want to be able to display the respective constributions of both
+        # components, therefore we keep they separately.
+        self.model_minmin.lumin *= (1. - gamma) / self.emissivity
+        self.model_minmax.lumin *= gamma / self.emissivity
+
+
+    def process(self, sed):
+        """Add the IR re-emission contributions
+
+        Parameters
+        ----------
+        sed  : pcigale.sed.SED object
+        parameters : dictionary containing the parameters
+
+        """
+
+        # Base name for adding information to the SED.
+        name = self.name or 'dl2007'
+
+        sed.add_module(name, self.parameters)
+        sed.add_info(name + '_qpah', self.parameters["qpah"])
+        sed.add_info(name + '_umin', self.parameters["umin"])
+        sed.add_info(name + '_umax', self.parameters["umax"])
+        sed.add_info(name + '_gamma', self.parameters["gamma"])
+
+        for attenuation in self.parameters['attenuation_value_names']:
+            sed.add_contribution(
+                name + '_Umin_Umin_' + attenuation,
+                self.model_minmin.wave,
+                sed.info[attenuation] * self.model_minmin.lumin
+            )
+            sed.add_contribution(
+                name + '_Umin_Umax_' + attenuation,
+                self.model_minmax.wave,
+                sed.info[attenuation] * self.model_minmax.lumin
+            )