allow several extinction laws for polar dust

pcigale/sed_modules/fritz2006.py

@@ -19,6 +19,56 @@ from . import SedModule
 
 import scipy.constants as cst
 
+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/1e3)**-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 = 1 / (wavelength/1e3)
+        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 Fritz2006(SedModule):
     """Fritz et al. (2006) AGN dust torus emission
 
@@ -85,6 +135,12 @@ class Fritz2006(SedModule):
             "AGN fraction.",
             0.1
         )),
+        ('law', (
+            'cigale_list(dtype=int, options=0 & 1 & 2)',
+            "The extinction law of polar dust: "
+            "0 (SMC), 1 (Calzetti 2000), or 2 (Gaskell et al. 2004)",
+            0
+        )),
         ('EBV', (
             'cigale_list(minvalue=0.)',
             "E(B-V) for extinction in polar direction",
@@ -97,7 +153,7 @@ class Fritz2006(SedModule):
         )),
         ("emissivity", (
             "cigale_list(minvalue=0.)",
-            "Emissivity index of the polar dust.",
+            "Emissivity index of the polar dust",
             1.6
         ))
    ])
@@ -111,6 +167,7 @@ class Fritz2006(SedModule):
         self.opening_angle = (180. - self.parameters["opening_angle"]) / 2.
         self.psy = float(self.parameters["psy"])
         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"])
@@ -143,7 +200,7 @@ class Fritz2006(SedModule):
         lambda_0 = 200e3
         # Calculate the extinction (SMC)
         # The analytical formula is from Bongiorno+2012
-        k_lam = 1.39*(self.fritz2006.wave/1e3)**-1.2
+        k_lam = k_ext(self.fritz2006.wave, self.law)
         A_lam = k_lam * self.EBV
         # The extinction factor, flux_1/flux_0
         ext_fac = 10**(A_lam/-2.5)

pcigale/sed_modules/skirtor2016.py

@@ -19,6 +19,56 @@ from . import SedModule
 
 import scipy.constants as cst
 
+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/1e3)**-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 = 1 / (wavelength/1e3)
+        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
 
@@ -84,6 +134,12 @@ class SKIRTOR2016(SedModule):
             "AGN fraction.",
             0.1
         )),
+        ('law', (
+            'cigale_list(dtype=int, options=0 & 1 & 2)',
+            "The extinction law of polar dust: "
+            "0 (SMC), 1 (Calzetti 2000), or 2 (Gaskell et al. 2004)",
+            0
+        )),
         ('EBV', (
             'cigale_list(minvalue=0.)',
             "E(B-V) for extinction in polar direction",
@@ -96,7 +152,7 @@ class SKIRTOR2016(SedModule):
         )),
         ("emissivity", (
             "cigale_list(minvalue=0.)",
-            "Emissivity index of the polar dust.",
+            "Emissivity index of the polar dust",
             1.6
         ))
     ])
@@ -111,6 +167,7 @@ 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"])
@@ -124,9 +181,9 @@ class SKIRTOR2016(SedModule):
         # We define various constants necessary to compute the model
         self.c = cst.c * 1e9
         lambda_0 = 200e3
-        # Calculate the extinction (SMC)
-        # The analytical formula is from Bongiorno+2012
-        k_lam = 1.39*(self.SKIRTOR2016.wave/1e3)**-1.2
+        # Calculate the extinction
+        k_lam = k_ext(self.SKIRTOR2016.wave, self.law)
+        1.39*(self.SKIRTOR2016.wave/1e3)**-1.2
         A_lam = k_lam * self.EBV
         # The extinction factor, flux_1/flux_0
         ext_fac = 10**(A_lam/-2.5)