Merge branch 'master' of https://github.com/cosimoNigro/agnpy

agnpy/absorption.py

@@ -1,8 +1,8 @@
 import numpy as np
-from astropy.constants import h, c, m_e, sigma_T
+from astropy.constants import c, G, h, m_e, M_sun, sigma_T
 import astropy.units as u
 from .compton import cos_psi, x_re_shell, x_re_ring, mu_star
-from .targets import SSDisk, SphericalShellBLR, RingDustTorus
+from .targets import PointSourceBehindJet, SSDisk, SphericalShellBLR, RingDustTorus
 
 
 mec2 = m_e.to("erg", equivalencies=u.mass_energy())
@@ -12,7 +12,7 @@
 ]
 
 
-__all__ = ["sigma", "Absorption"]
+__all__ = ["sigma", "tau_disk_finke_2016", "Absorption"]
 
 
 def sigma(s):
@@ -26,6 +26,38 @@ def sigma(s):
     return values
 
 
+def tau_disk_finke_2016(nu, blob, disk, r):
+    """Eq. 63 in Finke 2016, **for testing purposes only**
+    It assumes the blob moves parallel to the jet axis, mu_s = 1 => cos_psi = mu
+    """
+    r_tilde = (r / disk.R_g).to_value("")
+    l_tilde = np.logspace(0, 5, 50) * r_tilde
+    epsilon_1 = nu.to("", equivalencies=epsilon_equivalency)
+    epsilon_1 *= 1 + blob.z
+    # axis 0: R_tilde
+    # axis 1: l_tilde
+    # axis 2: epsilon
+    _R_tilde = np.reshape(disk.R_tilde, (disk.R_tilde.size, 1, 1))
+    _l_tilde = np.reshape(l_tilde, (1, l_tilde.size, 1))
+    _epsilon_1 = np.reshape(epsilon_1, (1, 1, epsilon_1.size))
+    _epsilon = disk.epsilon(_R_tilde)
+    _phi_disk = disk.phi_disk(_R_tilde)
+    _mu = np.sqrt(1 / (1 + _R_tilde ** 2 / _l_tilde ** 2))
+    _s = epsilon_1 * _epsilon * (1 - _mu) / 2
+    integrand = (
+        _l_tilde ** (-2)
+        * _R_tilde ** (-5 / 4)
+        * _phi_disk
+        * (1 + _R_tilde ** 2 / _l_tilde ** 2) ** (-3 / 2)
+        * (sigma(_s) / sigma_T).to_value("")
+        * (1 - _mu)
+    )
+    integral_R_tilde = np.trapz(integrand, disk.R_tilde, axis=0)
+    integral_l_tilde = np.trapz(integral_R_tilde, l_tilde, axis=0)
+    prefactor = 1e7 * disk.l_Edd ** (3 / 4) * disk.M_8 ** (1 / 4) * disk.eta ** (-3 / 4)
+    return prefactor * integral_l_tilde
+
+
 class Absorption:
     """class to compute the absorption due to gamma-gamma pair production
 
@@ -55,75 +87,88 @@ def set_phi(self, phi_size=50):
         self.phi_size = phi_size
         self.phi = np.linspace(0, 2 * np.pi, self.phi_size)
 
-    def set_l(self, l_size=10):
-        """set the range of integration for the distance
-        """
-        # integrate up 3000 pc
-        self.l_size = l_size
-        l_max = 3000 * u.pc
+    def set_l(self, l_size=50):
+        """set an array of distances over which to integrate
+        integrate up to 100 kpc"""
+        max_l = 100 * u.kpc
         self.l = (
             np.logspace(
-                np.log10(self.r.to_value("cm")),
-                np.log10(l_max.to_value("cm")),
-                self.l_size,
+                np.log10(self.r.to_value("cm")), np.log10(max_l.to_value("cm")), l_size
             )
             * u.cm
         )
 
-    def _opacity_disk(self, nu):
-        """opacity generated by a Shakura Sunyaev disk
+    def _opacity_point_source(self, nu):
+        """opacity generated by a point source behind the jet
 
         Parameters
         ----------
         nu : `~astropy.units.Quantity`
-            array of frequencies, in Hz, to compute the sed, **note** these are
+            array of frequencies, in Hz, to compute the sed, **note** these are 
             observed frequencies (observer frame).
         """
         # define the dimensionless energy
         epsilon_1 = nu.to("", equivalencies=epsilon_equivalency)
         # transform to BH frame
         epsilon_1 *= 1 + self.blob.z
+        # multidimensional integration (only in solid angle)
+        s = self.target.epsilon_0 * epsilon_1 * (1 - self.blob.mu_s) / 2
+        integral = (1 - self.blob.mu_s) * sigma(s) / self.r
+        prefactor = self.target.L_0 / (4 * np.pi * self.target.epsilon_0 * m_e * c ** 3)
+        tau = prefactor * integral
+        return tau.to_value("")
+
+    def _opacity_disk(self, nu):
+        """opacity generated by a Shakura Sunyaev disk
+        Parameters
+        ----------
+        nu : `~astropy.units.Quantity`
+            array of frequencies, in Hz, to compute the absorption, **note** 
+            these are observed frequencies (observer frame).
+        """
+        # define the dimensionless energy
+        epsilon_1 = nu.to("", equivalencies=epsilon_equivalency)
+        # transform to BH frame
+        epsilon_1 *= 1 + self.blob.z
         # each value of l, distance from the BH, defines a different range of
-        # cosine integration, we have to break the integration as the array of
-        # mu takes different values at each distance
-        integral = np.empty(len(epsilon_1))
-        for i, _epsilon_1 in enumerate(epsilon_1):
-            integrand_l = np.empty(len(self.l))
-            for j, _l in enumerate(self.l):
-                print("j: ", j)
-                print("l: ", _l)
-                l_tilde = (_l / self.target.R_g).to_value("")
-                # for the multidimensional integration on the angles only
-                # axis 0 : phi
-                # axis 1 : mu
-                _phi = np.reshape(self.phi, (1, self.phi.size))
-                mu = self.target.mu_from_r_tilde(l_tilde)
-                print("mu: ", mu)
-                print("-------------")
+        # cosine integration, we cannot use direct integration of multidimensional
+        # numpy arrays but we have to use an explicit loop over the distances
+        # as mu takes a different array of values at each distance
+        l_tilde = np.logspace(0, 5, 50) * (self.r / self.target.R_g).to_value("")
+        tau_epsilon_1 = np.ndarray(0)
+        for _epsilon_1 in epsilon_1:
+            integrand_l_tilde = np.ndarray(0)
+            for _l_tilde in l_tilde:
+                # multidimensional integration (only in solid angle)
+                # axis 0: mu
+                # axis 1: phi
+                mu = self.target.mu_from_r_tilde(_l_tilde)
                 _mu = np.reshape(mu, (mu.size, 1))
-                # epsilon and phi disk have the same dimension as mu
-                _epsilon = self.target.epsilon_mu(_mu, l_tilde)
-                _phi_disk_mu = self.target.phi_disk_mu(_mu, l_tilde)
+                _phi = np.reshape(self.phi, (1, self.phi.size))
+                _epsilon = self.target.epsilon_mu(_mu, _l_tilde)
+                _phi_disk = self.target.phi_disk_mu(_mu, _l_tilde)
+                # angle between the photons
                 _cos_psi = cos_psi(self.blob.mu_s, _mu, _phi)
                 _s = _epsilon_1 * _epsilon * (1 - _cos_psi) / 2
-                _integrand_mu = _phi_disk_mu / (
-                    _epsilon
-                    * np.power(_l, 3)
-                    * _mu
-                    * np.power(np.power(_mu, -2) - 1, 3 / 2)
+                _cross_section = sigma(_s).to_value("cm2")
+                _integrand = (
+                    _phi_disk
+                    / (_epsilon * _mu * _l_tilde ** 3 * (_mu ** (-2) - 1) ** (-3 / 2))
+                    * (1 - _cos_psi)
+                    * _cross_section
                 )
-                _integrand = _integrand_mu * (1 - _cos_psi) * sigma(_s)
-                # integrate over mu and phi
-                integral_mu = np.trapz(_integrand, mu, axis=0)
-                integral_phi = np.trapz(integral_mu, self.phi, axis=0)
-                integrand_l[i] = integral_phi.to_value("cm-1")
-            # integrate over l
-            integral[j] = np.trapz(integrand_l, self.l, axis=0).to_value("cm")
-
-        prefactor_num = 3 * self.target.L_disk * self.target.R_g
-        prefactor_denum = 16 * np.pi * self.target.eta * m_e * np.power(c, 3)
-        tau = prefactor_num / prefactor_denum * integral
-        return tau.to("")
+                # integrate over the solid angle
+                _integral_mu = np.trapz(_integrand, mu, axis=0)
+                _integral_phi = np.trapz(_integral_mu, self.phi, axis=0)
+                integrand_l_tilde = np.append(integrand_l_tilde, _integral_phi)
+            # integrate over the distance
+            integral_l_tilde = np.trapz(integrand_l_tilde, l_tilde, axis=0)
+            tau_epsilon_1 = np.append(tau_epsilon_1, integral_l_tilde)
+        prefactor = (3 * self.target.L_disk) / (
+            (4 * np.pi) ** 2 * self.target.eta * m_e * c ** 3 * self.target.R_g
+        )
+        tau_epsilon_1 *= prefactor.to_value("cm-2")
+        return tau_epsilon_1
 
     def _opacity_shell_blr(self, nu):
         """opacity generated by a spherical shell Broad Line Region
@@ -148,24 +193,20 @@ def _opacity_shell_blr(self, nu):
         _phi = np.reshape(self.phi, (1, self.phi.size, 1, 1))
         _l = np.reshape(self.l, (1, 1, self.l.size, 1))
         _epsilon_1 = np.reshape(epsilon_1, (1, 1, 1, epsilon_1.size))
-        # define integrating function
         _x = x_re_shell(_mu, self.target.R_line, _l)
         _mu_star = mu_star(_mu, self.target.R_line, _l)
-
         _cos_psi = cos_psi(self.blob.mu_s, _mu_star, _phi)
         _s = _epsilon_1 * self.target.epsilon_line * (1 - _cos_psi) / 2
-        _integrand = (1 - _cos_psi) * np.power(_x, -2) * sigma(_s)
-
-        prefactor_num = self.target.xi_line * self.target.L_disk
-        prefactor_denum = (
-            np.power(4 * np.pi, 2) * self.target.epsilon_line * m_e * np.power(c, 3)
-        )
-
+        # define integrating function
+        _integrand = (1 - _cos_psi) / _x ** 2 * sigma(_s)
+        # integrate
         integral_mu = np.trapz(_integrand, self.mu, axis=0)
         integral_phi = np.trapz(integral_mu, self.phi, axis=0)
         integral = np.trapz(integral_phi, self.l, axis=0)
-
-        tau = prefactor_num / prefactor_denum * integral
+        prefactor = (self.target.xi_line * self.target.L_disk) / (
+            (4 * np.pi) ** 2 * self.target.epsilon_line * m_e * c ** 3
+        )
+        tau = prefactor * integral
         return tau.to_value("")
 
     def _opacity_ring_torus(self, nu):
@@ -191,20 +232,18 @@ def _opacity_ring_torus(self, nu):
         _epsilon_1 = np.reshape(epsilon_1, (1, 1, epsilon_1.size))
         _x = x_re_ring(self.target.R_dt, _l)
         _mu = _l / _x
-
+        # define the
         _cos_psi = cos_psi(self.blob.mu_s, _mu, _phi)
         _s = _epsilon_1 * self.target.epsilon_dt * (1 - _cos_psi) / 2
-        _integrand = (1 - _cos_psi) * np.power(_x, -2) * sigma(_s)
-
-        prefactor_num = self.target.xi_dt * self.target.L_disk
-        prefactor_denum = (
-            np.power(4 * np.pi, 2) * self.target.epsilon_dt * m_e * np.power(c, 3)
-        )
-
+        # define integrating function
+        _integrand = (1 - _cos_psi) / _x ** 2 * sigma(_s)
+        # integrate
         integral_phi = np.trapz(_integrand, self.phi, axis=0)
         integral = np.trapz(integral_phi, self.l, axis=0)
-
-        tau = prefactor_num / prefactor_denum * integral
+        prefactor = (self.target.xi_dt * self.target.L_disk) / (
+            8 * np.pi ** 2 * self.target.epsilon_dt * m_e * c ** 3
+        )
+        tau = prefactor * integral
         return tau.to_value("")
 
     def tau(self, nu):
@@ -219,6 +258,8 @@ def tau(self, nu):
             array of frequencies, in Hz, to compute the opacity, **note** these are
             observed frequencies (observer frame).
         """
+        if isinstance(self.target, PointSourceBehindJet):
+            return self._opacity_point_source(nu)
         if isinstance(self.target, SSDisk):
             return self._opacity_disk(nu)
         if isinstance(self.target, SphericalShellBLR):

agnpy/compton.py

@@ -335,7 +335,7 @@ def _sed_flux_cmb(self, nu):
         sed = prefactor_num / prefactor_denom * integral_phi
         return sed.to("erg cm-2 s-1")
 
-    def _sed_flux_point_like(self, nu):
+    def _sed_flux_point_source(self, nu):
         """SED flux for EC on a point like source behind the jet
         
         Parameters
@@ -556,7 +556,7 @@ def sed_flux(self, nu):
         if isinstance(self.target, CMB):
             return self._sed_flux_cmb(nu)
         if isinstance(self.target, PointSourceBehindJet):
-            return self._sed_flux_point_like(nu)
+            return self._sed_flux_point_source(nu)
         if isinstance(self.target, SSDisk):
             return self._sed_flux_disk(nu)
         if isinstance(self.target, SphericalShellBLR):