phys_552.bessel
===============

.. py:module:: phys_552.bessel

.. autoapi-nested-parse::

   Some utilities for computing properties of the Bessel functions for the DVR
   basis.



Functions
---------

.. autoapisummary::

   phys_552.bessel.sinc
   phys_552.bessel.J
   phys_552.bessel.j_root
   phys_552.bessel.J_sqrt_pole


Module Contents
---------------

.. py:function:: sinc(x, d=0)

   Return the `d`'th derivative of `sinc(x) = sin(x)/x`.

   :param x: Argument
   :type x: {float, array}
   :param d: Order of derivative.
   :type d: int, optional

   .. rubric:: Examples

   >>> N = 3
   >>> np.allclose(sinc(2.0), np.sin(2.0)/2.0)
   True
   >>> print("skip ...");from mmfutils.math.differentiate import differentiate
   skip ...
   >>> x = np.array([-1,0,0.5,1])
   >>> np.allclose(differentiate(lambda x:sinc(x), x),
   ...             sinc(x, d=1))
   True


.. py:function:: J(nu, d=0)

   Return the `d`'th derivative of the bessel functions
   :math:`J_{\nu}(z)`.


   :param nu: Order.
   :type nu: float
   :param d: Compute the `d`'th derivative.
   :type d: int

   .. rubric:: Examples

   >>> J0 = J(0.5); J1 = J(1.5); J2 = J(2.5);
   >>> z = 2.5; nu = 1.5
   >>> abs(J0(z) + J2(z) - 2*nu/z*J1(z)) < _EPS
   True

   .. todo:: Fix tolerances so that these are computed to machine precision.


.. py:function:: j_root(nu, N, rel_tol=2 * _EPS)

   Return the first N positive roots of the bessel function
   `J_nu(x)`.

   :param nu: Order of bessel function
   :type nu: float
   :param N: Number of roots
   :type N: int
   :param rel_tol: Desired relative tolerance for roots.
   :type rel_tol: float, optional

   .. admonition:: Notes

      The general method is to first estimate the roots with a
      bisection/secant method, and then polish them using Newton's
      method.
      
      We start by estimating the lower bound for the first for
      non-negative :math:`\nu`
      
      .. math::
         \nu + \nu^{1/3} < j_{\nu,1} < \nu + 1.85575 \nu^{1/3} + \pi
      
      Then, using the fact that the roots are spaced by at least
      :math:`\pi`, we step through the sign changes to bracket all of
      the desired roots.
      
      lowest root using the following
      heuristics
      
      .. math::
         j_{\nu,s} \approx \begin{cases}
            2\sqrt{\nu + 1}
               & -1 < \nu < -0.8\\
            \left(\frac{\nu}{2} + \frac{3}{4}\right)\pi
               & -0.8 < \nu < 2.5\\
            \nu + 1.85575 \nu^{1/3}
               & 2.5 \leq \nu
         \end{cases}

   .. rubric:: Examples

   >>> nu = 2.5
   >>> j_ = j_root(nu, 2000)
   >>> J_ = J(nu)(j_)

   These are roots!

   >>> np.max(abs(J_/j_)) < _EPS
   True

   They are also distinct

   >>> pi < min(np.diff(j_))
   True

   And the spacing is decreasing, meaning we have not skipped any.

   >>> np.max(np.diff(np.diff(j_))) < 0
   True


.. py:function:: J_sqrt_pole(nu, zn, d=0)

   Return a function that computes the `d`'th derivative of
   `sqrt(z)*J(nu,z)/(z - zn)` where `zn` is a root: `J(nu, zn) = 0`.
   This combination appears in the Bessel function DVR basis.

   :param nu: Order
   :type nu: float
   :param zn: Root of `J(nu, z)`
   :type zn: float
   :param d: Order of derivative to take.
   :type d: int

   .. admonition:: Notes

      .. math::
         \frac{\sqrt{z}J_{\nu}(z)}{z - z_{n}}
      
      As :math:`z` approaches :math:`z_n`, this has the form of `0/0`,
      so one can apply a form of l'Hopital's rule to reduce the
      round-off error.  The specified form of the function has been
      chosen for special properties of the Bessel functions.  Express
      the function as
      
      .. math::
         F(z) &= \frac{f(z)}{z - z_n}  = \frac{\sqrt{z}J_{\nu}(z)}{z - z_n}\\
         F'(z) &= \frac{f'(z)}{z - z_n} - \frac{f(z)}{(z - z_n)^2}

   Let :math:`\delta = z - z_n`.  Close to the singular point we use
   the Taylor series:

   .. math::
      \sum_{m=0}^{\infty}\frac{a_m\delta^{m}}{m!}

   .. math::
      F(z) &= f'(z_n)
           + \sum_{m=3}^{\infty}\frac{f^{(m)}(z_n)\delta^{m-1}}{m!}
           = \sum_{m=0}^{\infty}\frac{f^{(m+1)}(z_n)\delta^{m}}{(m+1)m!}\\
       a_m &= \frac{f^{(m+1)}}{m+1}\\
      F'(z) &= \sum_{m=3}^{\infty}\frac{(m-1)f^{(m)}(z_n)\delta^{m-2}}{m!}
            = \sum_{m=1}^{\infty}\frac{f^{(m+2)}(z_n)\delta^{m}}{(m+2)m!}\\
        a_m &= \frac{f^{(m+2)}}{m+2}\\

   The first few derivatives are presented here:

   .. math::
      f(z) &= \sqrt{z}J_{\nu}(z)\\
      f'(z) &= \frac{J_{\nu}(z)}{2\sqrt{z}} + \sqrt{z}J'_{\nu}(z)
             = \frac{f(z)}{2z} + \sqrt{z}J'_{\nu}(z)\\
      f''(z) &= \sqrt{z}J_{\nu}(z)\left(
                    \frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)
              = f(z)\left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)\\
      f'''(z) &= f'(z)\left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)
              - 2f(z)\frac{\nu^2 - \tfrac{1}{4}}{z^3}\\
      f^{(4)}(z) &=
      f(z)\left[
         \left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)^2
         + 6\frac{\nu^2 - \tfrac{1}{4}}{z^4}\right]
      - 4f'(z)\frac{\nu^2 - \tfrac{1}{4}}{z^3}
   .. Checked with Maple.

   Evaluated at the root :math:`z=z_n` these become:

   .. math::
      f(z_{n}) &= 0\\
      f'(z_{n}) &= \sqrt{z_{n}}J'_{\nu}(z_{n})\\
      f''(z_{n}) &= 0\\
      f'''(z_{n}) &= f'(z_{n})\left(
                       \frac{\nu^2 - \tfrac{1}{4}}{z_{n}^2} - 1\right)\\
      f^{(4)}(z_{n}) &= - 4f'(z_{n})\frac{\nu^2 - \tfrac{1}{4}}{z_{n}^3}

   with both the function and the second derivative vanishing.

   To determine where to use this formula, we match the estimate
   roundoff error with the truncation error.  The Bessel functions
   are of order unity and are typically calculated to an absolute
   accuracy of :math:`\epsilon`.  The round-off error in the
   numerator is :math:`\epsilon f(z)` and :math:`\epsilon \sqrt{2}
   z_n` in the denominator.  The roundoff errors in the denominator
   dominate both cases:

   .. math::
      \delta F(z) &\sim \epsilon \frac{\sqrt{2}z_n F(z)}{\delta}
                   \sim \frac{\sqrt{2}\epsilon z_n f(z)}{\delta^2}
                   \sim \frac{\sqrt{2}\epsilon z_n f'(z_n)}{\delta}\\
      \delta F'(z) &\sim \frac{2\epsilon z_n f(z)}{\delta^3}
                    \sim \frac{2\epsilon z_n f'(z_n)}{\delta^2}

   To choose the appropriate transition point, we equate half of this
   with the truncation error to transition points:

   .. math::
      \delta_c &\sim \left(
         \frac{72\epsilon z_n f'(z_n)}{\sqrt{2}f^{(4)}(z_n)}
         \right)^{1/4}
      \sim \left(\frac{72\epsilon z_n}{\sqrt{2}} \right)^{1/4}\\
      \delta_c' &\sim \left(120\epsilon z_n\right)^{1/5}

   the fact that :\math:`f(z)` behaves
   asymptotically as a :math:`\sqrt{2/\pi}\cos(z + \phi)` and so all
   derivatives have essentially the same magnitude.

   .. rubric:: Examples

   >>> nu = 5.5
   >>> zn = j_root(nu,21)[-1]
   >>> abs(zn - 73.62361318251753391646) < 1e-16
   True
   >>> float(J_sqrt_pole(nu,zn)(zn))   # doctest: +ELLIPSIS
   -0.796778576780013...

   -0.796778576780013129760

   You can also use a vector of `zn`, but only if it is commensurate
   with the argument:

   >>> zn = j_root(nu,21)
   >>> float(J_sqrt_pole(nu, zn)(zn[-1])[20]) # doctest: +ELLIPSIS
   -0.796778576780013...