change disk SED shape in SKIRTOR; move some process to data-reading stage to improve efficiency

database_builder/__init__.py

@@ -733,31 +733,90 @@ def build_skirtor2016(base):
 
     # 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))
+    wl_uncut, intrin_disk_unnormed = np.genfromtxt(filename, unpack=True, usecols=(0, 6))
+    # Convert to nm units
+    wl_uncut *= 1e3
+    # Convert from F_lam*lam to F_lam
+    intrin_disk_unnormed /= wl_uncut
+
+    # Define the new disk SED shape
+    # Modify the disk SED to Feltre et al. (2012), following the method in
+    # https://sites.google.com/site/skirtorus/sed-library
+    # Set the broken power law shape
+    # 1st part: < 50 nm
+    disk_new_unnormed = wl_uncut.copy()
+    # 2nd part: 50-125 nm
+    wave_idxs = (wl_uncut>=50) & (wl_uncut<125)
+    # Here, the factor 50**0.2 is to re-nomalize for continuity
+    norm_fac = 50**0.2
+    disk_new_unnormed[wave_idxs] = disk_new_unnormed[wave_idxs]**0.8 * norm_fac
+    # 3rd part: 125-1e4 nm
+    norm_fac *= 125**2.3
+    wave_idxs = (wl_uncut>=125) & (wl_uncut<1e4)
+    disk_new_unnormed[wave_idxs] = disk_new_unnormed[wave_idxs]**-1.5 * norm_fac
+    # 4th part: >1e4 nm
+    norm_fac *= 1e4**2.5
+    wave_idxs = wl_uncut>=1e4
+    disk_new_unnormed[wave_idxs] = disk_new_unnormed[wave_idxs]**-4 * norm_fac
+    # Normalize the new disk SED, so that its total energy is 1
+    disk_new_unnormed /= np.trapz(disk_new_unnormed, x=wl_uncut)
+
+    # Change the SED shape of intrinsic disk luminosity
+    intrin_disk_unnormed = disk_new_unnormed * np.trapz(intrin_disk_unnormed, x=wl_uncut)
     # Convert anisotropic to the isotropic disk luminosity
-    intrin_iso_disk_unnormed = intrin_disk_unnormed*7/18
+    intrin_disk_unnormed *= 7/18
 
     for params in iter_params:
         filename = skirtor2016_dir + \
                 "t{}_p{}_q{}_oa{}_R{}_Mcl{}_i{}_sed.dat".format(*params)
         print("Importing {}...".format(filename))
+        # Convert some useful parameters to float
+        i_float, oa_float = float(params[6]), float(params[3])
 
-        wl, disk, scatt, dust = np.genfromtxt(filename, unpack=True,
-                                              usecols=(0, 2, 3, 4))
-        wl *= 1e3
-        disk += scatt
-        disk /= wl
-        dust /= wl
+        disk, scatt, dust = np.genfromtxt(filename, unpack=True, usecols=(2, 3, 4))
+        disk /= wl_uncut
+        scatt /= wl_uncut
+        dust /= wl_uncut
 
         # Normalization of the lumin_therm to 1W
-        norm = np.trapz(dust, x=wl)
+        norm = np.trapz(dust, x=wl_uncut)
         disk /= norm
+        scatt /= norm
         dust /= norm
-        intrin_iso_disk = intrin_iso_disk_unnormed / (norm*wl)
+        intrin_disk = intrin_disk_unnormed / norm
+
+        if i_float<=(90-oa_float):
+            # For type-1 AGN, update the disk emission
+            # Re-normalize the entire new disk SED to keep energy conservation
+            disk = np.trapz(disk, x=wl_uncut) * disk_new_unnormed
+
+        # Add the scatter component to disk for simplicity
+        disk += scatt
+
+        # Apply wavelength cut to avoid X-ray wavelength
+        lam_cut = 5.
+        lam_idxs = wl_uncut>=lam_cut
+        # Calculate the re-normalization factor to keep energy conservation
+        norm_fac = np.trapz(intrin_disk, x=wl_uncut) / \
+                   np.trapz(intrin_disk[lam_idxs], x=wl_uncut[lam_idxs])
+        # Perform the cut
+        wl = wl_uncut[lam_idxs]
+        disk = disk[lam_idxs]*norm_fac
+        dust = dust[lam_idxs]
+        intrin_disk = intrin_disk[lam_idxs]*norm_fac
+
+        # Modify the model so that the central source is isotropic
+        # Calculate the coverting factors disk emission
+        # (dust emission is normalized to unity, so it remains same)
+        disk_fac = 7/(12*np.cos(i_float*np.pi/180)**2 + 6*np.cos(i_float*np.pi/180))
+        dust_fac = 7/(4*np.sin(oa_float*np.pi/180)**2 + 3*np.sin(oa_float*np.pi/180))
+        # Only apply disk_fac to type 1 (for type 2 it is zero)
+        if i_float<=(90-oa_float): disk*= disk_fac/dust_fac
+        intrin_disk /= dust_fac
 
         models.append(SKIRTOR2016(params[0], params[1], params[2], params[3],
                                   params[4], params[5], params[6], wl, disk,
-                                  dust, intrin_iso_disk))
+                                  dust, intrin_disk))
 
     base.add_skirtor2016(models)
 

pcigale/data/__init__.py

@@ -220,7 +220,7 @@ class _SKIRTOR2016(BASE):
     wave = Column(PickleType)
     disk = Column(PickleType)
     dust = Column(PickleType)
-    intrin_iso_disk = Column(PickleType)
+    intrin_disk = Column(PickleType)
 
     def __init__(self, agn):
         self.t = agn.t
@@ -233,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
+        self.intrin_disk = agn.intrin_disk
 
 class _NebularLines(BASE):
     """Storage for line templates
@@ -826,7 +826,7 @@ 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
+        intrin_disk: array of float
             Luminosity of the intrinsic isotropic disk emission in W/nm
 
         Returns
@@ -851,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.intrin_iso_disk)
+                               result.disk, result.dust, result.intrin_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, intrin_iso_disk):
+    def __init__(self, t, pl, q, oa, R, Mcl, i, wave, disk, dust, intrin_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,7 @@ 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
+        intrin_disk: array of float
             Luminosity of the intrinsic isotropic disk emission in W/nm
         """
 
@@ -51,4 +51,4 @@ class SKIRTOR2016(object):
         self.wave = wave
         self.disk = disk
         self.dust = dust
-        self.intrin_iso_disk = intrin_iso_disk
+        self.intrin_disk = intrin_disk

pcigale/sed_modules/skirtor2016.py

@@ -120,51 +120,6 @@ class SKIRTOR2016(SedModule):
                                                     self.oa, self.R, self.Mcl,
                                                     self.i)
 
-        # Modify the disk SED to Feltre et al. (2012), following the method in
-        # https://sites.google.com/site/skirtorus/sed-library
-        # Set the broken power law shape
-        # 1st part: < 50 nm
-        disk_new = self.SKIRTOR2016.wave.copy()
-        # 2nd part: 50-125 nm
-        wave_idxs = (self.SKIRTOR2016.wave>=50) & (self.SKIRTOR2016.wave<125)
-        # Here, the factor 50**0.2 is to re-nomalize for continuity
-        norm_fac = 50**0.2
-        disk_new[wave_idxs] = disk_new[wave_idxs]**0.8 * norm_fac
-        # 3rd part: 125-1e4 nm
-        norm_fac *= 125**2.3
-        wave_idxs = (self.SKIRTOR2016.wave>=125) & (self.SKIRTOR2016.wave<1e4)
-        disk_new[wave_idxs] = disk_new[wave_idxs]**-1.5 * norm_fac
-        # 4th part: >1e4 nm
-        norm_fac *= 1e4**2.5
-        wave_idxs = self.SKIRTOR2016.wave>=1e4
-        disk_new[wave_idxs] = disk_new[wave_idxs]**-4 * norm_fac
-        # Re-normalize the entire disk SED to keep energy conservation
-        disk_new *= np.trapz(self.SKIRTOR2016.disk, x=self.SKIRTOR2016.wave) / \
-                    np.trapz(disk_new, x=self.SKIRTOR2016.wave)
-        self.SKIRTOR2016.disk = disk_new
-
-        # Apply wavelength cut to avoid X-ray wavelength
-        lam_cut = 5.
-        lam_idxs = self.SKIRTOR2016.wave>=lam_cut
-        # Calculate the re-normalization factor to keep energy conservation
-        norm_fac = np.trapz(self.SKIRTOR2016.intrin_iso_disk, x=self.SKIRTOR2016.wave) /\
-            np.trapz(self.SKIRTOR2016.intrin_iso_disk[lam_idxs], x=self.SKIRTOR2016.wave[lam_idxs])
-        # Perform the cut
-        self.SKIRTOR2016.wave = self.SKIRTOR2016.wave[lam_idxs]
-        self.SKIRTOR2016.disk = self.SKIRTOR2016.disk[lam_idxs]*norm_fac
-        self.SKIRTOR2016.dust = self.SKIRTOR2016.dust[lam_idxs]
-        self.SKIRTOR2016.intrin_iso_disk = self.SKIRTOR2016.intrin_iso_disk[lam_idxs]*norm_fac
-
-        # Modify the model so that the central source is isotropic
-        # Calculate the coverting factors disk emission
-        # (dust emission is normalized to unity, so it remains same)
-        disk_fac = 7/(12*np.cos(self.i*np.pi/180)**2 + 6*np.cos(self.i *np.pi/180))
-        dust_fac = 7/(4*np.sin(self.oa*np.pi/180)**2 + 3*np.sin(self.oa*np.pi/180))
-        # Only apply disk_fac to type 1 (for type 2 it is zero)
-        if self.i<=(90-self.oa):
-            self.SKIRTOR2016.disk*= disk_fac/dust_fac
-        self.SKIRTOR2016.intrin_iso_disk /= dust_fac
-
         # Apply polar-dust obscuration
         # We define various constants necessary to compute the model
         self.c = cst.c * 1e9
@@ -184,7 +139,7 @@ class SKIRTOR2016(SedModule):
             # 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),
+        l_ext = np.trapz(self.SKIRTOR2016.intrin_disk * (1-ext_fac),
                          x=self.SKIRTOR2016.wave) * \
                         (1 - np.sin( np.deg2rad(self.oa) ))
         # Casey (2012) modified black body model
@@ -202,15 +157,15 @@ class SKIRTOR2016(SedModule):
         # Update SKIRTOR model SED
         self.SKIRTOR2016.dust = dust_new
         self.SKIRTOR2016.disk = disk_new
-        self.SKIRTOR2016.intrin_iso_disk /= norm
+        self.SKIRTOR2016.intrin_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,
+        self.lumin_intrin_disk = np.trapz(self.SKIRTOR2016.intrin_disk,
                                           x=self.SKIRTOR2016.wave)
         # Calculate L_lam(2500A)
-        self.l_agn_2500A = np.interp(250, self.SKIRTOR2016.wave, self.SKIRTOR2016.intrin_iso_disk)
+        self.l_agn_2500A = np.interp(250, self.SKIRTOR2016.wave, self.SKIRTOR2016.intrin_disk)
         # Convert L_lam to L_nu
         self.l_agn_2500A *= 250**2/self.c