.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/fftlogtest.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_fftlogtest.py: FFTLog-Test =========== This example is a translation of `fftlogtest.f` from the Fortran package `FFTLog`, which was presented in Appendix B of [Hami00]_ and published at https://jila.colorado.edu/~ajsh/FFTLog. It serves as an example for the python package `pyfftlog` (which is a Python version of `FFTLog`), in the same manner as the original file `fftlogtest.f` serves as an example for Fortran package `FFTLog`. **What follows is the original documentation from the file `fftlogtest.f`:** This is fftlogtest.f -------------------- This is a simple test program to illustrate how `FFTLog` works. The test transform is: .. math:: :label: hamtest1 \int^\infty_0 r^{\mu+1} \exp\left(-\frac{r^2}{2} \right)\ J_\mu(k, r)\ k\ {\rm d}r = k^{\mu+1} \exp\left(-\frac{k^2}{2} \right) Disclaimer ''''''''''' `FFTLog` does NOT claim to provide the most accurate possible solution of the continuous transform (which is the stated aim of some other codes). Rather, `FFTLog` claims to solve the exact discrete transform of a logarithmically-spaced periodic sequence. If the periodic interval is wide enough, the resolution high enough, and the function well enough behaved outside the periodic interval, then `FFTLog` may yield a satisfactory approximation to the continuous transform. Observe: 1. How the result improves as the periodic interval is enlarged. With the normal FFT, one is not used to ranges orders of magnitude wide, but this is how `FFTLog` prefers it. 2. How the result improves as the resolution is increased. Because the function is rather smooth, modest resolution actually works quite well here. 3. That the central part of the transform is more reliable than the outer parts. Experience suggests that a good general strategy is to double the periodic interval over which the input function is defined, and then to discard the outer half of the transform. 4. That the best bias exponent seems to be :math:`q = 0`. 5. That for the critical index :math:`\mu = -1`, the result seems to be offset by a constant from the 'correct' answer. 6. That the result grows progressively worse as mu decreases below -1. The analytic integral above fails for :math:`\mu \le -1`, but `FFTLog` still returns answers. Namely, `FFTLog` returns the analytic continuation of the discrete transform. Because of ambiguity in the path of integration around poles, this analytic continuation is liable to differ, for :math:`\mu \le -1`, by a constant from the 'correct' continuation given by the above equation. `FFTLog` begins to have serious difficulties with aliasing as :math:`\mu` decreases below :math:`-1`, because then :math:`r^{\mu+1} \exp(-r^2/2)` is far from resembling a periodic function. You might have thought that it would help to introduce a bias exponent :math:`q = \mu`, or perhaps :math:`q = \mu+1`, or more, to make the function :math:`a(r) = A(r) r^{-q}` input to `fhtq` more nearly periodic. In practice a nonzero :math:`q` makes things worse. A symmetry argument lends support to the notion that the best exponent here should be :math:`q = 0,` as empirically appears to be true. The symmetry argument is that the function :math:`r^{\mu+1} \exp(-r^2/2)` happens to be the same as its transform :math:`k^{\mu+1} \exp(-k^2/2)`. If the best bias exponent were q in the forward transform, then the best exponent would be :math:`-q` that in the backward transform; but the two transforms happen to be the same in this case, suggesting :math:`q = -q`, hence :math:`q = 0`. This example illustrates that you cannot always tell just by looking at a function what the best bias exponent :math:`q` should be. You also have to look at its transform. The best exponent :math:`q` is, in a sense, the one that makes both the function and its transform look most nearly periodic. .. GENERATED FROM PYTHON SOURCE LINES 82-87 .. code-block:: Python import pyfftlog import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 88-92 Define the parameters you wish to use ------------------------------------- The presets are the *Reasonable choices of parameters* from `fftlogtest.f`. .. GENERATED FROM PYTHON SOURCE LINES 92-117 .. code-block:: Python # Range of periodic interval logrmin = -4 logrmax = 4 # Number of points (Max 4096) n = 64 # Order mu of Bessel function mu = 0 # Bias exponent: q = 0 is unbiased q = 0 # Sensible approximate choice of k_c r_c kr = 1 # Tell fhti to change kr to low-ringing value # WARNING: kropt = 3 will fail, as interaction is not supported kropt = 1 # Forward transform (changed from dir to tdir, as dir is a python fct) tdir = 1 .. GENERATED FROM PYTHON SOURCE LINES 118-120 Computation related to the logarithmic spacing ---------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 120-134 .. code-block:: Python # Central point log10(r_c) of periodic interval logrc = (logrmin + logrmax)/2 print(f"Central point of periodic interval at log10(r_c) = {logrc}") # Central index (1/2 integral if n is even) nc = (n + 1)/2.0 # Log-spacing of points dlogr = (logrmax - logrmin)/n dlnr = dlogr*np.log(10.0) .. rst-class:: sphx-glr-script-out .. code-block:: none Central point of periodic interval at log10(r_c) = 0.0 .. GENERATED FROM PYTHON SOURCE LINES 135-137 Compute input function: :math:`r^{\mu+1}\exp\left(-\frac{r^2}{2}\right)` ------------------------------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 137-142 .. code-block:: Python r = 10**(logrc + (np.arange(1, n+1) - nc)*dlogr) ar = r**(mu + 1)*np.exp(-r**2/2.0) .. GENERATED FROM PYTHON SOURCE LINES 143-145 Initialize FFTLog transform - note fhti resets `kr` --------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 145-149 .. code-block:: Python kr, xsave = pyfftlog.fhti(n, mu, dlnr, q, kr, kropt) print(f"pyfftlog.fhti: new kr = {kr}") .. rst-class:: sphx-glr-script-out .. code-block:: none pyfftlog.fhti: new kr = 0.9535389675791917 .. GENERATED FROM PYTHON SOURCE LINES 150-152 Call `pyfftlog.fht` (or `pyfftlog.fhtl`) ---------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 152-163 .. code-block:: Python logkc = np.log10(kr) - logrc print(f"Central point in k-space at log10(k_c) = {logkc}") # rk = r_c/k_c rk = 10**(logrc - logkc) # Transform # ak = pyfftlog.fftl(ar.copy(), xsave, rk, tdir) ak = pyfftlog.fht(ar.copy(), xsave, tdir) .. rst-class:: sphx-glr-script-out .. code-block:: none Central point in k-space at log10(k_c) = -0.020661554260541743 .. GENERATED FROM PYTHON SOURCE LINES 164-166 Compute Output function: :math:`k^{\mu+1}\exp\left(-\frac{k^2}{2}\right)` ------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 166-170 .. code-block:: Python k = 10**(logkc + (np.arange(1, n+1) - nc)*dlogr) theo = k**(mu + 1)*np.exp(-k**2/2.0) .. GENERATED FROM PYTHON SOURCE LINES 171-173 Plot result ----------- .. GENERATED FROM PYTHON SOURCE LINES 173-213 .. code-block:: Python plt.figure() # Input ax1 = plt.subplot(121) plt.title(r'$r^{\mu+1}\ \exp(-r^2/2)$') plt.xlabel('r') plt.loglog(r, ar, 'k', lw=2) plt.grid(axis='y', c='0.9') # Transformed result ax2 = plt.subplot(122, sharey=ax1) plt.title(r'$k^{\mu+1} \exp(-k^2/2)$') plt.xlabel('k') plt.loglog(k, theo, 'k', lw=2, label='Theoretical') plt.loglog(k, ak, 'r--', lw=2, label='FFTLog') plt.legend() plt.ylim([1e-8, 1e1]) ax2.yaxis.tick_right() ax2.yaxis.set_label_position("right") plt.grid(axis='y', c='0.9') # Switch off spines ax1.spines['top'].set_visible(False) ax1.spines['right'].set_visible(False) ax2.spines['top'].set_visible(False) ax2.spines['left'].set_visible(False) plt.tight_layout(rect=[0, 0, 1, .9]) # Main title plt.suptitle(r"$\int_0^\infty r^{\mu+1}\ \exp(-r^2/2)\ J_\mu(k,r)\ " + r"k\ {\rm d}r = k^{\mu+1} \exp(-k^2/2)$", y=0.98) plt.show() .. image-sg:: /examples/images/sphx_glr_fftlogtest_001.png :alt: $\int_0^\infty r^{\mu+1}\ \exp(-r^2/2)\ J_\mu(k,r)\ k\ {\rm d}r = k^{\mu+1} \exp(-k^2/2)$, $r^{\mu+1}\ \exp(-r^2/2)$, $k^{\mu+1} \exp(-k^2/2)$ :srcset: /examples/images/sphx_glr_fftlogtest_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 214-216 Print values ------------ .. GENERATED FROM PYTHON SOURCE LINES 216-221 .. code-block:: Python print(' k a(k) k^(mu+1) exp(-k^2/2)') print('----------------------------------------------------------------') for i in range(n): print(f"{k[i]:18.6e} {ak[i]:18.6e} {theo[i]:18.6e}") .. rst-class:: sphx-glr-script-out .. code-block:: none k a(k) k^(mu+1) exp(-k^2/2) ---------------------------------------------------------------- 1.101130e-04 6.332603e-05 1.101130e-04 1.468380e-04 9.168618e-05 1.468380e-04 1.958116e-04 1.374282e-04 1.958116e-04 2.611190e-04 2.131954e-04 2.611190e-04 3.482078e-04 3.318802e-04 3.482077e-04 4.643425e-04 4.923984e-04 4.643425e-04 6.192107e-04 6.460278e-04 6.192106e-04 8.257307e-04 7.968931e-04 8.257304e-04 1.101130e-03 1.113736e-03 1.101129e-03 1.468380e-03 1.464233e-03 1.468378e-03 1.958116e-03 1.959475e-03 1.958112e-03 2.611190e-03 2.610678e-03 2.611181e-03 3.482078e-03 3.482260e-03 3.482056e-03 4.643425e-03 4.643299e-03 4.643375e-03 6.192107e-03 6.191999e-03 6.191988e-03 8.257307e-03 8.257056e-03 8.257026e-03 1.101130e-02 1.101057e-02 1.101063e-02 1.468380e-02 1.468230e-02 1.468222e-02 1.958116e-02 1.957729e-02 1.957741e-02 2.611190e-02 2.610314e-02 2.610300e-02 3.482078e-02 3.479950e-02 3.479967e-02 4.643425e-02 4.638444e-02 4.638422e-02 6.192107e-02 6.180220e-02 6.180247e-02 8.257307e-02 8.229239e-02 8.229205e-02 1.101130e-01 1.094470e-01 1.094474e-01 1.468380e-01 1.452640e-01 1.452635e-01 1.958116e-01 1.920928e-01 1.920934e-01 2.611190e-01 2.523680e-01 2.523671e-01 3.482078e-01 3.277241e-01 3.277250e-01 4.643425e-01 4.168889e-01 4.168871e-01 6.192107e-01 5.111853e-01 5.111866e-01 8.257307e-01 5.871956e-01 5.871927e-01 1.101130e+00 6.005500e-01 6.005516e-01 1.468380e+00 4.996049e-01 4.996187e-01 1.958116e+00 2.879340e-01 2.879045e-01 2.611190e+00 8.632888e-02 8.634968e-02 3.482078e+00 8.102022e-03 8.108819e-03 4.643425e+00 1.180344e-04 9.656964e-05 6.192107e+00 -1.553139e-05 2.923736e-08 8.257307e+00 7.225353e-06 1.291402e-14 1.101130e+01 -2.588950e-06 5.164519e-26 1.468380e+01 7.719794e-07 2.222595e-46 1.958116e+01 1.586977e-07 1.078537e-82 2.611190e+01 -1.874092e-07 2.285953e-147 3.482078e+01 5.576689e-07 1.793712e-262 4.643425e+01 -1.317041e-07 0.000000e+00 6.192107e+01 6.415736e-07 0.000000e+00 8.257307e+01 1.351283e-07 0.000000e+00 1.101130e+02 7.997181e-07 0.000000e+00 1.468380e+02 5.394094e-07 0.000000e+00 1.958116e+02 1.165867e-06 0.000000e+00 2.611190e+02 1.176786e-06 0.000000e+00 3.482078e+02 1.889416e-06 0.000000e+00 4.643425e+02 2.248731e-06 0.000000e+00 6.192107e+02 3.228937e-06 0.000000e+00 8.257307e+02 4.113223e-06 0.000000e+00 1.101130e+03 5.651921e-06 0.000000e+00 1.468380e+03 7.408687e-06 0.000000e+00 1.958116e+03 1.001142e-05 0.000000e+00 2.611190e+03 1.330606e-05 0.000000e+00 3.482078e+03 1.792186e-05 0.000000e+00 4.643425e+03 2.410633e-05 0.000000e+00 6.192107e+03 3.277422e-05 0.000000e+00 8.257307e+03 4.510046e-05 0.000000e+00 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.541 seconds) .. _sphx_glr_download_examples_fftlogtest.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: fftlogtest.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: fftlogtest.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: fftlogtest.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_