{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goal of this notebook:\n",
    "\n",
    "- establish the infinite depth wave equation and dispersion relation\n",
    "- compute the particle trajectories\n",
    "- plot the result of a superposition of a wave coming from the left and a wave coming from the right to imitate the Van Dyke illustration of page 111"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a solution for the over-pressure due to the surface wave, which is valid for the full domain under the hypothesis that the depth of the propagation domain is very large compared to the wavelength of the wave:\n",
    "\n",
    "$$\n",
    "p(x, z, t) = \\Re \\left(\\exp(kz) \\exp(i(kx - \\omega t) \\right)\n",
    "$$\n",
    "\n",
    "If we choose a wavenumber, then we can deduce the frequency from the dispersion relation:\n",
    "$$\n",
    "D(\\omega, k) = \\omega^2 - gk = 0 \\Leftrightarrow \\omega = \\sqrt{ gk } \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First study : a wave propagating from left to right "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given these equation, we can model our first case study: a sinusoidal pressure wave proapgating from left to right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pressure plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first plot the pressure as a function of space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import holoviews as hv\n",
    "hv.extension('matplotlib')\n",
    "\n",
    "qm_opts = hv.opts.QuadMesh(fig_size=300, aspect=2.5, cmap='seismic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-6.25, 6.25, num=400).reshape(1, -1)\n",
    "z = np.linspace(0, -5, num=100).reshape(-1, 1)\n",
    "X, Z = np.meshgrid(x, z)\n",
    "k = 2 * np.pi / 4.\n",
    "rho = 1000\n",
    "\n",
    "def omega_of_k(k):\n",
    "    return np.sqrt(9.81 * k)\n",
    "\n",
    "p = np.real(np.exp(k * z) * np.exp(1j * (k * x - omega_of_k(k) * 0)))\n",
    "\n",
    "hv.QuadMesh((X, Z, p)).opts(qm_opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily vary the value of the wavenumber in this plot to see different behaviours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_of_k(k):\n",
    "    p = np.real(np.exp(k * z) * np.exp(1j * (k * x - omega_of_k(k) * 0)))\n",
    "    return hv.QuadMesh((X, Z, p)).opts(qm_opts)\n",
    "\n",
    "dmap_p = hv.DynamicMap(p_of_k, kdims=['k']).redim.range(k=(1, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(dmap_p[1.] + dmap_p[6]).cols(1).opts(fig_inches=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen above, low wavenumbers (large wavelengths) go deep, while large wavenumbers (small wavelengths) stay at the surface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pressure plot over time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, this pressure wave propagates from left to right due to the change of time. Let's do an animation of the pressure field as a function of its wavenumber."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def animate_pressure(k, Nframes=10):\n",
    "    omega = omega_of_k(k)\n",
    "    T = 2 * np.pi / omega\n",
    "    ts = np.arange(Nframes) / Nframes * T\n",
    "    ps = [np.real(np.exp(k * z) * np.exp(1j * (k * x - omega_of_k(k) * t))) for t in ts]\n",
    "    return {t: hv.QuadMesh((X, Z, p)).opts(qm_opts) for p,t in zip(ps, ts)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber' fig='auto'\n",
    "hmap_p_lowk = hv.HoloMap(animate_pressure(1))\n",
    "hmap_p_lowk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Displacement field plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the question is: what are the two displacement components due to the pressure field? We already know that the surface displacement is proportional to the value of the pressure on the surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def surface_displacement(k):\n",
    "    return hv.Curve((x, np.real(np.exp(1j * (k * x - omega_of_k(k) * 0)))))\n",
    "\n",
    "dmap_y = hv.DynamicMap(surface_displacement, kdims=['k']).redim.range(k=dmap_p.range('k'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(dmap_p * dmap_y)[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "((dmap_p * dmap_y)[1] + (dmap_p * dmap_y)[6]).cols(1).opts(fig_inches=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this explains the displacement at the surface. But what about inside the fluid? It turns out that we can use the momentum balance equations with the known pressure to derive the displacement components and simple \"Fourier maths\":"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left \\{ \\begin{aligned}\n",
    "\\rho \\dfrac{\\partial^2  u}{\\partial t^2 } = - \\dfrac{\\partial p}{\\partial x}\\\\\n",
    "\\rho \\dfrac{\\partial^2  v}{\\partial t^2 } = - \\dfrac{\\partial p}{\\partial z}\n",
    "\\end{aligned} \\right .\n",
    "\\Leftrightarrow\n",
    "\\left \\{ \\begin{aligned}\n",
    "-\\rho \\omega^2 u = - i k p\\\\\n",
    "-\\rho \\omega^2  v = - k p \n",
    "\\end{aligned} \\right .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_solution(k, t=0, nx=28, nz=14):\n",
    "    x = np.linspace(-6.25, 6.25, num=nx).reshape(1, -1)\n",
    "    z = np.linspace(0, -5, num=nz).reshape(-1, 1)\n",
    "    X, Z = np.meshgrid(x, z)\n",
    "    p = np.exp(k * Z) * np.exp(1j * (k * X - omega_of_k(k) * t))\n",
    "    U = 1/(rho * omega_of_k(k) ** 2) * k * p\n",
    "    V = 1/(rho * omega_of_k(k) ** 2) * 1j * k * p\n",
    "    return X, Z, np.real(p), np.real(U), np.real(V)\n",
    "\n",
    "X, Z, p, U, V = compute_solution(k, t=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can turn this computation into a vector field quite easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vector_opts = hv.opts.VectorField(color='Magnitude', \n",
    "                                  magnitude=hv.dim('Magnitude').norm(), \n",
    "                                  rescale_lengths=False,\n",
    "                                 fig_size=300, aspect=2.5)\n",
    "\n",
    "def make_vector_field(X, Z, U, V, vector_opts=None):\n",
    "    mag = np.sqrt(U**2 + V**2)\n",
    "    angle = (np.pi/2.) - np.arctan2(U/mag, V/mag)\n",
    "    field = hv.VectorField((X, Z, angle, mag))\n",
    "    if vector_opts is not None:\n",
    "        field = field.opts(vector_opts)\n",
    "    return field\n",
    "\n",
    "make_vector_field(X, Z, U, V, vector_opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can even animate this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def animate_vector_field(k, Nframes=10):\n",
    "    omega = omega_of_k(k)\n",
    "    T = 2 * np.pi / omega\n",
    "    ts = np.arange(Nframes) / Nframes * T\n",
    "    frames = {}\n",
    "    for t in ts:\n",
    "        X, Z, p, U, V = compute_solution(k, t)\n",
    "        frames[t] = make_vector_field(X, Z, U, V, vector_opts)\n",
    "    return frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does this look at low wavenumbers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber'\n",
    "hv.HoloMap(animate_vector_field(.2, Nframes=20))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Something that is not completely apparent in the above animation is that the trajectories of the particles are little circles. Let's try to superpose the points obtained over several frames on a still image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_trajectory_points(k, Nframes=10, viz_amp=3000, **kwargs):\n",
    "    omega = omega_of_k(k)\n",
    "    T = 2 * np.pi / omega\n",
    "    ts = np.arange(Nframes) / Nframes * T\n",
    "    positions = []\n",
    "    for t in ts:\n",
    "        X, Z, p, U, V = compute_solution(k, t, **kwargs)\n",
    "        positions.append(np.c_[(X + U * viz_amp).flat, (Z + V * viz_amp).flat])\n",
    "    return np.concatenate(positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_opts = hv.opts.Points(color='k', fig_size=300, aspect=2.5, s=np.array([2]).reshape(-1, 1))\n",
    "\n",
    "(hv.Points(compute_trajectory_points(0.1, nx=14, nz=7, Nframes=50)).opts(point_opts) + \\\n",
    "    hv.Points(compute_trajectory_points(1., nx=14, nz=7, Nframes=50)).opts(point_opts)).cols(1).opts(fig_inches=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we see here (for a low wavenumber wave) are the particle displacements plotted over one period of time, and they're all little circles. This does not change for high wavenumbers, however the circle amplitude gets smaller and smaller due to depth."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now put everything together and do an animation of the wave, its particles and the circles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def animate_everything(k, Nframes=20):\n",
    "    trajectory = hv.Points(compute_trajectory_points(k, Nframes=50)).opts(point_opts)\n",
    "    omega = omega_of_k(k)\n",
    "    T = 2 * np.pi / omega\n",
    "    ts = np.arange(Nframes) / Nframes * T\n",
    "    frames = {}\n",
    "    for t in ts:\n",
    "        X, Z, p, U, V = compute_solution(k, t)\n",
    "        frames[t] = hv.QuadMesh((X, Z, p)).opts(qm_opts) * trajectory * make_vector_field(X, Z, U, V, vector_opts)\n",
    "    return frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber'\n",
    "hv.HoloMap(animate_everything(0.4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber'\n",
    "hv.HoloMap(animate_everything(1.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Second study: two waves propagating in opposite directions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now conclude with a recreation of pictures to be found in Milton Van Dyke's *Album of fluid motion* on the chapter on surface waves. I copied three of these together to get the following image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image(filename='files/VanDykeSurfaceWaves.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the previous code to visualize the effect of two waves of opposite directions propagating at the same time.\n",
    "We will trace the resulting pressure field, the vectors as well as the trajectories all on top of each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_solution_two_waves(reflected_part, k, t=0, nx=28, nz=14):\n",
    "    x = np.linspace(-6.25, 6.25, num=nx).reshape(1, -1)\n",
    "    z = np.linspace(0, -5, num=nz).reshape(-1, 1)\n",
    "    X, Z = np.meshgrid(x, z)\n",
    "    p_to_left = np.exp(k * Z) * np.exp(1j * (k * X - omega_of_k(k) * t))\n",
    "    p_to_right = reflected_part * np.exp(k * Z) * np.exp(1j * (k * X + omega_of_k(k) * t))\n",
    "    p = p_to_left + p_to_right\n",
    "    U = 1/(rho * omega_of_k(k) ** 2) * k * p\n",
    "    V = 1/(rho * omega_of_k(k) ** 2) * 1j * k * p\n",
    "    return X, Z, np.real(p), np.real(U), np.real(V)\n",
    "\n",
    "def compute_trajectory_points_two_waves(reflection_part, k, Nframes=10, viz_amp=3000, **kwargs):\n",
    "    omega = omega_of_k(k)\n",
    "    T = 2 * np.pi / omega\n",
    "    ts = np.arange(Nframes) / Nframes * T\n",
    "    positions = []\n",
    "    for t in ts:\n",
    "        X, Z, p, U, V = compute_solution_two_waves(reflection_part, k, t, **kwargs)\n",
    "        positions.append(np.c_[(X + U * viz_amp).flat, (Z + V * viz_amp).flat])\n",
    "    return np.concatenate(positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now do a plot of all this animated information together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def animate_two_waves(k, reflected_part, Nframes=10, viz_amp=4000):\n",
    "    omega = omega_of_k(k)\n",
    "    trajectories = hv.Points(compute_trajectory_points_two_waves(reflected_part, \n",
    "                                                                 k, \n",
    "                                                                 Nframes=100, \n",
    "                                                                 viz_amp=viz_amp)).opts(point_opts).opts(alpha=.3)\n",
    "    T = 2 * np.pi / omega\n",
    "    ts = np.arange(Nframes) / Nframes * T\n",
    "    frames = {}\n",
    "    circle_indices = None\n",
    "    for t in ts:\n",
    "        X, Z, p, U, V = compute_solution_two_waves(reflected_part, k, t)\n",
    "        frames[t] = hv.QuadMesh((X, Z, p)).opts(cmap='plasma', alpha=.2) * \\\n",
    "                            make_vector_field(X, Z, U, V) * \\\n",
    "                            hv.Points(((X + U * viz_amp).flat, \n",
    "                                       (Z + V * viz_amp).flat)).opts(color='black', \n",
    "                                                                     s=np.array([20]).reshape(-1, 1)) * \\\n",
    "                            trajectories\n",
    "    return frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber'\n",
    "hv.HoloMap(animate_two_waves(k=.62, reflected_part=0, Nframes=15), kdims='t').opts(show_frame=False, show_title=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber'\n",
    "hv.HoloMap(animate_two_waves(k=.62, reflected_part=.53, Nframes=15), kdims='t').opts(show_frame=False, show_title=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output holomap='scrubber'\n",
    "hv.HoloMap(animate_two_waves(k=.62, reflected_part=.999, Nframes=15), kdims='t').opts(show_frame=False, show_title=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The comparison to the still pictures I mentioned above are interesting:\n",
    "\n",
    "- circles --> lines that show a standing wave\n",
    "- it's not the same phenomenon: we're using the infinite depth approximation, not the real solution to the experiment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# But there's more ! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They're is always more to learn. Here, my principal source of learning was the MOOC Fundamentals of waves and vibrations, for free on Coursera.\n",
    "Feynman's chapter on waves is also fantastic (http://www.feynmanlectures.caltech.edu/I_51.html)\n",
    "\n",
    "http://courses.washington.edu/mengr543/handouts/Album-Fluid-Motion-Van-Dyke.pdf"
   ]
  }
 ],
 "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.1"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "119px",
    "width": "253px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}