modify skirtor from anisotorpic to isotropic; add polar dust for skirtor;...

modify skirtor from anisotorpic to isotropic; add polar dust for skirtor; connect xray with skirtor; bug fixed for fritz

database_builder/__init__.py

@@ -663,7 +663,7 @@ def build_fritz2006(base):
         data = datafile.readlines()
         datafile.close()
 
-        # Get intrinsic (de-reddened) AGN luminosity, i.e., 
+        # Get intrinsic (de-reddened) AGN luminosity, i.e.,
         # central AGN luminosity at psy=89.990
         n = len(psy)-1
         block = data[nskip + blocksize * n + 4 * (n + 1) - 1:
@@ -731,6 +731,12 @@ def build_skirtor2016(base):
                    for p6 in Mcl
                    for p7 in i)
 
+    # Extract the intrinsic AGN disk luminosity
+    filename = skirtor2016_dir+"t9_p1_q1_oa80_R30_Mcl0.97_i0_sed.dat"
+    intrin_disk_unnormed = np.genfromtxt(filename, unpack=True, usecols=(6))
+    # Convert anisotropic to the isotropic disk luminosity
+    intrin_iso_disk_unnormed = intrin_disk_unnormed*7/18
+
     for params in iter_params:
         filename = skirtor2016_dir + \
                 "t{}_p{}_q{}_oa{}_R{}_Mcl{}_i{}_sed.dat".format(*params)
@@ -747,10 +753,11 @@ def build_skirtor2016(base):
         norm = np.trapz(dust, x=wl)
         disk /= norm
         dust /= norm
+        intrin_iso_disk = intrin_iso_disk_unnormed / (norm*wl)
 
         models.append(SKIRTOR2016(params[0], params[1], params[2], params[3],
                                   params[4], params[5], params[6], wl, disk,
-                                  dust))
+                                  dust, intrin_iso_disk))
 
     base.add_skirtor2016(models)
 

pcigale/data/__init__.py

@@ -220,6 +220,7 @@ class _SKIRTOR2016(BASE):
     wave = Column(PickleType)
     disk = Column(PickleType)
     dust = Column(PickleType)
+    intrin_iso_disk = Column(PickleType)
 
     def __init__(self, agn):
         self.t = agn.t
@@ -232,7 +233,7 @@ class _SKIRTOR2016(BASE):
         self.wave = agn.wave
         self.disk = agn.disk
         self.dust = agn.dust
-
+        self.intrin_iso_disk = agn.intrin_iso_disk
 
 class _NebularLines(BASE):
     """Storage for line templates
@@ -731,7 +732,7 @@ class Database(object):
         lumin_agn: array of float
             Luminosity density of the central AGN at each wavelength in W/nm.
         lumin_inrin_agn: array of float
-            Luminosity density of the intrinsic (de-reddened) central AGN 
+            Luminosity density of the intrinsic (de-reddened) central AGN
             at each wavelength in W/nm.
 
 
@@ -825,6 +826,8 @@ class Database(object):
             Luminosity of the accretion disk in W/nm
         dust: array of float
             Luminosity of the dust in W/nm
+        intrin_iso_disk: array of float
+            Luminosity of the intrinsic isotropic disk emission in W/nm
 
         Returns
         -------
@@ -848,7 +851,7 @@ class Database(object):
         if result:
             return SKIRTOR2016(result.t, result.pl, result.q, result.oa,
                                result.R, result.Mcl, result.i, result.wave,
-                               result.disk, result.dust)
+                               result.disk, result.dust, result.intrin_iso_disk)
         else:
             raise DatabaseLookupError(
                 "The SKIRTOR2016 model is not in the database.")

pcigale/data/skirtor2016.py

@@ -6,7 +6,7 @@ class SKIRTOR2016(object):
 
     """
 
-    def __init__(self, t, pl, q, oa, R, Mcl, i, wave, disk, dust):
+    def __init__(self, t, pl, q, oa, R, Mcl, i, wave, disk, dust, intrin_iso_disk):
         """Create a new AGN model. The descriptions of the parameters are taken
            directly from https://sites.google.com/site/skirtorus/sed-library.
 
@@ -37,7 +37,8 @@ class SKIRTOR2016(object):
             Luminosity of the accretion disk in W/nm
         dust: array of float
             Luminosity of the dust in W/nm
-
+        intrin_iso_disk: array of float
+            Luminosity of the intrinsic isotropic disk emission in W/nm
         """
 
         self.t = t
@@ -50,3 +51,4 @@ class SKIRTOR2016(object):
         self.wave = wave
         self.disk = disk
         self.dust = dust
+        self.intrin_iso_disk = intrin_iso_disk

pcigale/sed_modules/fritz2006.py

@@ -166,6 +166,7 @@ class Fritz2006(SedModule):
         self.fritz2006.lumin_therm = lumin_therm_new
         self.fritz2006.lumin_scatt = lumin_scatt_new
         self.fritz2006.lumin_agn = lumin_agn_new
+        self.fritz2006.lumin_intrin_agn /= norm
 
         # Integrate AGN luminosity for different components
         self.l_agn_scatt = np.trapz(self.fritz2006.lumin_scatt, x=self.fritz2006.wave)

pcigale/sed_modules/skirtor2016.py

@@ -17,6 +17,7 @@ import numpy as np
 from pcigale.data import Database
 from . import SedModule
 
+import scipy.constants as cst
 
 class SKIRTOR2016(SedModule):
     """SKIRTOR 2016 (Stalevski et al., 2016) AGN dust torus emission
@@ -82,6 +83,21 @@ class SKIRTOR2016(SedModule):
             'cigale_list(minvalue=0., maxvalue=1.)',
             "AGN fraction.",
             0.1
+        )),
+        ('EBV', (
+            'cigale_list(minvalue=0.)',
+            "E(B-V) for extinction in polar direction",
+            0.1
+        )),
+        ('temperature', (
+            'cigale_list(minvalue=0.)',
+            "Temperature of the dust in K",
+            100.
+        )),
+        ("emissivity", (
+            "cigale_list(minvalue=0.)",
+            "Emissivity index of the dust.",
+            1.6
         ))
     ])
 
@@ -95,6 +111,9 @@ class SKIRTOR2016(SedModule):
         self.Mcl = float(self.parameters["Mcl"])
         self.i = int(self.parameters["i"])
         self.fracAGN = float(self.parameters["fracAGN"])
+        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,
@@ -109,6 +128,54 @@ class SKIRTOR2016(SedModule):
             self.SKIRTOR2016.disk*= (4*np.sin(self.oa*np.pi/180)**2 + 3*np.sin(self.oa*np.pi/180)) / \
                                     (12*np.cos(self.i*np.pi/180)**2 + 6*np.cos(self.i *np.pi/180))
 
+        # Apply polar-dust obscuration
+        # We define various constants necessary to compute the model
+        self.c = cst.c * 1e9
+        lambda_0 = 200e3
+        # Calculate the extinction (SMC)
+        k_lam = 1.39*(self.SKIRTOR2016.wave/1e3)**-1.2-0.38
+        A_lam = k_lam * self.EBV
+        # Set negative value to zero
+        A_lam[A_lam<0] = 0
+        # The extinction factor, flux_1/flux_0
+        ext_fac = 10**(A_lam/-2.5)
+        # Calculate the new AGN SED shape after extinction
+        if self.i <= (90-self.oa):
+            # The direct and scattered components (line-of-sight) are extincted for type-1 AGN
+            disk_new = self.SKIRTOR2016.disk * ext_fac
+        else:
+            # Keep the direct and scatter components for type-2
+            disk_new = self.SKIRTOR2016.disk
+        # Calculate the total extincted luminosity averaged over all directions
+        l_ext = np.trapz(self.SKIRTOR2016.intrin_iso_disk * (1-ext_fac),
+                         x=self.SKIRTOR2016.wave) * \
+                        (1 - np.sin( np.deg2rad(self.oa) ))
+        # Casey (2012) modified black body model
+        conv = self.c / (self.SKIRTOR2016.wave * self.SKIRTOR2016.wave)
+        blackbody = (conv * (1. - np.exp(-(lambda_0 / self.SKIRTOR2016.wave)
+                    ** self.emissivity)) * (self.c / self.SKIRTOR2016.wave) ** 3. / (np.exp(
+                    cst.h * self.c / (self.SKIRTOR2016.wave * cst.k * self.temperature)) - 1.))
+        blackbody *= l_ext / np.trapz(blackbody, x=self.SKIRTOR2016.wave)
+        # Add the black body to dust thermal emission
+        dust_new = self.SKIRTOR2016.dust + blackbody
+        # Normalize direct, scatter, and thermal components
+        norm = np.trapz(dust_new, x=self.SKIRTOR2016.wave)
+        dust_new /= norm
+        disk_new /= norm
+        # Update SKIRTOR model SED
+        self.SKIRTOR2016.dust = dust_new
+        self.SKIRTOR2016.disk = disk_new
+        self.SKIRTOR2016.intrin_iso_disk /= norm
+
+        # Integrate AGN luminosity for different components
+        self.lumin_disk = np.trapz(self.SKIRTOR2016.disk, x=self.SKIRTOR2016.wave)
+        # Intrinsic (de-reddened) AGN luminosity from the central source
+        self.lumin_intrin_disk = np.trapz(self.SKIRTOR2016.intrin_iso_disk,
+                                          x=self.SKIRTOR2016.wave)
+        # Calculate L_lam(2500A)
+        self.l_agn_2500A = np.interp(250, self.SKIRTOR2016.wave, self.SKIRTOR2016.intrin_iso_disk)
+        # Convert L_lam to L_nu
+        self.l_agn_2500A *= 250**2/self.c
 
     def process(self, sed):
         """Add the IR re-emission contributions
@@ -133,14 +200,18 @@ class SKIRTOR2016(SedModule):
         sed.add_info('agn.Mcl', self.Mcl)
         sed.add_info('agn.i', self.i)
         sed.add_info('agn.fracAGN', self.fracAGN)
+        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.:
             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_disk = agn_power * self.lumin_disk
+            lumin_intrin_disk = agn_power * self.lumin_intrin_disk
             lumin = lumin_dust + lumin_disk
+            l_agn_2500A = agn_power * self.l_agn_2500A
         else:
             raise Exception("AGN fraction is exactly 1. Behaviour "
                             "undefined.")
@@ -148,6 +219,8 @@ class SKIRTOR2016(SedModule):
         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_info('agn.intrin_disk_luminosity', lumin_intrin_disk, True)
+        sed.add_info('agn.agn_intrin_Lnu_2500A', l_agn_2500A, True)
 
         sed.add_contribution('agn.SKIRTOR2016_dust', self.SKIRTOR2016.wave,
                              agn_power * self.SKIRTOR2016.dust)