_notebooks/2022-04-28-utd-transition-function.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "95d57ecc",
   "metadata": {},
   "source": [
    "---\n",
    "title: \"Uniform Theory of Diffraction: how to compute the transition function\"\n",
    "math: true\n",
    "categories:\n",
    "  - Research\n",
    "tags:\n",
    "  - propagation\n",
    "  - programming\n",
    "image:\n",
    "    path: /assets/images/misc/McNamara_UTD_book_cover.png\n",
    "    alt: \"Introduction to The Uniform Geometrical Theory of Diffraction - book cover\"\n",
    "---\n",
    "\n",
    "In the Uniform Theory of Diffraction (UTD), diffraction coefficients require evaluating a transition function, \\\\(F(x)\\\\), who has no closed form. However, we will see that its evaluation can be done quite easily with modern scientific libraries.\n",
    "\n",
    "<!--more-->\n",
    "\n",
    "## Context\n",
    "\n",
    "In 1990, McNamara et al. published *Introduction to the Uniform Geometrical Theory of Diffraction* {% cite mcnamaraIntroductionUniformGeometrical1990 %} in which they describe with details how to apply the UTD and how to compute its coefficients. I will not go into details here, but one of the pillars of the UTD is the use of a so-called *transition function*, whose expression is given by\n",
    "\n",
    "$$ F(x) = 2j \\sqrt{x} e^{j x} \\int\\limits_\\sqrt{x}^\\infty e^{-j u^2} \\text{d}u,$$\n",
    "\n",
    "with \\\\(j^2 = -1\\\\) (see Page 184).\n",
    "\n",
    "On the same page, the authors propose approximations for \\\\(F(x)\\\\) for \\\\(x \\le 0.3 \\\\) and \\\\(x \\ge 5.5 \\\\), saying we could interpolate when \\\\(0.3 \\lt x \\lt 5.5 \\\\).\n",
    "\n",
    "However, with modern scientific libraries, we should be able to avoid these approximations and take a more generic approach.\n",
    "\n",
    "### Link to Fresnel integral\n",
    "\n",
    "As mentioned in the book, the integral part in the transition function resembles a Fresnel integral\n",
    "\n",
    "$$\\int\\limits_0^\\infty e^{-j u^2} \\text{d}u = \\sqrt{\\frac{\\pi}{2}} \\frac{1 - j}{2}.$$\n",
    "\n",
    "Then, \\\\(F(x)\\\\) can be rewritten as\n",
    "\n",
    "$$ F(x) = 2j \\sqrt{x} e^{j x} \\Big( \\underbrace{\\int\\limits_0^\\infty e^{-j u^2} \\text{d}u}_{\\sqrt{\\frac{\\pi}{2}} \\frac{1 - j}{2}} - \\int\\limits_0^\\sqrt{x} e^{-j u^2} \\text{d}u\\Big).$$\n",
    "\n",
    "Finally, using [Euler's formula](https://en.wikipedia.org/wiki/Euler%27s_formula), we can further expand\n",
    "\n",
    "$$\\int\\limits_0^\\sqrt{x} e^{-j u^2} \\text{d}u = \\int\\limits_0^\\sqrt{x} \\cos(u^2) -j\\sin(u^2) \\text{d}u = \\underbrace{\\int\\limits_0^\\sqrt{x} \\cos(u^2)\\text{d}u}_{C(\\sqrt{x})} - j \\underbrace{\\int\\limits_0^\\sqrt{x} \\sin(u^2)\\text{d}u}_{S(\\sqrt{x})},$$\n",
    "\n",
    "where \\\\(C(x)\\\\) and \\\\(S(x)\\\\) are also [Fresnel integrals](https://dlmf.nist.gov/7.2#iii), available in many (free) scientific libraries, such as in [`scipy.special.fresnel`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.fresnel.html). If one cannot find an implementation of those functions, it is worth mentioning that \\\\(C(x)\\\\) and \\\\(S(x)\\\\) can be expressed using the [error function](https://en.wikipedia.org/wiki/Fresnel_integral#Properties), which is very likely to be freely available in most programming languages.\n",
    "\n",
    "### Transition function via Fresnel integrals\n",
    "\n",
    "As a summary, \\\\(F(x)\\\\) can be rewritten as\n",
    "\n",
    "$$ 2j \\sqrt{x} e^{j x} \\Big( \\sqrt{\\frac{\\pi}{2}} \\frac{1 - j}{2} - C(\\sqrt{x}) + j S(\\sqrt{x})\\Big).$$\n",
    "\n",
    "## Implementation\n",
    "\n",
    "Here, I will detail how one can implement the transition function in a few lines with Python, using the very popular open source packages [NumPy](https://numpy.org/) and [SciPy](https://scipy.org/).\n",
    "\n",
    "The advantage of using those libraries is that they allow for both scalar and vector inputs, and handle both real and complex numbers.\n",
    "\n",
    "Finally, the plotting will be performed by [matplotlib](https://matplotlib.org/), another very popular open source package.\n",
    "\n",
    "**Note:** SciPy defines \\\\(C(x)\\\\) and \\\\(S(x)\\\\) integrals using \\\\(\\frac{\\pi}{2}u^2\\\\) instead of \\\\(u^2\\\\). As such, the result of the integrals must be multiplied by some **factor** (\\\\(\\sqrt{\\frac{\\pi}{2}}\\\\)) and the argument of \\\\(C(\\sqrt{x})\\\\) and \\\\(S(\\sqrt{x})\\\\) must be divided by this factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "097238ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Package imports\n",
    "\n",
    "import numpy as np\n",
    "import scipy.special as sc\n",
    "from plotting import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfa47b58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transition function\n",
    "\n",
    "\n",
    "def F(x):\n",
    "    factor = np.sqrt(np.pi / 2)\n",
    "    sqrtx = np.sqrt(x)\n",
    "\n",
    "    S, C = sc.fresnel(sqrtx / factor)\n",
    "\n",
    "    return (\n",
    "        2j * sqrtx * np.exp(1j * x) * (factor * ((1 - 1j) / 2 - C + 1j * S))\n",
    "        # We changed the parenthesis so that\n",
    "        # \\sqrt{pi/2} now multiplies C and S\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b338d194",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting\n",
    "\n",
    "x = np.logspace(-3, 1, 100)\n",
    "y = F(x)\n",
    "\n",
    "A = np.abs(y)  # Amplitude of F(x)\n",
    "P = np.angle(y, deg=True)  # Phase (in deg.) of F(x)\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(10, 6.5))\n",
    "\n",
    "ax1.semilogx(x, A, \"k-\")\n",
    "ax1.set_ylabel(\"Magnitude - solid line\")\n",
    "ax2 = plt.twinx()\n",
    "ax2.semilogx(x, P, \"k--\")\n",
    "ax2.set_ylabel(\"Phase (°) - dashed line\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30a89a2a",
   "metadata": {},
   "source": [
    "As a verification, we obtain the same results as in the reference book {% cite mcnamaraIntroductionUniformGeometrical1990 -L figure -l 4.16 %}."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7fdaf87",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "{% bibliography --cited %}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}