Update docs

Author: bgmeulem

docs/_static/_images/isofluxlines.png

@@ -9,4 +9,17 @@ @ARTICLE{Luminet_1979
         pages = {228-235},
        adsurl = {https://ui.adsabs.harvard.edu/abs/1979A&A....75..228L},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
-}
+}
+
+@ARTICLE{Page_1974,
+  title     = "Disk-accretion onto a black hole. Time-averaged structure of
+               accretion disk",
+  author    = "Page, Don N and Thorne, Kip S",
+  journal   = "Astrophys. J.",
+  publisher = "American Astronomical Society",
+  volume    =  191,
+  pages     = "499",
+  month     =  jul,
+  year      =  1974,
+  language  = "en"
+}

docs/source/bibliography.bib

@@ -270,9 +270,11 @@ def plot_isoradials(
 
         Example::
 
-        >>> direct_irs = [6, 10, 15, 20]
-        >>> ghost_irs = [6, 20, 50, 100]  # ghost_r can go to infinity
-        >>> ax = bh.plot_isoradials(direct_irs, ghost_irs, lw=1, colors='white')
+            from luminet.black_hole import BlackHole
+
+            direct_irs = [6, 10, 15, 20]
+            ghost_irs = [6, 20, 50, 100]  # ghost_r can go to infinity
+            ax = bh.plot_isoradials(direct_irs, ghost_irs, lw=1, colors='white')
 
         .. image:: /../_static/_images/isoradials.png
            :align: center
@@ -319,8 +321,10 @@ def plot(self, **kwargs) -> Axes:
 
         Example::
 
-            >>> bh = BlackHole()
-            >>> bh.plot()
+            from luminet.black_hole import BlackHole
+
+            bh = BlackHole()
+            bh.plot()
 
         .. image:: /../_static/_images/bh.png
            :align: center
@@ -345,10 +349,12 @@ def plot_isoredshifts(self, redshifts=None, order=0, ax=None, **kwargs) -> Axes:
 
         Example::
 
-            >>> bh = BlackHole()
-            >>> redshifts = [-.2, -.1, 0., .1, .2, .3, .4]
-            >>> ax = bh.plot_isoredshifts(redshifts, c='white')
-            >>> ax = bh.disk_apparent_inner_edge.plot(ax=ax, c='white')
+            from luminet.black_hole import BlackHole
+
+            bh = BlackHole()
+            redshifts = [-.2, -.1, 0., .1, .2, .3, .4]
+            ax = bh.plot_isoredshifts(redshifts, c='white')
+            ax = bh.disk_apparent_inner_edge.plot(ax=ax, c='white')
 
         .. image:: /../_static/_images/isoredshifts.png
            :align: center
@@ -362,14 +368,17 @@ def plot_isoredshifts(self, redshifts=None, order=0, ax=None, **kwargs) -> Axes:
             ax = isoredshift.plot(ax, **kwargs)
         return ax
 
-    def plot_isofluxlines(self, mask_inner=True, normalize=True, order=0, ax=None, **kwargs) -> Axes:
+    def plot_isofluxlines(self, mask_inner=True, mask_outer=True, normalize=True, order=0, ax=None, **kwargs) -> Axes:
         """Plot lines of equal flux.
 
         Args:
             normalize (bool): Whether to normalize the fluxlines by the maximum flux or not. Defaults to True.
             mask_inner (bool): 
                 Whether to place a mask over the apparent inner edge, where the direct image produces no flux. 
-                Useful to mitigate matplotlib tricontour artifacts.
+                Useful to mitigate matplotlib tricontour artifacts. Default is ``True``
+            mask_outer (bool): 
+                Whether to place a mask over the apparent outer edge, where we are not capturing photons from.
+                Useful to mitigate matplotlib tricontour artifacts. Default is ``True``.
             order (int): The order of the image to plot siofluxlines for. Default is :math:`0`.
             ax (:class:`~matplotlib.axes.Axes`, optional): Axes object to plot on.
                 Useful for when you want to plot multiple things one a single canvas.
@@ -383,9 +392,11 @@ def plot_isofluxlines(self, mask_inner=True, normalize=True, order=0, ax=None, *
 
         Example::
 
-            >>> bh = BlackHole(incl=1.4, radial_resolution=200)
-            >>> levels = np.logspace(-1.3, 0, 10)
-            >>> ax = bh.plot_isofluxlines(colors='white', levels=levels, linewidths=1)
+            from luminet.black_hole import BlackHole
+
+            bh = BlackHole(incl=1.4, radial_resolution=200)
+            levels = [.05, .1, .15, .2, .25, .3, .6, .9, 1.2, 1.5, 1.8, 2.1]
+            ax = bh.plot_isofluxlines(colors='white', levels=levels, linewidths=1)
 
         .. image:: /../_static/_images/isofluxlines.png
            :align: center
@@ -426,6 +437,15 @@ def plot_isofluxlines(self, mask_inner=True, normalize=True, order=0, ax=None, *
                 color='k',
                 zorder=len(contour.levels) + 1
                 )
+        if mask_outer:
+            max_r = ax.get_ylim()[-1]
+            ax.fill_between(
+                self.disk_apparent_outer_edge.angles,
+                self.disk_apparent_outer_edge.impact_parameters,
+                max_r, 
+                color='k',
+                zorder=len(contour.levels) + 1
+                )
         return ax
 
     def sample_photons(self, n_points=1000) -> Tuple[pd.DataFrame]:

luminet/black_hole.py

@@ -47,7 +47,8 @@ def calc_b_from_perigee(p: float, bh_mass: float) -> float:
         The fraction on the right hand side equals :math:`b^2`, not :math:`b`.
         You can verify this by filling in :math:`u_2` in Equation 3.
         Only this way do the limits :math:`P -> 3M` and :math:`P >> M` hold true,
-        as well as the value for :math:`b_c`
+        as well as the value for :math:`b_c`.  The resulting images of the paper are correct though.
+
 
     Returns:
         float: Impact parameter :math:`b`
@@ -74,7 +75,8 @@ def calc_k(periastron: float, bh_mass: float) -> float:
         float: Modulus of the elliptic integral
 
     Attention:
-        Mind the typo in :cite:t:`Luminet_1979`. The numerator should be in brackets.
+        Mind the typo in :cite:t:`Luminet_1979`. The numerator should be in brackets. The resulting images of the paper are correct though.
+
     """
     q = calc_q(periastron, bh_mass)
     if q is np.nan:
@@ -468,32 +470,73 @@ def ellipse(r, a, incl) -> float:
 
 
 def calc_Z1(bh_mass, a):
-    """Calculate :math:`Z1` for Kerr black holes.
+    r"""Calculate :math:`Z1` for Kerr black holes.
+
+    The variable :math:`Z1` is used to calculate the innermost orbit for Kerr black holes.
+
+    .. math::
+    
+       Z_1 \equiv 1 + \sqrt[3]{1-a_*^2}\left[ \sqrt[3]{1+a_*} + \sqrt[3]{1-a_*} \right]
+
+    Args:
+        bh_mass (float): Mass of the black hole.
+        a (float): Specific angular momentum of the black hole. Should always be between :math:`-1` and :math:`1`. :math:`a > 0` if the accretion disk orbits in the same direction as the hole rotates; :math:`a < 0` if it orbits in the opposite direction.
 
     See also:
      :cite:t:`Page_1974` Equation 15l
+
+    See also:
+        :meth:`calc_innermost_orbit` for the calculation of the innermost orbit of Kerr black holes.
     """
     a_ = a/bh_mass
     return 1 + (1-a_**2)**(1/3)*((1+a_)**(1/3) + (1 - a_)**(1/3))
 
+
 def calc_Z2(bh_mass, a):
-    """Calculate :math:`Z2`
+    r"""Calculate :math:`Z2` for Kerr black holes.
+
+    The variable :math:`Z2` is used to calculate the innermost orbit for Kerr black holes.
+
+    .. math::
+
+       Z_2 \equiv \sqrt{3a_*^2+Z_1^2}
+
+    Args:
+        bh_mass (float): Mass of the black hole.
+        a (float): Specific angular momentum of the black hole. Should always be between :math:`-1` and :math:`1`. :math:`a > 0` if the accretion disk orbits in the same direction as the hole rotates; :math:`a < 0` if it orbits in the opposite direction.
 
     See also:
         :cite:t:`Page_1974` Equation 15m
+
+    See also:
+        :meth:`calc_innermost_orbit` for the calculation of the innermost orbit of Kerr black holes.
     """
     Z1 = calc_Z1(bh_mass, a)
     a_ = a/bh_mass
     return np.sqrt(3*a_**2 + Z1**2)
 
 
-def calc_innermost_orbit(bh_mass, a):
-    """Calculcate the innermost orbit :math:`r_{ms}` for a Kerr black hole.
+def calc_innermost_stable_orbit(bh_mass, a):
+    r"""Calculcate the innermost stable orbit :math:`r_{ms}` for a Kerr black hole.
 
-    A larger angular momentum :math:`a` will yield innermost orbits closer to the black hole.
+    A larger specific angular momentum :math:`a` will yield innermost orbits closer to the black hole.
+
+    .. math::
+
+       \begin{align*}
+       r_{ms} &= Mx_0^2 \\
+       x_0^2 &= 3 + Z_2 - sgn(a_*)\sqrt{(3-Z_1)(3+Z_1+2Z_2) }
+       \end{align*}
+
+    Args:
+        bh_mass (float): Mass of the black hole.
+        a (float): Specific angular momentum of the black hole. Should always be between :math:`-1` and :math:`1`. :math:`a > 0` if the accretion disk orbits in the same direction as the hole rotates; :math:`a < 0` if it orbits in the opposite direction.
 
     See also:
         :cite:t:`Page_1974`
+
+    See also:
+        :meth:`calc_Z1` and :meth:`calc_Z2`.
     """
     Z1 = calc_Z1(bh_mass, a)
     Z2 = calc_Z2(bh_mass, a)
@@ -502,18 +545,25 @@ def calc_innermost_orbit(bh_mass, a):
         3 + Z2 - np.sign(a_)*((3-Z1)*(3+Z1+2*Z2))**.5
     )
 
+
 def calc_x0(bh_mass, a):
     r"""Calculcate :math:`x_0` for Kerr black holes.
 
     .. math::
 
-       x_0 = \sqrt{\frac{r_{ms}}/{M}}
+       x_0 = \sqrt{\frac{r_{ms}}{M}}
 
+    Args:
+        bh_mass (float): Mass of the black hole.
+        a (float): Specific angular momentum of the black hole. Should always be between :math:`-1` and :math:`1`. :math:`a > 0` if the accretion disk orbits in the same direction as the hole rotates; :math:`a < 0` if it orbits in the opposite direction.
     
     See also:
         :cite:t:`Page_1974` Equation 15k
+
+    See also:
+        meth:`calc_innermost_orbit` for the calculation of :math:`r_{ms}`
     """
-    rms = calc_innermost_orbit(bh_mass, a)
+    rms = calc_innermost_stable_orbit(bh_mass, a)
     return np.sqrt(rms/bh_mass)
 
 
@@ -526,7 +576,7 @@ def calc_f_kerr(bh_mass, a, r):
 
        F_s(r) = \frac{\dot{M}_0}{4\pi}e^{-(\nu + \psi + \mu)}f
 
-    Here, :math:`nu`, :math:`psi` and :math:`mu` are metric coefficients (functions of r), including Kerr metric. :math:`\dot{M}_0` is the radius-independent, time-averaged rate at which mass flows inward. Defining the innermost stable orbit as :math:`r_{ms}`, :math:`x=\sqrt{r/M}`, :math:`x_0=\sqrt{r_{ms}/M}` and :math:`a^*=a/M`, the :math:`f`-function is defined as:
+    Here, :math:`\nu`, :math:`\psi` and :math:`\mu` are metric coefficients (functions of :math:`r`) of the Kerr metric. :math:`\dot{M}_0` is the radius-independent, time-averaged rate at which mass flows inward. Defining the innermost stable orbit as :math:`r_{ms}`, :math:`x=\sqrt{r/M}=\sqrt{r^*}`, :math:`x_0=\sqrt{r_{ms}/M}` and :math:`a^*=a/M`, the :math:`f`-function is defined as:
 
     .. math::
 
@@ -537,6 +587,43 @@ def calc_f_kerr(bh_mass, a, r):
          &- \frac{3(x_3 - a^*)^2}{x_3(x_3-x_1)(x_3-x_2)}\ln\left(\frac{x-x_3}{x_0-x_3}\right) \Bigg]
         \end{align*}
 
+    , where
+
+    .. math::
+
+       \begin{align*}
+       x_1 &= 2\cos(\frac{1}{3}\cos^{-1}(a_*) - \frac{\pi}{3}) \\
+       x_2 &= 2\cos(\frac{1}{3}\cos^{-1}(a_*) + \frac{\pi}{3}) \\
+       x_3 &= -2\cos(\frac{1}{3}\cos^{-1}(a_*)) \\
+       \end{align*}
+
+    For a Swarzschild black hole, :math:`a=0` and these simplify to:
+
+    .. math::
+
+       \begin{align*}
+       x_1 &= \sqrt{3} \\
+       x_2 &= 0 \\
+       x_3 &= - \sqrt{3} \\
+       f &= \frac{3}{2M}\frac{1}{{r^{*}}^{1.5}(r^*-3)}\left[x - x_0 + \frac{\sqrt{3}}{2}\ln\left( \frac{(\sqrt{6} - \sqrt{3})(\sqrt{r^*}+\sqrt{3})}{(\sqrt{6} + \sqrt{3})(\sqrt{r^*} - \sqrt{3})}  \right)  \right]
+       \end{align*}
+
+    Args:
+        bh_mass (float): Mass of the black hole.
+        a (float): Specific angular momentum of the black hole. Should always be between :math:`-1` and :math:`1`. :math:`a > 0` if the accretion disk orbits in the same direction as the hole rotates; :math:`a < 0` if it orbits in the opposite direction.
+        r (float): Radius of the orbit
+
+    Attention:
+        :cite:t:`Luminet_1979` has a mistake in Equation 15. The factor in fromt of the :math:`log` should be :math:`\sqrt{3}/2` instead of :math:`\sqrt{3}/3`. This can be verified by solving :cite:t:`Page_1974` Equation 15n.  The resulting images of the paper are correct though.
+
+    See also:
+        :cite:t:`Page_1974` for more information.
+
+    See also:
+        :meth:`calc_flux_intrinsic_kerr` for the calculation of the intrinsic flux.
+
+    See also:
+       :meth:`calc_innermost_stable_orbit` for the calculation of :math:`r_{ms}`
     """
     a_ = a/bh_mass
     x = np.sqrt(r/bh_mass)
@@ -555,7 +642,35 @@ def calc_f_kerr(bh_mass, a, r):
 
 
 def calc_flux_intrinsic_kerr(bh_mass, a, r, acc):
-    """Calculate the intrinsic flux of the accretion disk of a Kerr black hole, in function of the accretion rate, specific angular momentum, and radius of emission."""
+    r"""Calculate the intrinsic flux of the accretion disk of a Kerr black hole, in function of the accretion rate, specific angular momentum, and radius of emission.
+    
+    The intrinsic flux is not redshift-corrected. Observed photons will have a flux that deviates from this by a factor of :math:`1/(1+z)^4`
+
+    The intrinsic flux in function of the radius is defined as:
+
+    .. math::
+
+       F_s(r) &= \frac{\dot{M_0}}{4\pi}e^{-(\nu+\psi+\mu)}f \\
+
+    where
+
+    .. math::
+
+       \begin{align*}
+       e^{\nu+\psi+\mu} &= r \\
+       f &= -\Omega_{,r}(E^{\dagger}-\Omega L^\dagger)^{-2}\int_{r_{ms}}^r(E^\dagger \ - \Omega L^\dagger)L^\dagger_{,r}dr
+       \end{align*}
+
+    Args:
+        bh_mass (float): Mass of the black hole.
+        a (float): Specific angular momentum of the black hole. Should always be between :math:`-1` and :math:`1`. :math:`a > 0` if the accretion disk orbits in the same direction as the hole rotates; :math:`a < 0` if it orbits in the opposite direction.
+        r (float): Radius of the orbit.
+        acc (float): (initial) accretion rate of the black hole :math:`\dot{M}_0`
+
+    See also:
+        :meth:`calc_f_kerr` for an algebraic expression of the :math:`f` function.
+
+    """
     f = calc_f_kerr(bh_mass=bh_mass, a=a, r=r)
     exp_nupsimu = r
     return acc * f / exp_nupsimu / 4 / np.pi
@@ -564,14 +679,14 @@ def calc_flux_intrinsic_kerr(bh_mass, a, r, acc):
 def calc_flux_intrinsic_swarzschild(bh_mass, r, acc):
     r"""Calculate the intrinsic flux of a photon.
     
-    The intrinsic flux is not redshift-corrected. Observed photons will have a flux
-    that deviates from this.
+    The intrinsic flux is not redshift-corrected. Observed photons will have a flux that deviates from this by a factor of :math:`1/(1+z)^4`
 
     .. math::
 
-        F_s = \frac{3 M \dot{M}}{8 \pi (r^* - 3) r^{5/2}} \left( \sqrt{r^*} - \sqrt{6} + \frac{1}{\sqrt{3}} \log \left( \frac{(\sqrt{r^*} + \sqrt{3})(\sqrt{6}-\sqrt{3})}{(\sqrt{6} + \sqrt{3})(\sqrt{r^*}-\sqrt{3})} \right) \right)
+       F_s = \frac{3 M \dot{M}}{8 \pi (r^* - 3) {r^*}^{5/2}} \left( \sqrt{r^*} - \sqrt{6} + \frac{\sqrt{3}}{2} \log \left( \frac{(\sqrt{r^*} + \sqrt{3})(\sqrt{6}-\sqrt{3})}{(\sqrt{6} + \sqrt{3})(\sqrt{r^*}-\sqrt{3})} \right) \right)
 
     where :math:`r^*=r/M`
+
     Args:
         r (float): radius on the accretion disk (BH frame)
         acc (float): accretion rate
@@ -581,7 +696,7 @@ def calc_flux_intrinsic_swarzschild(bh_mass, r, acc):
         float: Intrinsic flux of the photon :math:`F_s`
 
     Attention:
-        :cite:t:`Luminet_1979` has a mistake in Equation 15. The factor in fromt of the :math:`log` should be :math:`\sqrt{3}/2` instead of :math:`sqrt{3}/3`. This can be verified (tediously) by solving :cite:t:`Page_1974` Equation 15n.
+        :cite:t:`Luminet_1979` has a mistake in Equation 15. The factor in fromt of the :math:`log` should be :math:`\sqrt{3}/2` instead of :math:`\sqrt{3}/3`. This can be verified by solving :cite:t:`Page_1974` Equation 15n. The resulting images of the paper are correct though.
     """
     r_ = r / bh_mass
     log_arg = (np.sqrt(r_) + np.sqrt(3)) * (np.sqrt(6) - np.sqrt(3)) / ((np.sqrt(r_) - np.sqrt(3)) * (np.sqrt(6) + np.sqrt(3)))
@@ -621,7 +736,7 @@ def calc_redshift_factor(radius, angle, incl, bh_mass, b):
 
     Attention:
         :cite:t:`Luminet_1979` does not have the correct equation for the redshift factor.
-        The correct formula is given above.
+        The correct formula is given above. The resulting images of the paper are correct though.
     """
     # gff = (radius * np.sin(incl) * np.sin(angle)) ** 2
     # gtt = - (1 - (2. * M) / radius)