solve conflict

.gitignore

@@ -6,5 +6,8 @@ pcigale/data/data.db
 pcigale.egg-info
 docs/_build
 docs/api
+<<<<<<< HEAD
 dist/
+=======
+>>>>>>> 530a13b71a7b3350e90220c17c299cec51dd966a
 *.swp

pcigale/analysis_modules/pdf_analysis/utils.py

@@ -266,6 +266,18 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
         nan_idxs = agn_idxs[ np.abs(det_alpha_ox) > models.intprop['xray.max_dev_alpha_ox'][wz][agn_idxs] ]
         chi2[nan_idxs] = np.nan
 
+    # Penalize det_alpha_ox which lie out of the user-set range
+    if 'xray' in models.params.modules:
+        # Get the model indice that have valid AGN component
+        agn_idxs = np.where( models.extprop['agn.intrin_Lnu_2500A'][wz]>0 )[0]
+        # Calculate expected alpha_ox from Lnu_2500 (Just et al. 2007)
+        exp_alpha_ox = -0.137*np.log10(models.extprop['agn.intrin_Lnu_2500A'][wz][agn_idxs]*1e7*scaling[agn_idxs])+2.638
+        # Calculate det_alpha_ox = alpha_ox - alpha_ox(Lnu_2500)
+        det_alpha_ox = models.intprop['xray.alpha_ox'][wz][agn_idxs] - exp_alpha_ox
+        # If det_alpha_ox out of range, set corresponding chi2 to nan
+        nan_idxs = agn_idxs[ np.abs(det_alpha_ox) > models.intprop['xray.max_dev_alpha_ox'][wz][agn_idxs] ]
+        chi2[nan_idxs] = np.nan
+
     # Computation of the χ² from intensive properties
     for name, prop in obs.intprop.items():
         model = models.intprop[name][wz]

pcigale/sed_modules/fritz2006.py

@@ -253,6 +253,78 @@ class Fritz2006(SedModule):
         # Convert L_lam to L_nu
         self.l_agn_2500A *= 250**2/self.c
 
+        # Apply wavelength cut to avoid X-ray wavelength
+        lam_cut = 10**0.9
+        lam_idxs = self.fritz2006.wave>=lam_cut
+        # Calculate the re-normalization factor to keep energy conservation
+        norm_fac = np.trapz(self.fritz2006.lumin_intrin_agn, x=self.fritz2006.wave) /\
+            np.trapz(self.fritz2006.lumin_intrin_agn[lam_idxs], x=self.fritz2006.wave[lam_idxs])
+        # Perform the cut
+        self.fritz2006.wave = self.fritz2006.wave[lam_idxs]
+        self.fritz2006.lumin_agn = self.fritz2006.lumin_agn[lam_idxs]*norm_fac
+        self.fritz2006.lumin_scatt = self.fritz2006.lumin_scatt[lam_idxs]*norm_fac
+        self.fritz2006.lumin_therm  = self.fritz2006.lumin_therm[lam_idxs]
+        self.fritz2006.lumin_intrin_agn = self.fritz2006.lumin_intrin_agn[lam_idxs]*norm_fac
+
+        # 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)
+        # The analytical formula is from Bongiorno+2012
+        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)
+        # Calculate the new AGN SED shape after extinction
+        if self.psy > (90-self.opening_angle):
+            # The direct and scattered components (line-of-sight) are extincted for type-1 AGN
+            lumin_agn_new = self.fritz2006.lumin_agn * ext_fac
+            lumin_scatt_new = self.fritz2006.lumin_scatt * ext_fac
+        else:
+            # Keep the direct and scatter components for type-2
+            lumin_agn_new = self.fritz2006.lumin_agn
+            lumin_scatt_new = self.fritz2006.lumin_scatt
+        # Calculate the total extincted luminosity averaged over all directions
+        # note that self.opening_angle is different from Fritz's definiting!
+        l_ext = np.trapz(self.fritz2006.lumin_intrin_agn * (1-ext_fac),
+                         x=self.fritz2006.wave) * \
+                        (1 - np.cos( np.deg2rad(self.opening_angle) ))
+        # Casey (2012) modified black body model
+        conv = self.c / (self.fritz2006.wave * self.fritz2006.wave)
+        # To avoid inf occurance in exponential, set blackbody=0 when h*c/lam*k*T is large
+        hc_lkt = cst.h * self.c / (self.fritz2006.wave * cst.k * self.temperature)
+        non_0_idxs = hc_lkt<100
+        # Generate the blackbody
+        lumin_blackbody = np.zeros(len(self.fritz2006.wave))
+        lumin_blackbody[non_0_idxs] = conv[non_0_idxs] * \
+            (1. - np.exp(-(lambda_0 / self.fritz2006.wave[non_0_idxs])** self.emissivity)) * \
+            (self.c / self.fritz2006.wave[non_0_idxs]) ** 3. / \
+            (np.exp(hc_lkt[non_0_idxs]) - 1.)
+        lumin_blackbody *= l_ext / np.trapz(lumin_blackbody, x=self.fritz2006.wave)
+        # Add the black body to dust thermal emission
+        lumin_therm_new = self.fritz2006.lumin_therm + lumin_blackbody
+        # Normalize direct, scatter, and thermal components
+        norm = np.trapz(lumin_therm_new, x=self.fritz2006.wave)
+        lumin_therm_new /= norm
+        lumin_scatt_new /= norm
+        lumin_agn_new /= norm
+        # Update fritz2006 lumin
+        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)
+        self.l_agn_agn = np.trapz(self.fritz2006.lumin_agn, x=self.fritz2006.wave)
+        # Intrinsic (de-reddened) AGN luminosity from the central source
+        self.l_agn_intrin_agn = np.trapz(self.fritz2006.lumin_intrin_agn, x=self.fritz2006.wave)
+        # Calculate L_lam(2500A)
+        self.l_agn_2500A = np.interp(250, self.fritz2006.wave, self.fritz2006.lumin_intrin_agn)
+        # Convert L_lam to L_nu
+        self.l_agn_2500A *= 250**2/self.c
+
     def process(self, sed):
         """Add the IR re-emission contributions