Initial implementation of the SKIRTOR 2016 models.

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, THEMIS)
+                          NebularContinuum, SKIRTOR2016, Schreiber2016, THEMIS)
 
 
 def read_bc03_ssp(filename):
@@ -694,6 +694,55 @@ def build_fritz2006(base):
 
     base.add_fritz2006(models)
 
+def build_skirtor2016(base):
+    models = []
+    skirtor2016_dir = os.path.join(os.path.dirname(__file__), 'skirtor2016/')
+
+    files = glob.glob(skirtor2016_dir + '/*')
+    files = [file.split('/')[-1] for file in files]
+    params = [f.split('_')[:-1] for f in files]
+
+    # Parameters of SKIRTOR 2016
+    t = list({param[0][1:] for param in params})
+    p = list({param[1][1:] for param in params})
+    q = list({param[2][1:] for param in params})
+    oa = list({param[3][2:] for param in params})
+    R = list({param[4][1:] for param in params})
+    Mcl = list({param[5][3:] for param in params})
+    i = list({param[6][1:] for param in params})
+
+    iter_params = ((p1, p2, p3, p4, p5, p6, p7)
+                   for p1 in t
+                   for p2 in p
+                   for p3 in q
+                   for p4 in oa
+                   for p5 in R
+                   for p6 in Mcl
+                   for p7 in i)
+
+    for params in iter_params:
+        filename = skirtor2016_dir + \
+                "t{}_p{}_q{}_oa{}_R{}_Mcl{}_i{}_sed.dat".format(*params)
+        print("Importing {}...".format(filename))
+
+        wl, disk, scatt, dust = np.genfromtxt(filename, unpack=True,
+                                              usecols=(0, 2, 3, 4))
+        wl *= 1e3
+        disk += scatt
+        disk /= wl
+        dust /= wl
+
+        # Normalization of the lumin_therm to 1W
+        norm = np.trapz(dust, x=wl)
+        disk /= norm
+        dust /= norm
+
+        models.append(SKIRTOR2016(params[0], params[1], params[2], params[3],
+                                  params[4], params[5], params[6], wl, disk,
+                                  dust))
+
+    base.add_skirtor2016(models)
+
 def build_nebular(base):
     models_lines = []
     models_cont = []
@@ -884,22 +933,27 @@ def build_base(bc03res='lr'):
     print("\nDONE\n")
     print('#' * 78)
 
-    print("7- Importing Dale et al (2014) templates\n")
+    print("7- Importing SKIRTOR 2016 models\n")
+    build_skirtor2016(base)
+    print("\nDONE\n")
+    print('#' * 78)
+
+    print("8- Importing Dale et al (2014) templates\n")
     build_dale2014(base)
     print("\nDONE\n")
     print('#' * 78)
 
-    print("8- Importing nebular lines and continuum\n")
+    print("9- Importing nebular lines and continuum\n")
     build_nebular(base)
     print("\nDONE\n")
     print('#' * 78)
 
-    print("9- Importing Schreiber et al (2016) models\n")
+    print("10- Importing Schreiber et al (2016) models\n")
     build_schreiber2016(base)
     print("\nDONE\n")
     print('#' * 78)
 
-    print("10- Importing Jones et al (2017) models)\n")
+    print("11- Importing Jones et al (2017) models)\n")
     build_themis(base)
     print("\nDONE\n")
     print('#' * 78)

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 .skirtor2016 import SKIRTOR2016
 from .themis import THEMIS
 
 DATABASE_FILE = pkg_resources.resource_filename(__name__, 'data.db')
@@ -202,6 +203,35 @@ class _Fritz2006(BASE):
         self.lumin_agn = agn.lumin_agn
 
 
+class _SKIRTOR2016(BASE):
+    """Storage for SKIRTOR 2016 models
+    """
+
+    __tablename__ = 'skirtor2016'
+    t = Column(Float, primary_key=True)
+    pl = Column(Float, primary_key=True)
+    q = Column(Float, primary_key=True)
+    oa = Column(Float, primary_key=True)
+    R = Column(Float, primary_key=True)
+    Mcl = Column(Float, primary_key=True)
+    i = Column(Float, primary_key=True)
+    wave = Column(PickleType)
+    disk = Column(PickleType)
+    dust = Column(PickleType)
+
+    def __init__(self, agn):
+        self.t = agn.t
+        self.pl = agn.pl
+        self.q = agn.q
+        self.oa = agn.oa
+        self.R = agn.R
+        self.Mcl = agn.Mcl
+        self.i = agn.i
+        self.wave = agn.wave
+        self.disk = agn.disk
+        self.dust = agn.dust
+
+
 class _NebularLines(BASE):
     """Storage for line templates
     """
@@ -737,6 +767,96 @@ class Database(object):
         """
         return self._get_parameters(_Fritz2006)
 
+    def add_skirtor2016(self, models):
+        """
+        Add a SKIRTOR 2016 (Stalevski et al., 2016) AGN model to the database.
+
+        Parameters
+        ----------
+        models: list of pcigale.data.SKIRTOR2016 objects
+
+        """
+        if self.is_writable:
+            for model in models:
+                self.session.add(_SKIRTOR2016(model))
+            try:
+                self.session.commit()
+            except exc.IntegrityError:
+                self.session.rollback()
+                raise DatabaseInsertError(
+                    'The agn model is already in the base.')
+        else:
+            raise Exception('The database is not writable.')
+
+    def get_skirtor2016(self, t, pl, q, oa, R, Mcl, i):
+        """
+        Get the SKIRTOR 2016 AGN model corresponding to a given set of
+        parameters.
+
+        Parameters
+        ----------
+        t: float
+            average edge-on optical depth at 9.7 micron; the actual one along
+            the line of sight may vary depending on the clumps distribution
+        pl: float
+            power-law exponent that sets radial gradient of dust density
+        q: float
+            index that sets dust density gradient with polar angle
+        oa: float
+            angle measured between the equatorial plan and edge of the torus.
+            Half-opening angle of the dust-free cone is 90-oa
+        R: float
+            ratio of outer to inner radius, R_out/R_in
+        Mcl: float
+            Angle between AGN axis and line of sight.
+        i: float
+            inclination, i.e. viewing angle, i.e. position of the instrument
+            w.r.t. the AGN axis. i=0: face-on, type 1 view; i=90: edge-on, type
+            2 view.
+        wave: array of float
+            Wavelength grid in nm.
+        disk: array of flaot
+            Luminosity of the accretion disk in W/nm
+        dust: array of float
+            Luminosity of the dust in W/nm
+
+        Returns
+        -------
+        agn: pcigale.data.SKIRTOR2016
+            The AGN model.
+
+        Raises
+        ------
+        DatabaseLookupError: if the requested template is not in the database.
+
+        """
+        result = (self.session.query(_SKIRTOR2016).
+                  filter(_SKIRTOR2016.t == t).
+                  filter(_SKIRTOR2016.pl == pl).
+                  filter(_SKIRTOR2016.q == q).
+                  filter(_SKIRTOR2016.oa == oa).
+                  filter(_SKIRTOR2016.R == R).
+                  filter(_SKIRTOR2016.Mcl == Mcl).
+                  filter(_SKIRTOR2016.i == i).
+                  first())
+        if result:
+            return SKIRTOR2016(result.t, result.pl, result.q, result.oa,
+                               result.R, result.Mcl, result.i, result.wave,
+                               result.disk, result.dust)
+        else:
+            raise DatabaseLookupError(
+                "The SKIRTOR2016 model is not in the database.")
+
+    def get_skirtor2016_parameters(self):
+        """Get parameters for the SKIRTOR 2016 AGN models.
+
+        Returns
+        -------
+        paramaters: dictionary
+            dictionary of parameters and their values
+        """
+        return self._get_parameters(_SKIRTOR2016)
+
     def add_nebular_lines(self, models):
         """
         Add ultraviolet and optical line templates to the database.

pcigale/data/skirtor2016.py

+class SKIRTOR2016(object):
+    """SKIRTOR (Stalevski et al., 2016) AGN dust torus emission model.
+
+    This class holds the UV-optical data associated with a SKIRTORN AGN model
+    (Stalevski et al., 2012, 2016).
+
+    """
+
+    def __init__(self, t, pl, q, oa, R, Mcl, i, wave, disk, dust):
+        """Create a new AGN model. The descriptions of the parameters are taken
+           directly from https://sites.google.com/site/skirtorus/sed-library.
+
+        Parameters
+        ----------
+        t: float
+            average edge-on optical depth at 9.7 micron; the actual one along
+            the line of sight may vary depending on the clumps distribution
+        pl: float
+            power-law exponent that sets radial gradient of dust density
+        q: float
+            index that sets dust density gradient with polar angle
+        oa: float
+            angle measured between the equatorial plan and edge of the torus.
+            Half-opening angle of the dust-free cone is 90-oa
+        R: float
+            ratio of outer to inner radius, R_out/R_in
+        Mcl: float
+            fraction of total dust mass inside clumps. 0.97 means 97% of total
+            mass is inside the clumps and 3% in the interclump dust.
+        i: float
+            inclination, i.e. viewing angle, i.e. position of the instrument
+            w.r.t. the AGN axis. i=0: face-on, type 1 view; i=90: edge-on, type
+            2 view.
+        wave: array of float
+            Wavelength grid in nm.
+        disk: array of flaot
+            Luminosity of the accretion disk in W/nm
+        dust: array of float
+            Luminosity of the dust in W/nm
+
+        """
+
+        self.t = t
+        self.pl = pl
+        self.q = q
+        self.oa = oa
+        self.R = R
+        self.Mcl = Mcl
+        self.i = i
+        self.wave = wave
+        self.disk = disk
+        self.dust = dust

pcigale/sed_modules/skirtor2016.py

+# -*- coding: utf-8 -*-
+# Copyright (C) 2013, 2014 Department of Physics, University of Crete
+# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
+# Author: Laure Ciesla
+
+"""
+SKIRTOR 2016 (Stalevski et al., 2016) AGN dust torus emission module
+==================================================
+
+This module implements the SKIRTOR 2016 models.
+
+"""
+from collections import OrderedDict
+
+import numpy as np
+
+from pcigale.data import Database
+from . import SedModule
+
+
+class SKIRTOR2016(SedModule):
+    """SKIRTOR 2016 (Stalevski et al., 2016) AGN dust torus emission
+
+
+    The relative normalization of these components is handled through a
+    parameter which is the fraction of the total IR luminosity due to the AGN
+    so that: L_AGN = fracAGN * L_IRTOT, where L_AGN is the AGN luminosity,
+    fracAGN is the contribution of the AGN to the total IR luminosity
+    (L_IRTOT), i.e. L_Starburst+L_AGN.
+
+    """
+
+    parameter_list = OrderedDict([
+        ('t', (
+            "cigale_list(options=3 & 5 & 7 & 9 & 11)",
+            "Average edge-on optical depth at 9.7 micron; the actual one along"
+            "the line of sight may vary depending on the clumps distribution. "
+            "Possible values are: 3, 5, 7, 8, and 11.",
+            3
+        )),
+        ('pl', (
+            "cigale_list(options=0. & .5 & 1. & 1.5)",
+            "Power-law exponent that sets radial gradient of dust density."
+            "Possible values are: 0., 0.5, 1., and 1.5.",
+            1.0
+        )),
+        ('q', (
+            "cigale_list(options=0. & .5 & 1. & 1.5)",
+            "Index that sets dust density gradient with polar angle."
+            "Possible values are:  0., 0.5, 1., and 1.5.",
+            1.0
+        )),
+        ('oa', (
+            'cigale_list(options=10 & 20 & 30 & 40 & 50 & 60 & 70 & 80)',
+            "Angle measured between the equatorial plan and edge of the torus. "
+            "Half-opening angle of the dust-free cone is 90-oa"
+            "Possible values are: 10, 20, 30, 40, 50, 60, 70, and 80",
+            40
+        )),
+        ('R', (
+            'cigale_list(options=10 & 20 & 30)',
+            "Ratio of outer to inner radius, R_out/R_in."
+            "Possible values are: 10, 20, and 30",
+            20
+        )),
+        ('Mcl', (
+            'cigale_list(options=0.97)',
+            "fraction of total dust mass inside clumps. 0.97 means 97% of "
+            "total mass is inside the clumps and 3% in the interclump dust. "
+            "Possible values are: 0.97.",
+            0.97
+        )),
+        ('i', (
+            'cigale_list(options=0 & 10 & 20 & 30 & 40 & 50 & 60 & 70 & 80 & 90)',
+            "inclination, i.e. viewing angle, i.e. position of the instrument "
+            "w.r.t. the AGN axis. i=0: face-on, type 1 view; i=90: edge-on, "
+            "type 2 view."
+            "Possible values are: 0, 10, 20, 30, 40, 50, 60, 70, 80, and 90.",
+            40
+        )),
+        ('fracAGN', (
+            'cigale_list(minvalue=0., maxvalue=1.)',
+            "AGN fraction.",
+            0.1
+        ))
+    ])
+
+    def _init_code(self):
+        """Get the template set out of the database"""
+        self.t = int(self.parameters["t"])
+        self.pl = float(self.parameters["pl"])
+        self.q = float(self.parameters["q"])
+        self.oa = int(self.parameters["oa"])
+        self.R = int(self.parameters["R"])
+        self.Mcl = float(self.parameters["Mcl"])
+        self.i = int(self.parameters["i"])
+        self.fracAGN = float(self.parameters["fracAGN"])
+
+        with Database() as base:
+            self.SKIRTOR2016 = base.get_skirtor2016(self.t, self.pl, self.q,
+                                                    self.oa, self.R, self.Mcl,
+                                                    self.i)
+
+    def process(self, sed):
+        """Add the IR re-emission contributions
+
+        Parameters
+        ----------
+        sed: pcigale.sed.SED object
+        parameters: dictionary containing the parameters
+
+        """
+
+        if 'dust.luminosity' not in sed.info:
+            sed.add_info('dust.luminosity', 1., True)
+        luminosity = sed.info['dust.luminosity']
+
+        sed.add_module(self.name, self.parameters)
+        sed.add_info('agn.t', self.t)
+        sed.add_info('agn.pl', self.pl)
+        sed.add_info('agn.q', self.q)
+        sed.add_info('agn.oa', self.oa)
+        sed.add_info('agn.R', self.R)
+        sed.add_info('agn.Mcl', self.Mcl)
+        sed.add_info('agn.i', self.i)
+        sed.add_info('agn.fracAGN', self.fracAGN)
+
+        # Compute the AGN luminosity
+        if self.fracAGN < 1.:
+            agn_power = luminosity * (1./(1.-self.fracAGN) - 1.)
+            lumin_dust = agn_power
+            lumin_disk = agn_power * np.trapz(self.SKIRTOR2016.disk,
+                                              x=self.SKIRTOR2016.wave)
+            lumin = lumin_dust + lumin_disk
+        else:
+            raise Exception("AGN fraction is exactly 1. Behaviour "
+                            "undefined.")
+
+        sed.add_info('agn.dust_luminosity', lumin_dust, True)
+        sed.add_info('agn.disk_luminosity', lumin_disk, True)
+        sed.add_info('agn.luminosity', lumin, True)
+
+        sed.add_contribution('agn.SKIRTOR2016_dust', self.SKIRTOR2016.wave,
+                             agn_power * self.SKIRTOR2016.dust)
+        sed.add_contribution('agn.SKIRTOR2016_disk', self.SKIRTOR2016.wave,
+                             agn_power * self.SKIRTOR2016.disk)
+
+
+# SedModule to be returned by get_module
+Module = SKIRTOR2016