pyextal.jaxMultipole.radial
===========================

.. py:module:: pyextal.jaxMultipole.radial

.. autoapi-nested-parse::

   Fourier-Bessel transforms of Slater-type radial functions.

   Every radial integral the model needs is one of

   .. math::

       \langle j_l \rangle_n(k; \alpha) = \frac{\alpha^{n+3}}{(n+2)!}
           \int_0^\infty r^{n+2} e^{-\alpha r} j_l(kr)\, dr

   -- the transform of a unit-normalised Slater density ``alpha^(n+3)/(n+2)!
   r^n e^(-alpha r)``. The spherical core and valence factors are sums of
   ``<j_0>`` over pairs of basis functions (`jaxMultipole.bank.pair_table`); the
   deformation terms are a single ``<j_l>`` each.

   Stable closed form
   ------------------

   The integral is a Gauss hypergeometric function,

   .. math::

       \langle j_l \rangle_n = C_{n,l}\, u^l\, {}_2F_1\!\left(\tfrac{n+l+3}{2},
           \tfrac{n+l+4}{2}; l+\tfrac32; -u^2\right), \qquad
       u = k/\alpha, \quad C_{n,l} = \frac{(n+l+2)!}{(2l+1)!!\,(n+2)!},

   and a Pfaff transformation turns it into ``(1+u^2)^(-p)`` times a *terminating*
   polynomial in ``t = u^2/(1+u^2)``, which lies in [0, 1):

   * ``n - l`` even: ``p = (n+l+4)/2``, degree ``(n-l)/2``,
     ``c_q = ((l-n)/2)_q ((n+l+4)/2)_q / ((l+3/2)_q q!)``;
   * ``n - l`` odd: ``p = (n+l+3)/2``, degree ``(n-l+1)/2``,
     ``c_q = ((n+l+3)/2)_q ((l-n-1)/2)_q / ((l+3/2)_q q!)``.

   Both need ``n >= l - 1``. The coefficients are exact rationals, the polynomial
   is evaluated by Horner's rule, and nothing cancels at small ``k`` -- unlike the
   textbook route through ``m!/(alpha - ik)^(m+1)`` and the trigonometric expansion
   of ``j_l``, which for ``n = l = 4`` is already 41 times off at ``k = 0.05`` in
   float64 and useless in float32. It is exact at ``k = 0`` (1 for ``l = 0``, 0
   otherwise) and smooth in ``alpha``, so derivatives in kappa are clean. For
   example ``<j_0>_2 = (1 - u^2)/(1 + u^2)^4``.

   Units: ``k`` and ``alpha`` must share a length unit. The banks give exponents in
   bohr^-1, so the model passes ``k = 4 pi s * BOHR``.



Functions
---------

.. autoapisummary::

   pyextal.jaxMultipole.radial.fb_coefficients
   pyextal.jaxMultipole.radial.sto_fb
   pyextal.jaxMultipole.radial.k_bohr
   pyextal.jaxMultipole.radial.spherical_f
   pyextal.jaxMultipole.radial.deformation_fb


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

.. py:function:: fb_coefficients(n, l)

   Exact constants of :func:`sto_fb`.

   :param n: Power of r in the Slater density.
   :type n: int
   :param l: Order of the Bessel function.
   :type l: int

   :returns: ``(C, p, coeffs)`` with ``C`` the prefactor, ``p`` the integer
             power of ``1/(1+u^2)`` and ``coeffs`` the polynomial in ``t``, lowest
             order first, all as floats.
   :rtype: tuple

   :raises ValueError: If ``n < l - 1``, where the integral diverges at r = 0.


.. py:function:: sto_fb(k, alpha, n, l)

   ``<j_l>_n(k; alpha)`` for a unit-normalised Slater density.

   :param k: Scattering vector magnitude ``4 pi s``, any shape.
   :type k: array-like
   :param alpha: Slater exponent(s), broadcasting against ``k``.
   :type alpha: array-like
   :param n: Power of r.
   :type n: int
   :param l: Bessel order.
   :type l: int

   :returns: The transform, broadcast shape, dtype of ``k * alpha``.
   :rtype: jnp.ndarray


.. py:function:: k_bohr(s)

   ``4 pi s`` in bohr^-1, from ``s = sin(theta)/lambda`` in Angstrom^-1.


.. py:function:: spherical_f(k, table, kappa=1.0)

   Scattering factor of a spherical shell, one electron's worth.

   :param k: ``4 pi s`` in bohr^-1.
   :type k: array-like
   :param table: ``N -> (w, alpha)`` from `jaxMultipole.bank.pair_table`.
   :type table: dict
   :param kappa: Expansion parameter; the density
                 ``kappa^3 rho(kappa r)`` scatters as ``f(k/kappa)``.
   :type kappa: float | array

   :returns: ``f(k)``, shape of ``k``; zero for an empty table.
   :rtype: jnp.ndarray


.. py:function:: deformation_fb(k, l, n_l, zeta_l, kappa_p=1.0)

   ``<j_l>`` of a deformation radial function ``R_l ~ r^n_l e^(-kappa' zeta_l r)``.

   :param k: ``4 pi s`` in bohr^-1.
   :type k: array-like
   :param l: Multipole order.
   :type l: int
   :param n_l: Power of r; ``n_l >= l`` is required for a density regular at
               the nucleus (the transform itself only needs ``n_l >= l - 1``).
   :type n_l: int
   :param zeta_l: Exponent, bohr^-1.
   :type zeta_l: float
   :param kappa_p: kappa'.
   :type kappa_p: float | array

   :returns: Same shape as ``k``.
   :rtype: jnp.ndarray