Add the themis data structures and database routines, and build the models into the database.

database_builder/__init__.py

@@ -22,7 +22,7 @@ import scipy.constants as cst
 from astropy.table import Table
 from pcigale.data import (Database, Filter, M2005, BC03, Fritz2006,
                           Dale2014, DL2007, DL2014, NebularLines,
-                          NebularContinuum, Schreiber2016)
+                          NebularContinuum, Schreiber2016, THEMIS)
 
 
 def read_bc03_ssp(filename):
@@ -769,6 +769,76 @@ def build_schreiber2016(base):
     base.add_schreiber2016(models)
 
 
+def build_themis(base):
+    models = []
+    themis_dir = os.path.join(os.path.dirname(__file__), 'themis/')
+
+    # Mass fraction of hydrocarbon solids i.e., a-C(:H) smaller than 1.5 nm,
+    # also known as HAC
+    qhac = {"000": 0.02, "010": 0.06, "020": 0.10, "030": 0.14, "040": 0.17,
+            "050": 0.20, "060": 0.24, "070": 0.28, "080": 0.32, "090": 0.36,
+            "100": 0.40}
+
+    uminimum = ["0.100", "0.120", "0.150", "0.170", "0.200", "0.250", "0.300",
+                "0.350", "0.400", "0.500", "0.600", "0.700", "0.800", "1.000",
+                "1.200", "1.500", "1.700", "2.000", "2.500", "3.000", "3.500",
+                "4.000", "5.000", "6.000", "7.000", "8.000", "10.00", "12.00",
+                "15.00", "17.00", "20.00", "25.00", "30.00", "35.00", "40.00",
+                "50.00", "80.00"]
+
+    alpha = ["1.0", "1.1", "1.2", "1.3", "1.4", "1.5", "1.6", "1.7", "1.8",
+             "1.9", "2.0", "2.1", "2.2", "2.3", "2.4", "2.5", "2.6", "2.7",
+             "2.8", "2.9", "3.0"]
+
+    # Mdust/MH used to retrieve the dust mass as models as given per atom of H
+    MdMH = {"000": 7.4e-3, "010": 7.4e-3, "020": 7.4e-3, "030": 7.4e-3,
+            "040": 7.4e-3, "050": 7.4e-3, "060": 7.4e-3, "070": 7.4e-3,
+            "080": 7.4e-3, "090": 7.4e-3, "100": 7.4e-3}
+
+    # Here we obtain the wavelength beforehand to avoid reading it each time.
+    datafile = open(themis_dir + "U{}_{}_MW3.1_{}/spec_1.0.dat"
+                    .format(uminimum[0], uminimum[0], "000"))
+
+    data = "".join(datafile.readlines()[-576:])
+    datafile.close()
+
+    wave = np.genfromtxt(io.BytesIO(data.encode()), usecols=(0))
+
+    # We convert wavelengths from μm to nm
+    wave *= 1000.
+
+    # Conversion factor from Jy cm² sr¯¹ H¯¹ to W nm¯¹ (kg of H)¯¹
+    conv = 4. * np.pi * 1e-30 / (cst.m_p+cst.m_e) * cst.c / (wave*wave) * 1e9
+
+    for model in sorted(qhac.keys()):
+        for umin in uminimum:
+            filename = (themis_dir + "U{}_{}_MW3.1_{}/spec_1.0.dat"
+                        .format(umin, umin, model))
+            print("Importing {}...".format(filename))
+            with open(filename) as datafile:
+                data = "".join(datafile.readlines()[-576:])
+            lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
+
+            # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
+            lumin *= conv / MdMH[model]
+
+            models.append(THEMIS(qhac[model], umin, umin, 1.0, wave, lumin))
+            for al in alpha:
+                filename = (themis_dir + "U{}_1e7_MW3.1_{}/spec_{}.dat"
+                            .format(umin, model, al))
+                print("Importing {}...".format(filename))
+                with open(filename) as datafile:
+                    data = "".join(datafile.readlines()[-576:])
+                lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
+
+                # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
+                lumin *= conv/MdMH[model]
+
+                models.append(THEMIS(qhac[model], umin, 1e7, al, wave, lumin))
+
+    base.add_themis(models)
+
+
 def build_base(bc03res='lr'):
     base = Database(writable=True)
     base.upgrade_base()
@@ -819,7 +889,12 @@ def build_base(bc03res='lr'):
     build_schreiber2016(base)
     print("\nDONE\n")
     print('#' * 78)
-    
+
+    print("10- Importing Jones et al (2017) models)\n")
+    build_themis(base)
+    print("\nDONE\n")
+    print('#' * 78)
+
     base.session.close_all()
 
 

pcigale/data/__init__.py

@@ -31,6 +31,7 @@ from .fritz2006 import Fritz2006
 from .nebular_continuum import NebularContinuum
 from .nebular_lines import NebularLines
 from .schreiber2016 import Schreiber2016
+from .themis import THEMIS
 
 DATABASE_FILE = pkg_resources.resource_filename(__name__, 'data.db')
 
@@ -252,6 +253,27 @@ class _Schreiber2016(BASE):
         self.lumin = ir.lumin
 
 
+class _THEMIS(BASE):
+    """Storage for the Jones et al (2017) IR models
+    """
+
+    __tablename__ = 'THEMIS_models'
+    qhac = Column(Float, primary_key=True)
+    umin = Column(Float, primary_key=True)
+    umax = Column(Float, primary_key=True)
+    alpha = Column(Float, primary_key=True)
+    wave = Column(PickleType)
+    lumin = Column(PickleType)
+
+    def __init__(self, model):
+        self.qhac = model.qhac
+        self.umin = model.umin
+        self.umax = model.umax
+        self.alpha = model.alpha
+        self.wave = model.wave
+        self.lumin = model.lumin
+
+
 class Database(object):
     """Object giving access to pcigale database."""
 
@@ -873,6 +895,69 @@ class Database(object):
         """
         return self._get_parameters(_Schreiber2016)
 
+    def add_themis(self, models):
+        """
+        Add a list of Jones et al (2017) models to the database.
+
+        Parameters
+        ----------
+        models: list of pcigale.data.THEMIS objects
+
+        """
+        if self.is_writable:
+            for model in models:
+                self.session.add(_THEMIS(model))
+            try:
+                self.session.commit()
+            except exc.IntegrityError:
+                self.session.rollback()
+                raise DatabaseInsertError(
+                    'Error.')
+        else:
+            raise Exception('The database is not writable.')
+
+    def get_themis(self, qhac, umin, umax, alpha):
+        """
+        Get the Jones et al (2017) model corresponding to the given set of
+        parameters.
+
+        Parameters
+        ----------
+        qhac: float
+            Mass fraction of hydrocarbon solids i.e., a-C(:H) smaller than
+        1.5 nm, also known as HAC
+        umin: float
+            Minimum radiation field
+        umin: float
+            Maximum radiation field
+        alpha: float
+            Powerlaw slope dU/dM∝U¯ᵅ
+
+        Returns
+        -------
+        model: pcigale.data.THEMIS
+            The Jones et al (2017) model.
+
+        Raises
+        ------
+        DatabaseLookupError: if the requested model is not in the database.
+
+        """
+        result = (self.session.query(_THEMIS).
+                  filter(_THEMIS.qhac == qhac).
+                  filter(_THEMIS.umin == umin).
+                  filter(_THEMIS.umax == umax).
+                  filter(_THEMIS.alpha == alpha).
+                  first())
+        if result:
+            return THEMIS(result.qhac, result.umin, result.umax, result.alpha,
+                          result.wave, result.lumin)
+        else:
+            raise DatabaseLookupError(
+                "The THEMIS model for qhac <{0}>, umin <{1}>, umax <{2}>, and "
+                "alpha <{3}> is not in the database."
+                .format(qhac, umin, umax, alpha))
+
     def _get_parameters(self, schema):
         """Generic function to get parameters from an arbitrary schema.
 
@@ -904,6 +989,16 @@ class Database(object):
         else:
             raise Exception('The database is not writable.')
 
+    def get_themis_parameters(self):
+        """Get parameters for the THEMIS models.
+
+        Returns
+        -------
+        paramaters: dictionary
+            dictionary of parameters and their values
+        """
+        return self._get_parameters(_THEMIS)
+
     def add_filters(self, pcigale_filters):
         """
         Add a list of filters to the pcigale database.

pcigale/data/themis.py

+# -*- coding: utf-8 -*-
+# DustPedia,  http://www.dustpedia.com/
+# Author: Angelos Nersesian, Frederic Galliano, Anthony Jones, Pieter De Vis
+
+
+class THEMIS(object):
+    """Jones et al (2017) dust models,
+
+    This class holds the data associated with the Jones et al (2017)
+    dust models.
+
+    """
+
+    def __init__(self, qhac, umin, umax, alpha, wave, lumin):
+        """Create a new IR model
+
+        Parameters
+        ----------
+        qhac: float
+            Mass fraction of hydrocarbon solids i.e., a-C(:H) smaller than
+        1.5 nm, also known as HAC
+        umin: float
+            Minimum radiation field illuminating the dust
+        umax: float
+            Maximum radiation field illuminating the dust
+        alpha: float
+            Powerlaw slope dU/dM∝U¯ᵅ
+        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.qhac = qhac
+        self.umin = umin
+        self.umax = umax
+        self.alpha = alpha
+        self.wave = wave
+        self.lumin = lumin