{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Demonstration of Hankel transform identities\n\nBelow we demonstrate a range of known Hankel transform pairs\nfrom various sources.\n\nFirst we demonstrate the Gaussian function from Pissens [#Pissens]_ and its\ninverse transform.\n\nThen we check the \"generalised top-hat\" and \"generalised jinc\"\nfunctions from Guizar-Sicairos and Guitierrez-Vega [#Guizar]_.\n\nFinally, we look at the function $f(r) = \\frac{1}{r^2 + a^2}$,\nthe Hankel transform of which is $K_0(av)$, where $K_0$\nis the modified Bessel function of the second kind of order 0. [#Pissens]_\n\n\n.. [#Pissens] *\u201cChapter 9: The Hankel Transform.\u201d* Piessens, R.\n  in The Transforms and Applications Handbook: Second Edition.\n  Ed. Alexander D. Poularikas\n  Boca Raton: CRC Press LLC, 2000\n.. [#Guizar] *\"Computation of quasi-discrete Hankel transforms of the integer\n  order for propagating optical wave fields\"*\n  Manuel Guizar-Sicairos and Julio C. Guitierrez-Vega\n  J. Opt. Soc. Am. A **21** (1) 53-58 (2004)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport scipy.special as scipy_bessel\nimport matplotlib.pyplot as plt\n\nfrom pyhank import qdht, iqdht, HankelTransform"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First we try a Gaussian function, the Hankel transform of which should also be Gaussian.\n\nNote the definition in Guizar-Sicairos [#Guizar]_ varies from that used by\nPissens [#Pissens]_ by a factor of $2\\pi$ in\nboth scaling of the argument (so we use ``HankelTransform.kr`` rather than\n``HankelTransform.v``) and also scaling of the magnitude.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = 3\nradius = np.linspace(0, 3, 1024)\nf = np.exp(-a ** 2 * radius ** 2)\nkr, actual_ht = qdht(radius, f)\nexpected_ht = 2*np.pi*(1 / (2 * a**2)) * np.exp(-kr**2 / (4 * a**2))\nassert np.allclose(expected_ht, actual_ht)\n\nplt.figure()\nplt.subplot(2, 1, 1)\nplt.title('Gaussian function')\nplt.plot(radius, f)\nplt.xlabel('Radius /$r$')\nplt.subplot(2, 1, 2)\nplt.plot(kr, expected_ht, label='Analytical')\nplt.plot(kr, actual_ht, marker='x', linestyle='None', label='QDHT')\nplt.title('Hankel transform - also Gaussian')\nplt.xlabel('Frequency /$v$')\nplt.xlim([0, 50])\nplt.legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we repeat for the inverse transform\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kr = np.linspace(0, 50, 1024)\nht = 2*np.pi*(1 / (2 * a**2)) * np.exp(-kr**2 / (4 * a**2))\nr, actual_f = iqdht(kr, ht)\nexpected_f = np.exp(-a ** 2 * r ** 2)\nassert np.allclose(expected_f, actual_f)\nplt.figure()\nplt.subplot(2, 1, 1)\nplt.title('Hankel transform - Gaussian function')\nplt.plot(kr, ht)\nplt.xlabel('Radius /$r$')\nplt.subplot(2, 1, 2)\nplt.plot(radius, expected_f, label='Analytical')\nplt.plot(radius, actual_f, marker='x', linestyle='None', label='QDHT')\nplt.title('Original function after IQDHT - also Gaussian')\nplt.xlabel('Frequency /$v$')\nplt.xlim([0, 0.2])\nplt.legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next we define functions to calculate the generalised top-hat and jinc\nfunctions, as defined by Guizar-Sicairos and Guitierrez-Vega [#Guizar]_.\n\nNote that for $p=0$ these become a standard top-hat and\n$\\textrm{jinc}(r) = \\frac{J_1(r)}{r}$ functions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def generalised_top_hat(r: np.ndarray, a: float, p: int) -> np.ndarray:\n    top_hat = np.zeros_like(r)\n    top_hat[r <= a] = 1\n    return r ** p * top_hat\n\n\ndef generalised_jinc(v: np.ndarray, a: float, p: int):\n    val = np.zeros_like(v)\n    val[v != 0] = a ** (p + 1) * scipy_bessel.jv(p + 1, 2 * np.pi * a * v[v != 0]) / v[v != 0]\n    if p == -1:\n        val[v == 0] = np.inf\n    elif p == -2:\n        val[v == 0] = -np.pi\n    elif p == 0:\n        val[v == 0] = np.pi * a ** 2\n    else:\n        val[v == 0] = 0\n    return val"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For demonstration, we choose $a = 0.5$ and run the code for\norders 0, 1 and 4 plotting and checking the mean absolute error each time.\nFirst check that the Hankel transform of the generalised jinc is calculated\ncorrectly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radius = np.linspace(0, 30, 1024)\na = 0.5\nfor order in [0, 1, 4]:\n    f = generalised_jinc(radius, a, order)\n    kr, actual_ht = qdht(radius, f, order=order)\n    v = kr / (2*np.pi)\n    expected_ht = generalised_top_hat(v, a, order)\n\n    plt.figure()\n    plt.subplot(2, 1, 1)\n    plt.title(f'Generalised jinc function, order = {order}')\n    plt.plot(radius, f)\n    plt.xlabel('Radius /$r$')\n    plt.subplot(2, 1, 2)\n    plt.plot(v, expected_ht, label='Analytical')\n    plt.plot(v, actual_ht, marker='x', linestyle='None', label='QDHT')\n    plt.title(f'Hankel transform - generalised top-hat, order = {order}')\n    plt.xlabel('Frequency /$v$')\n    plt.xlim([0, 1.5])\n    plt.legend()\n    plt.tight_layout()\n\n    error = np.mean(np.abs(expected_ht-actual_ht))\n    assert error < 1e-3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we repeat but the other way round: the Hankel transform of the top-hat\nfunction should be the jinc function.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radius = np.linspace(0, 2, 1024)\nfor order in [0, 1, 4]:\n    transformer = HankelTransform(order=order, radial_grid=radius)\n    f = generalised_top_hat(transformer.r, a, order)\n    actual_ht = transformer.qdht(f)\n    expected_ht = generalised_jinc(transformer.v, a, order)\n\n    plt.figure()\n    plt.subplot(2, 1, 1)\n    plt.title(f'Generalised top-hat function, order = {order}')\n    plt.plot(radius, f)\n    plt.xlabel('Radius /$r$')\n    plt.subplot(2, 1, 2)\n    plt.plot(v, expected_ht, label='Analytical')\n    plt.plot(v, actual_ht, marker='x', linestyle='None', label='QDHT')\n    plt.title(f'Hankel transform - generalised jinc, order = {order}')\n    plt.xlabel('Frequency /$v$')\n    plt.xlim([0, 1.5])\n    plt.legend()\n    plt.tight_layout()\n\n    error = np.mean(np.abs(expected_ht - actual_ht))\n    assert error < 1e-3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we investigate the function $f(r) = \\\\frac{1}{r^2 + a^2}$,\nthe Hankel transform of which is $K_0(av)$.\n\nNote again the scaling factor of $2\\pi$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = 1\nradius = np.linspace(0, 50, 1024)\ntransformer = HankelTransform(order=0, radial_grid=radius)\nf = 1 / (transformer.r**2 + a**2)\nactual_ht = transformer.qdht(f)\nexpected_ht = 2 * np.pi * scipy_bessel.kn(0, a * transformer.kr)\n\nplt.figure()\nplt.subplot(2, 1, 1)\nplt.title('$\\\\frac{1}{r^2 + a^2}$')\nplt.plot(radius, f)\nplt.xlabel('Radius /$r$')\nplt.xlim([0, 20])\n\nplt.subplot(2, 1, 2)\nplt.plot(kr, expected_ht, label='Analytical')\nplt.plot(kr, actual_ht, marker='x', linestyle='None', label='QDHT')\nplt.title(r'Hankel transform - $2 \\pi K_0(ak)$')\nplt.xlabel('Frequency /$v$')\nplt.xlim([0, 8])\nplt.legend()\nplt.tight_layout()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}