{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Comparison against literature\n\nThis example is a comparison of PyHank against results from the\noriginal publication of the method.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.special as scipybessel\n\nfrom pyhank import HankelTransform"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First we will reproduce figure 1 of\n\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\n.. |Guizar| replace:: Guizar-Sicairos & Guitierrez-Vega\n\nFirst define python functions to calculate the sinc function and its transform\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def sinc(x):\n    return np.sin(x) / x"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Equation 12 of |Guizar|\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def hankel_transform_of_sinc(v):\n    ht = np.zeros_like(v)\n    ht[v < gamma] = (v[v < gamma] ** p * np.cos(p * np.pi / 2)\n                     / (2 * np.pi * gamma * np.sqrt(gamma ** 2 - v[v < gamma] ** 2)\n                        * (gamma + np.sqrt(gamma ** 2 - v[v < gamma] ** 2)) ** p))\n    ht[v >= gamma] = (np.sin(p * np.arcsin(gamma / v[v >= gamma]))\n                      / (2 * np.pi * gamma * np.sqrt(v[v >= gamma] ** 2 - gamma ** 2)))\n    return ht"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now plot the values of the hankel transform and the dynamical error as in figure 1 of |Guizar| `Guizar`_\nfor order 1 and 4\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for p in [1, 4]:\n    transformer = HankelTransform(p, max_radius=3, n_points=256)\n    gamma = 5\n    func = sinc(2 * np.pi * gamma * transformer.r)\n    expected_ht = hankel_transform_of_sinc(transformer.v)\n    ht = transformer.qdht(func)\n    dynamical_error = 20 * np.log10(np.abs(expected_ht - ht) / np.max(ht))\n    not_near_gamma = np.logical_or(transformer.v > gamma * 1.25,\n                                   transformer.v < gamma * 0.75)\n\n    plt.figure()\n    plt.subplot(2, 1, 1)\n    plt.plot(transformer.v, expected_ht, label='Analytical')\n    plt.plot(transformer.v, ht, marker='+', linestyle='None', label='QDHT')\n    plt.title(f'Hankel Transform, p={p}')\n    plt.legend()\n\n    plt.subplot(2, 1, 2)\n    plt.plot(transformer.v, dynamical_error)\n    plt.title('Dynamical error')\n    plt.tight_layout()\n\n    # Check that the error is low, as they do in the paper. Numbers are estimated from their\n    # graphs as they do not quote any for this part\n    assert np.all(dynamical_error < -10)\n    assert np.all(dynamical_error[not_near_gamma] < -35)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we will reproduce figure 3 and confirm we can replicate\nthe errors in the top half of table 1.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "p = 4\na = 1\ntransformer = HankelTransform(order=p, max_radius=2, n_points=1024)\ntop_hat = np.zeros_like(transformer.r)\ntop_hat[transformer.r <= a] = 1\nfunc = transformer.r ** p * top_hat\nexpected_ht = a ** (p + 1) * scipybessel.jv(p + 1, 2 * np.pi * a * transformer.v) / transformer.v\nht = transformer.qdht(func)\n\nretrieved_func = transformer.iqdht(ht)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the overlay as in figure 3 of |Guizar| `Guizar`_\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure()\nplt.subplot(2, 1, 1)\nplt.plot(transformer.v, expected_ht, label='Analytical')\nplt.plot(transformer.v, ht, marker='x', linestyle='None', label='QDHT')\nplt.title(f'Hankel transform $f_2(v)$, order {p}')\nplt.xlabel('Frequency /$v$')\nplt.xlim([0, 10])\nplt.legend()\n\nplt.subplot(2, 1, 2)\nplt.title('Round-trip QDHT vs analytical function')\nplt.plot(transformer.r, func, label='Analytical')\nplt.plot(transformer.r, retrieved_func, marker='x', linestyle='--', label='QDHT+iQDHT')\nplt.xlabel('Radius /$r$')\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now check that the error is the same as that given in Table 1\nof |Guizar| `Guizar`_\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# First calculate e_1 and e_2\nerror_2 = np.mean(np.abs(expected_ht-ht))\nerror_1 = np.mean(np.abs(func-retrieved_func))\nprint(f'Error in Hankel transform is {error_2:.2e}')\nprint(f'Error in reconstructed function is {error_1:.2e}')\n# Note we used 1024 points first\nassert np.isclose(error_2, 4.8e-5, atol=1e-6)\n# Note that Guizar-Sicairos & Guitierrez-Vega got 2.7e-14, so ours is slightly lower\nassert np.isclose(error_1, 2.15e-14, atol=1e-15)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now repeat for 512 points\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "transformer = HankelTransform(order=p, max_radius=2, n_points=512)\ntop_hat = np.zeros_like(transformer.r)\ntop_hat[transformer.r <= a] = 1\nfunc = transformer.r ** p * top_hat\nexpected_ht = a ** (p + 1) * scipybessel.jv(p + 1, 2 * np.pi * a * transformer.v) / transformer.v\nht = transformer.qdht(func)\n\nretrieved_func = transformer.iqdht(ht)\n\nerror_2 = np.mean(np.abs(expected_ht-ht))\nerror_1 = np.mean(np.abs(func-retrieved_func))\nprint(f'Error in Hankel transform is {error_2:.2e}')\nprint(f'Error in reconstructed function is {error_1:.2e}')\n# Note the below is 10 times smaller than\n# #uizar-Sicairos & Guitierrez-Vega (1.3e-3)\nassert np.isclose(error_2, 1.3e-4, atol=1e-5)\nassert np.isclose(error_1, 2.2e-13, atol=1e-14)"
      ]
    }
  ],
  "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.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}