9=eg.dZddlZddlmZddlmZmZddlm Z m Z gdZ ej dZ dd Zdd Zdd Zdd ZddZdS)zFast Hankel transforms using the FFTLog algorithm. The implementation closely follows the Fortran code of Hamilton (2000). added: 14/11/2020 Nicolas Tessore N)warn)rfftirfft)loggammapoch)fhtifht fhtoffsetc\tj|d}|dkr=|dz dz }tj|}|tj| ||z z|zz}t |||||}t ||} |dkr$| tj| ||z |z|zzz} | S)acCompute the fast Hankel transform. Computes the discrete Hankel transform of a logarithmically spaced periodic sequence using the FFTLog algorithm [1]_, [2]_. Parameters ---------- a : array_like (..., n) Real periodic input array, uniformly logarithmically spaced. For multidimensional input, the transform is performed over the last axis. dln : float Uniform logarithmic spacing of the input array. mu : float Order of the Hankel transform, any positive or negative real number. offset : float, optional Offset of the uniform logarithmic spacing of the output array. bias : float, optional Exponent of power law bias, any positive or negative real number. Returns ------- A : array_like (..., n) The transformed output array, which is real, periodic, uniformly logarithmically spaced, and of the same shape as the input array. See Also -------- ifht : The inverse of `fht`. fhtoffset : Return an optimal offset for `fht`. Notes ----- This function computes a discrete version of the Hankel transform .. math:: A(k) = \int_{0}^{\infty} \! a(r) \, J_\mu(kr) \, k \, dr \;, where :math:`J_\mu` is the Bessel function of order :math:`\mu`. The index :math:`\mu` may be any real number, positive or negative. The input array `a` is a periodic sequence of length :math:`n`, uniformly logarithmically spaced with spacing `dln`, .. math:: a_j = a(r_j) \;, \quad r_j = r_c \exp[(j-j_c) \, \mathtt{dln}] centred about the point :math:`r_c`. Note that the central index :math:`j_c = (n-1)/2` is half-integral if :math:`n` is even, so that :math:`r_c` falls between two input elements. Similarly, the output array `A` is a periodic sequence of length :math:`n`, also uniformly logarithmically spaced with spacing `dln` .. math:: A_j = A(k_j) \;, \quad k_j = k_c \exp[(j-j_c) \, \mathtt{dln}] centred about the point :math:`k_c`. The centre points :math:`r_c` and :math:`k_c` of the periodic intervals may be chosen arbitrarily, but it would be usual to choose the product :math:`k_c r_c = k_j r_{n-1-j} = k_{n-1-j} r_j` to be unity. This can be changed using the `offset` parameter, which controls the logarithmic offset :math:`\log(k_c) = \mathtt{offset} - \log(r_c)` of the output array. Choosing an optimal value for `offset` may reduce ringing of the discrete Hankel transform. If the `bias` parameter is nonzero, this function computes a discrete version of the biased Hankel transform .. math:: A(k) = \int_{0}^{\infty} \! a_q(r) \, (kr)^q \, J_\mu(kr) \, k \, dr where :math:`q` is the value of `bias`, and a power law bias :math:`a_q(r) = a(r) \, (kr)^{-q}` is applied to the input sequence. Biasing the transform can help approximate the continuous transform of :math:`a(r)` if there is a value :math:`q` such that :math:`a_q(r)` is close to a periodic sequence, in which case the resulting :math:`A(k)` will be close to the continuous transform. References ---------- .. [1] Talman J. D., 1978, J. Comp. Phys., 29, 35 .. [2] Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191) Examples -------- This example is the adapted version of ``fftlogtest.f`` which is provided in [2]_. It evaluates the integral .. math:: \int^\infty_0 r^{\mu+1} \exp(-r^2/2) J_\mu(k, r) k dr = k^{\mu+1} \exp(-k^2/2) . >>> import numpy as np >>> from scipy import fft >>> import matplotlib.pyplot as plt Parameters for the transform. >>> mu = 0.0 # Order mu of Bessel function >>> r = np.logspace(-7, 1, 128) # Input evaluation points >>> dln = np.log(r[1]/r[0]) # Step size >>> offset = fft.fhtoffset(dln, initial=-6*np.log(10), mu=mu) >>> k = np.exp(offset)/r[::-1] # Output evaluation points Define the analytical function. >>> def f(x, mu): ... """Analytical function: x^(mu+1) exp(-x^2/2).""" ... return x**(mu + 1)*np.exp(-x**2/2) Evaluate the function at ``r`` and compute the corresponding values at ``k`` using FFTLog. >>> a_r = f(r, mu) >>> fht = fft.fht(a_r, dln, mu=mu, offset=offset) For this example we can actually compute the analytical response (which in this case is the same as the input function) for comparison and compute the relative error. >>> a_k = f(k, mu) >>> rel_err = abs((fht-a_k)/a_k) Plot the result. >>> figargs = {'sharex': True, 'sharey': True, 'constrained_layout': True} >>> fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), **figargs) >>> ax1.set_title(r'$r^{\mu+1}\ \exp(-r^2/2)$') >>> ax1.loglog(r, a_r, 'k', lw=2) >>> ax1.set_xlabel('r') >>> ax2.set_title(r'$k^{\mu+1} \exp(-k^2/2)$') >>> ax2.loglog(k, a_k, 'k', lw=2, label='Analytical') >>> ax2.loglog(k, fht, 'C3--', lw=2, label='FFTLog') >>> ax2.set_xlabel('k') >>> ax2.legend(loc=3, framealpha=1) >>> ax2.set_ylim([1e-10, 1e1]) >>> ax2b = ax2.twinx() >>> ax2b.loglog(k, rel_err, 'C0', label='Rel. Error (-)') >>> ax2b.set_ylabel('Rel. Error (-)', color='C0') >>> ax2b.tick_params(axis='y', labelcolor='C0') >>> ax2b.legend(loc=4, framealpha=1) >>> ax2b.set_ylim([1e-9, 1e-3]) >>> plt.show() rrroffsetbiasnpshapearangeexpfhtcoeff_fhtq) adlnmurrnj_cjuAs 1lib/python3.11/site-packages/scipy/fft/_fftlog.pyr r sx  BA qyysAg IaLL ua#gs*++ + CF666A a A qyy RVTEAGS=612 3 33 Hc^tj|d}|dkr?|dz dz }tj|}|tj|||z |z|zzz}t |||||}t ||d} |dkr!| tj| ||z z|zz} | S)aCompute the inverse fast Hankel transform. Computes the discrete inverse Hankel transform of a logarithmically spaced periodic sequence. This is the inverse operation to `fht`. Parameters ---------- A : array_like (..., n) Real periodic input array, uniformly logarithmically spaced. For multidimensional input, the transform is performed over the last axis. dln : float Uniform logarithmic spacing of the input array. mu : float Order of the Hankel transform, any positive or negative real number. offset : float, optional Offset of the uniform logarithmic spacing of the output array. bias : float, optional Exponent of power law bias, any positive or negative real number. Returns ------- a : array_like (..., n) The transformed output array, which is real, periodic, uniformly logarithmically spaced, and of the same shape as the input array. See Also -------- fht : Definition of the fast Hankel transform. fhtoffset : Return an optimal offset for `ifht`. Notes ----- This function computes a discrete version of the Hankel transform .. math:: a(r) = \int_{0}^{\infty} \! A(k) \, J_\mu(kr) \, r \, dk \;, where :math:`J_\mu` is the Bessel function of order :math:`\mu`. The index :math:`\mu` may be any real number, positive or negative. See `fht` for further details. rrrrrT)inverser) r!rrrrrrrr rs r"r r s^  BA qyysAg IaLL ta#gs]V3455 5 CF666A aD!!!A qyy RVTE1s7OC' ( (( Hr#c||}}|dz|zdz }|dz|z dz }tjdtj|dzz||zz |dzdz} tj|dzdzt} tj|dzdzt} | | jdd<|| jdd<t| | || jdd<t| | | dt|z zz} | xj| jzc_| xjt|zz c_| xj| jz c_| xj| z c_tj | | d| jd<tj | dsd|zt|||z z| d<| S)z?Compute the coefficient array for a fast Hankel transform. rrr)dtypeN)outr) rlinspacepiemptycompleximagrealrLN_2risfiniter ) rrrrrlnkrqxpxmyr vs r"rrsd!D Q$q&!B Q$q&!B Aruad|QsU+QT!V44A Aaw'''A Aaw'''AAF111IAF111I QAAF111I QAD4KAFFafFFFFd1fFFFFafFFFFaKFFF1!AF2J ;qt  &!td2r"uoo%! Hr#c<||}}|dz|zdz }|dz|z dz }tjd|zz }t|d|zz} t|d|zz} t|z |z | j| jztjz z} || tj| z |zzS)aReturn optimal offset for a fast Hankel transform. Returns an offset close to `initial` that fulfils the low-ringing condition of [1]_ for the fast Hankel transform `fht` with logarithmic spacing `dln`, order `mu` and bias `bias`. Parameters ---------- dln : float Uniform logarithmic spacing of the transform. mu : float Order of the Hankel transform, any positive or negative real number. initial : float, optional Initial value for the offset. Returns the closest value that fulfils the low-ringing condition. bias : float, optional Exponent of power law bias, any positive or negative real number. Returns ------- offset : float Optimal offset of the uniform logarithmic spacing of the transform that fulfils a low-ringing condition. See Also -------- fht : Definition of the fast Hankel transform. References ---------- .. [1] Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191) rry?)rr*rr/r-round) rrinitialrr1r2r3r4r5zpzmargs r"r r 8sFt!D Q$q&!B Q$q&!B qu A "r!t)  B "r!t)  B $; rw0"%7 7C 3#&+ ++r#Fctj|d}tj|dr+|s)td|}d|d<n@|ddkr4|r2td|}tj|d<t |d}|s||z}n||z}t||d}|ddddf}|S)zUCompute the biased fast Hankel transform. This is the basic FFTLog routine. rrz.singular transform; consider changing the biasz6singular inverse transform; consider changing the bias)axis.N) rrisinfrcopyinfrconjr)rr r%rr!s r"rrfs  BA x!~~ g  =>>> FFHH! 1w EFFF FFHHv! QRA  Q QVVXX  aA #ttt) A Hr#)r r )F)__doc__numpyrwarningsr_basicrrspecialrr __all__logr/r r rr rr#r"rKs$$$$$$$$   rvayyp p p p fC C C C L% % % % P+,+,+,+,\       r#