Integrate the polar dust model of X-CIGALE (Yang et al., 2020) for the skirtor2016 module.

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 ### Added
+- The polar dust model of X-CIGALE (Yang et al., 2020) for the `skirtor2016` module has been integrated into the regular version. (Médéric Boquien, based on the initial work of Guang Yang)
 ### Changed
 ### Fixed
 - Ensure that `pcigale-plots` correctly detects the `skirtor2016` AGN models. (Médéric Boquien, reported by Guang Yang)

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
+# Authors: Médéric Boquien, Laure Ciesla, Guang Yang
 
 """
 SKIRTOR 2016 (Stalevski et al., 2016) AGN dust torus emission module
@@ -13,10 +12,60 @@ This module implements the SKIRTOR 2016 models.
 from collections import OrderedDict
 
 import numpy as np
+import scipy.constants as cst
 
 from pcigale.data import Database
 from . import SedModule
 
+def k_ext(wavelength, ext_law):
+    """
+    Compute k(λ)=A(λ)/E(B-V) for a specified extinction law
+
+    Parameters
+    ----------
+    wavelength: array of floats
+        Wavelength grid in nm.
+    ext_law: the extinction law
+             0=SMC, 1=Calzetti2000, 2=Gaskell2004
+
+    Returns
+    -------
+    a numpy array of floats
+
+    """
+    if ext_law == 0:
+        # SMC, from Bongiorno+2012
+        return 1.39 * (wavelength * 1e-3) ** -1.2
+    elif ext_law == 1:
+        # Calzetti2000, from dustatt_calzleit.py
+        result = np.zeros(len(wavelength))
+        # Attenuation between 120 nm and 630 nm
+        mask = (wavelength < 630)
+        result[mask] = 2.659 * (-2.156 + 1.509e3 / wavelength[mask] -
+                                0.198e6 / wavelength[mask] ** 2 +
+                                0.011e9 / wavelength[mask] ** 3) + 4.05
+        # Attenuation between 630 nm and 2200 nm
+        mask = (wavelength >= 630)
+        result[mask] = 2.659 * (-1.857 + 1.040e3 / wavelength[mask]) + 4.05
+        return result
+    elif ext_law == 2:
+        # Gaskell+2004, from the appendix of that paper
+        x = 1e3 / wavelength
+        Alam_Av = np.zeros(len(wavelength))
+        # Attenuation for x = 1.6 -- 3.69
+        mask = (x < 3.69)
+        Alam_Av[mask] = -0.8175 + 1.5848*x[mask] - 0.3774*x[mask]**2 + 0.0296*x[mask]**3
+        # Attenuation for x = 3.69 -- 8
+        mask = (x >= 3.69)
+        Alam_Av[mask] = 1.3468 + 0.0087*x[mask]
+        # Set negative values to zero
+        Alam_Av[Alam_Av < 0.] = 0.
+        # Convert A(λ)/A(V) to A(λ)/E(B-V)
+        # assuming A(B)/A(V) = 1.182 (Table 3 of Gaskell+2004)
+        return Alam_Av / 0.182
+    else:
+        raise KeyError("Extinction law is different from the expected ones")
+
 
 class SKIRTOR2016(SedModule):
     """SKIRTOR 2016 (Stalevski et al., 2016) AGN dust torus emission
@@ -82,6 +131,27 @@ class SKIRTOR2016(SedModule):
             'cigale_list(minvalue=0., maxvalue=1.)',
             "AGN fraction.",
             0.1
+        )),
+        ('law', (
+            'cigale_list(dtype=int, options=0 & 1 & 2)',
+            "Extinction law of the polar dust: "
+            "0 (SMC), 1 (Calzetti 2000), or 2 (Gaskell et al. 2004)",
+            0
+        )),
+        ('EBV', (
+            'cigale_list(minvalue=0.)',
+            "E(B-V) for the extinction in the polar direction in magnitudes.",
+            0.1
+        )),
+        ('temperature', (
+            'cigale_list(minvalue=0.)',
+            "Temperature of the polar dust in K.",
+            100.
+        )),
+        ("emissivity", (
+            "cigale_list(minvalue=0.)",
+            "Emissivity index of the polar dust.",
+            1.6
         ))
     ])
 
@@ -95,11 +165,69 @@ class SKIRTOR2016(SedModule):
         self.Mcl = float(self.parameters["Mcl"])
         self.i = int(self.parameters["i"])
         self.fracAGN = float(self.parameters["fracAGN"])
+        self.law = int(self.parameters["law"])
+        self.EBV = float(self.parameters["EBV"])
+        self.temperature = float(self.parameters["temperature"])
+        self.emissivity = float(self.parameters["emissivity"])
 
         with Database() as base:
             self.SKIRTOR2016 = base.get_skirtor2016(self.t, self.pl, self.q,
                                                     self.oa, self.R, self.Mcl,
                                                     self.i)
+            self.AGN1 = base.get_skirtor2016(self.t, self.pl, self.q, self.oa,
+                                             self.R, self.Mcl, 0.)
+
+        # Calculate the extinction
+        ext_fac = 10 ** (-.4*k_ext(self.SKIRTOR2016.wave, self.law) * self.EBV)
+
+        # Calculate the new AGN SED shape after extinction
+        # The direct and scattered components (line-of-sight) are extincted for
+        # type-1 AGN
+        # Keep the direct and scatter components for type-2
+        if self.i <= (90. - self.oa):
+            self.SKIRTOR2016.disk *= ext_fac
+
+        # Calculate the total extincted luminosity averaged over all directions
+        # The computation is non-trivial as the disk emission is anisotropic:
+        # L(θ, λ) = A(λ)×cosθ×(1+2×cosθ). Given that L(θ=0, λ) = A(λ)×3, we get
+        # A = L(θ=0, λ)/3. Then Lpolar = 2×∭L(θ, λ)×(1-ext_fact(λ))×sinθ dφdθdλ,
+        # woth φ between 0 and 2π, θ between 0 and π/2-OA, and λ the wavelengths
+        # Integrating over φ,
+        # Lpolar = 4π/3×∬L(θ=0, λ)×(1-ext_fact(λ))×cosθ×(1+2×cosθ)×sinθ dθdλ.
+        # Now doing the usual integration over θ,
+        # Lpolar = 4π/3×∫L(θ=0, λ)×(1-ext_fact(λ))×[7/6-1/2×sin²OA-2/3×sin³OA] dλ.
+        # Now a critical point is that the SKIRTOR models are provided in flux
+        # and are multiplied by 4πd². Because the emission is anisotropic, we
+        # need to redivide by 4π to get the correct luminosity for a given θ,
+        # hence Lpolar = [7/18-1/6×sin²OA-2/9×sin³OA]×∫L(θ=0, λ)×(1-ext_fact(λ)) dλ.
+        # Integrating over λ gives the bolometric luminosity
+        sin_oa = np.sin(np.deg2rad(self.oa))
+        l_ext = (7./18. - sin_oa**2/6. - sin_oa**3*2./9.) * \
+                np.trapz(self.AGN1.disk * (1. - ext_fac),
+                         x=self.SKIRTOR2016.wave)
+
+        # Casey (2012) modified black body model
+        c = cst.c * 1e9
+        lambda_0 = 200e3
+        conv = c / self.SKIRTOR2016.wave ** 2.
+        hc_lkt = cst.h * c / (self.SKIRTOR2016.wave * cst.k * self.temperature)
+        err_settings = np.seterr(over='ignore')  # ignore exp overflow
+        blackbody = conv * \
+            (1. - np.exp(-(lambda_0 / self.SKIRTOR2016.wave) ** self.emissivity)) * \
+            (c / self.SKIRTOR2016.wave) ** 3. / (np.exp(hc_lkt) - 1.)
+        np.seterr(**err_settings)  # Restore the previous settings
+        blackbody *= l_ext / np.trapz(blackbody, x=self.SKIRTOR2016.wave)
+
+        # Add the black body to dust thermal emission
+        self.SKIRTOR2016.dust += blackbody
+
+        # Normalize direct, scatter, and thermal components
+        norm = 1. / np.trapz(self.SKIRTOR2016.dust, x=self.SKIRTOR2016.wave)
+        self.SKIRTOR2016.dust *= norm
+        self.SKIRTOR2016.disk *= norm
+
+        # Integrate AGN luminosity for different components
+        self.lumin_disk = np.trapz(self.SKIRTOR2016.disk, x=self.SKIRTOR2016.wave)
 
     def process(self, sed):
         """Add the IR re-emission contributions
@@ -124,6 +252,10 @@ class SKIRTOR2016(SedModule):
         sed.add_info('agn.Mcl', self.Mcl)
         sed.add_info('agn.i', self.i, unit='deg')
         sed.add_info('agn.fracAGN', self.fracAGN)
+        sed.add_info('agn.law', self.law)
+        sed.add_info('agn.EBV', self.EBV)
+        sed.add_info('agn.temperature', self.temperature)
+        sed.add_info('agn.emissivity', self.emissivity)
 
         # Compute the AGN luminosity
         if self.fracAGN < 1.: