{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cylindrical waves on a grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In its simplest form, we can use a Green's function approach and just evaluate the solution of a single line source. Based on the analytical solution for the 3D Helmholtz equation (as can be read in [this reference](https://users.flatironinstitute.org/~ahb/notes/waveequation.pdf)), we can write\n",
    "\n",
    "$$\n",
    "f(r) = \\frac{e^{ikr}}{4πr} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write some straightforward code to evaluate this on a grid, that I choose to be the same than in the tweet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "wavelength = 1.\n",
    "n_lambda = 20.\n",
    "n_points = 401\n",
    "dx = n_lambda * wavelength / (n_points - 1)\n",
    "k = 2 * np.pi / wavelength\n",
    "print(f\"number of points per wavelenght: {wavelength/dx}\")\n",
    "x = np.linspace(-n_lambda//2 * wavelength, n_lambda//2 * wavelength, num=n_points)\n",
    "y = np.linspace(0, n_lambda * wavelength, num=n_points)\n",
    "X, Y = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's write a function that allows us to compute the field amplitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba as nb\n",
    "\n",
    "@nb.jit()\n",
    "def compute_amplitude_from_source_point(source_location, X, Y, k, phase=0):\n",
    "    x0, y0 = source_location\n",
    "    R = np.sqrt((X - x0)**2 + (Y - y0)**2)\n",
    "    amp = np.exp(1j * (k * R + phase)) \n",
    "    return np.real(amp)\n",
    "\n",
    "amp = compute_amplitude_from_source_point((0., 0.), X, Y, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "hv.extension('matplotlib', logo=False)\n",
    "hv.output(holomap='scrubber')\n",
    "FIG_OPTS = hv.opts(aspect=1.05, fig_inches=5)\n",
    "\n",
    "def plot(amp):\n",
    "    return hv.Image(amp[::-1], bounds=[-10, 0, 10, 20]).opts(colorbar=True, cmap='seismic').redim.label(x='x (wavelengths)', y='y (wavelengths)', z='amplitude').opts(FIG_OPTS)\n",
    "\n",
    "plot(amp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A line of cylindrical sources with phases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the field produced by a set of sources on a grid, we can just add up the resulting fields. Let's write some helper functions to do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_symmetric_point_source(n_one_side, dx, wavelength_solid = 2.5 / 1.5):\n",
    "    pos = np.arange(-n_one_side, n_one_side + 1) * dx\n",
    "    phases = 2 * np.pi / wavelength_solid * pos\n",
    "    return np.c_[pos, np.zeros_like(pos)], phases\n",
    "\n",
    "line_source, phases = make_symmetric_point_source(5, wavelength / 5.)\n",
    "line_source, np.cos(phases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@nb.jit()\n",
    "def compute_amplitude_from_several_points(source_locations, phases, X, Y, k):\n",
    "    field = np.zeros_like(X)\n",
    "    for i in range(len(source_locations)):\n",
    "        source_location, phase = source_locations[i], phases[i]\n",
    "        field += compute_amplitude_from_source_point(source_location, X, Y, k, phase)\n",
    "    return field\n",
    "\n",
    "field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_field(field, line_source, phases):\n",
    "    SCATTER_OPTS = hv.opts(color='z', cmap='gray', padding=0.05, s=100)\n",
    "    return (plot(field) * hv.Scatter(np.c_[line_source, np.cos(phases)], vdims=['y', 'z']).opts(SCATTER_OPTS)).opts(FIG_OPTS)\n",
    "\n",
    "plot_field(field, line_source, phases)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what this looks like with several source points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz_data = {}\n",
    "for n_points in [5, 10, 15, 20, 25, 30, 35, 40, 45]:\n",
    "    line_source, phases = make_symmetric_point_source(n_points, wavelength / 5.)\n",
    "    field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)\n",
    "    viz_data[n_points] = plot_field(field, line_source, phases)\n",
    "    \n",
    "hv.HoloMap(viz_data, kdims='number of points')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The influence of the phase velocity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz_data = {}\n",
    "for c_phi in np.linspace(2, 3, num=10):\n",
    "    wavelength_solid = c_phi / 1.5\n",
    "    line_source, phases = make_symmetric_point_source(80, wavelength / 10., wavelength_solid=wavelength_solid)\n",
    "    field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)\n",
    "    viz_data[c_phi] = plot_field(field, line_source, phases)\n",
    "    \n",
    "hv.HoloMap(viz_data, kdims='phase velocity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Directivity diagram "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A way to show how this works is to make a polar diagram as a function of phase velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive exploration with a DynamicMap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use one of the features of `holoviews`, a `DynamicMap` to explore the parameters interactively: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_with_options(c_phi, n_points, normalize=False):\n",
    "    wavelength_solid = c_phi / 1.5\n",
    "    line_source, phases = make_symmetric_point_source(n_points, wavelength / 10., wavelength_solid=wavelength_solid)\n",
    "    field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)\n",
    "    if normalize:\n",
    "        field /= field.max()\n",
    "    return plot_field(field, line_source, phases)\n",
    "\n",
    "dmap = hv.DynamicMap(plot_with_options, kdims=['c_phi', 'n_points']).redim.range(c_phi=(1.5, 6), \n",
    "                                                                                n_points=(0, 100)).redim.default(n_points=15)\n",
    "# uncomment this for interactive exploration\n",
    "hv.output(dmap, holomap='widgets')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An animation like the one on Twitter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above function returns a plot for each call specifying the focal distance, angle and number of source points. Using repeated calls to the function as well as some interpolation between frames, we can try to remake the animation from the initial tweet with a little guesswork as to what the exact parameters are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "states = [(10, 0), # c_phi,  n_points for first state\n",
    "                  20, # number of frames to interpolate\n",
    "         (10, 100), # c_phi,  n_points for second state\n",
    "                  20, # ...\n",
    "         (1.5, 100),\n",
    "                  20,\n",
    "          (1.5, 2),\n",
    "                  5,\n",
    "          (2.5, 2),\n",
    "                 20,\n",
    "          (2.5, 100),\n",
    "                 20,\n",
    "          (10., 100),\n",
    "                 10,\n",
    "          (10, 0.)\n",
    "         ]\n",
    "\n",
    "def interp(start, end, nframes, frame, isint=False):\n",
    "    val =  frame / (nframes - 1) * (end - start) + start\n",
    "    if isint:\n",
    "        val = int(val)\n",
    "    return val\n",
    "\n",
    "hmap_data = []\n",
    "frame_index = 0\n",
    "for start_state, nframes, end_state in zip(states[::2], states[1::2], states[2::2]):\n",
    "    for frame in tqdm(range(nframes)):\n",
    "        c_phi = interp(start_state[0], end_state[0], nframes, frame)\n",
    "        n_points = interp(start_state[1], end_state[1], nframes, frame, isint=True)\n",
    "        hmap_data.append([frame_index, plot_with_options(c_phi, n_points, normalize=True)])\n",
    "        frame_index += 1\n",
    "\n",
    "# optional: save animation locally\n",
    "hv.save(hv.HoloMap(hmap_data), r'c:\\Temp\\holomap.gif', fps=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}