m-thesis-introduction/simulations/01_travelling_1d_sine.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D Travelling Signal (sine with linear dispersion)\n",
    "\n",
    "\n",
    "$$\n",
    "u(x, t) = \\sin\\left( \\frac{2 \\pi x}{\\lambda} -2\\pi f t\\right)\n",
    "\\\\\n",
    "= \\sin\\left( 2 \\pi f \\, ( \\frac{x}{v(f)} -t) \\right)\n",
    "$$\n",
    "\n",
    "$\\lambda = \\dfrac{v}{f}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.fft as ft\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import axes3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TravellingSine:\n",
    "    \"\"\"\n",
    "    Scalar Travelling Sine wave with a constant dispersion relation per frequency\n",
    "    \"\"\"\n",
    "    def __init__(self, x_0, t_0, frequencies):\n",
    "        self.x_0 = x_0\n",
    "        self.t_0 = t_0\n",
    "        self.freqs = np.array(frequencies).reshape(-1)\n",
    "\n",
    "    def __call__(self, x_f, t_f, phase_velocities = None):\n",
    "        x_f = np.array(x_f)\n",
    "        t_f = np.array(t_f)\n",
    "        return self._signal(self.x_0, x_f, self.t_0, t_f, self.freqs, phase_velocities)\n",
    "\n",
    "    def _signal(self, x_0, x_f, t_0, t_f, freqs, phase_velocities = None, debug = 0):\n",
    "\n",
    "        # broadcasting between x_f and t_f\n",
    "        if np.ndim(x_f) != 0 or np.ndim(t_f) != 0:\n",
    "            if np.ndim(x_f) == 0:\n",
    "                x_f = x_f * np.ones(t_f.shape)\n",
    "            if np.ndim(t_f) == 0:\n",
    "                t_f = t_f * np.ones(x_f.shape)\n",
    "        \n",
    "        assert len(x_f) == len(t_f), \"x_f and t_f MUST have equal length: {} and {}\".format(len(x_f), len(t_f))\n",
    "        \n",
    "        # broadcasting between frequencies and phase velocities\n",
    "        if phase_velocities is None:\n",
    "            phase_velocities = 1\n",
    "        \n",
    "        if callable(phase_velocities):\n",
    "            phase_velocities = phase_velocities(freqs)\n",
    "\n",
    "        if np.ndim(phase_velocities) == 0:\n",
    "            phase_velocities = phase_velocities * np.ones(self.freqs.shape)\n",
    "        \n",
    "        assert len(freqs) == len(phase_velocities), \"freqs and phase_velocities MUST have equal length: {} and {}\".format(len(freqs), len(phase_velocities))\n",
    "        \n",
    "        # get the wavelengths for each frequency\n",
    "        wavelengths = phase_velocities / freqs\n",
    "        \n",
    "        # make sure these are 1d arrays\n",
    "        time_diff  = np.array(t_0 - t_f).reshape(-1)        \n",
    "        space_diff = np.array(x_0 - x_f).reshape(-1)       \n",
    "    \n",
    "        # ndarray of phase components\n",
    "        time_phase  = np.outer(2*np.pi*freqs, time_diff)\n",
    "        space_phase = np.outer(2*np.pi/wavelengths, space_diff)\n",
    "        \n",
    "        if debug:\n",
    "            if debug > 1:\n",
    "                #print('time_diff', time_diff)\n",
    "                print('frequencies', freqs)\n",
    "                #print('time_phase', time_phase)\n",
    "\n",
    "                #print('space_diff', space_diff)\n",
    "                print('wavelengths', wavelengths)\n",
    "                #print('space_phase', space_phase)\n",
    "            \n",
    "            print('time_phase shape', time_phase.shape)\n",
    "            print('space_phase shape', space_phase.shape)\n",
    "    \n",
    "        signal = np.sin( time_phase + space_phase )\n",
    "\n",
    "        if freqs.shape:\n",
    "            signal = signal.sum(axis=0)\n",
    " \n",
    "        return signal\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAEECAYAAACLCeeIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eZgkVZ3u/zmRW+VSe3V1bb1V7/RKQ7fdbOLMT4ZhBJSrIqMXUGFklLnMzBUv97lXBB3nOstVRxllVBRHR8AR5+EioI4OLTBCb9JNQ2/VtXTXvlfue5zfH5ERlXtlVlVXdzXxPk8/XZkZeeJEZMQb3/Oe7/c9QkqJCRMmTJhYGCjnuwMmTJgw8XaCSbomTJgwsYAwSdeECRMmFhAm6ZowYcLEAsIkXRMmTJhYQJika8KECRMLCJN0TcwKQoi/EkI8nvp7jRBCpn32SyHEhxe4PyeFEFcv5D5NmJgNrOe7AyZKhxAikPbSBUSBZOr1J6SU/7LwvcqFlPK687DP9Qu9TxMmZgOTdBcRpJQe/W8hRA9wl5TyV4W2F0JYpZSJheibCRMmSoMpL1xESA35nxJCPCGE8AMfEULsEUK8JoSYEkIMCiG+JoSwpbb/jhDiS1ltPCeE+G+pv9uEEP8mhBgVQnQLIT5VYj9eEULcmfr7LiHEb4QQX0n1oUsIcV3atqtT2/tTssQ3ddkiT7uNQojnU+1MCCFeSvusTwhxbdp5eEII8cNUu28KIXakbVvScQkh1gohxoUQIvX6cSHEQNrnTwoh7k07zuOp/XUKIe5K265DCHF92mt7qv9bU6+vTPuNDgshrinlPJtYnDBJ9+LD+4AfAdXAU0ACuA9oAK4Ergc+kdr2R8CH0kilHvg94CkhhAX4GXAAaAXeDdwvhPj9WfTpCuAoUA98BXgs7bMngP9MffZXwEeKtHM/0AUsAZqAzxbZ9r3AD4Aa4AXgawDlHJeUsgOIAFtTb10NRIUQa1OvrwF+k/p7GPgjoAq4G/i6TqqpY7wtrek/BAaklG8IIZYB/w/4HFAHPAD8NPVbmLgIYZLuxYdXpJTPSilVKWVYSnlASrlPSpmQUnYB3wLemdp2L2AD9qRefxB4WUo5DOwGqqSUfy2ljEkpT6OR5Ydm0adOKeV3pZRJ4PtAmxCiQQjRDmwDHkrt4yXguSLtxIEWYHlq+98U2fY3UspfpPb5A2B76v1yj+sl4J1CiFY0Av631Ou1gAN4EyB1zrukhv8Afo1G0qA93N4rhKhIvf7j1HsAtwP/L9VXVUr5c+AI2sPRxEUIk3QvPvSmvxBCbEhJBkNCCB/webSoFymlihYN61HYHwP6ZNwKYHlqyDslhJgCPoMWYZaLobS/Q6n/PWgEOi6lDBfqfxa+BJwBfp0awt9fxj7dqb/LPa7fANeiPaheQntQvVN/LVOOUUKI9wgh9qVkgyngOqbP8wmgE/gjIYQHeA/TpLsCuC2rP7tT58bERQhzIu3iQ7Zt3D8BrwG3SikDQohPo930Op4AfiaE+DKwAy2SA438OqSUG89hXweBeiFEhZQyknpvGVpEmQMppQ/4C+AvhBBbgBeFEPtniHizUe5x/QZN9hgFXgReJiVVpD5DCOEEfoIWLT8npYwLIX4GiLR2dInBBRyWUvak9ed7Uso/LeMYTCximJHuxY9KwAsEhRAbmdZzAZBSHkh9/i3g+RSxAbwKxIQQ/10IUSGEsAghtgghLpuvjkkpO9G03s+lJpeuQtNF80IIcWNq4k2k+pxkOmWuVJR1XFLK46l9fAgtsp0EJoGbmdZzHYAdjZiTQoj3ANka8RNoWu6fMB3lgiZ9vE8I8e5UXyqEEO8SQpiR7kUKk3Qvfvx34A7Ajxb1PpVnmyeA/480Mkilmt0A7AJ6gLHU96vmuX+3oU1IjaNNJj2Fln+cD+uB/wACaJNv/yClfKWcnc3yuF4CRqSUeubCbwAVTXtFSjmFFoH/GzABvB9tsi59v33AQTTp4Mdp7/egTX5+Fo20z6L9Zua9eZFCmCbmJi4kCCGeRht+f+F898WEiXMB82lq4rxCCLFLCLFKCKEIIW5A05ufOd/9MmHiXMGcSDNxvtECPI2Wo9oH3C2lfOP8dsmEiXMHU14wYcKEiQWEKS+YMGHCxALCJF0TJkyYWECYpGvChAkTCwiTdE2YMGFiAWGSrgkTJkwsIEzSNWHChIkFhEm6JkyYMLGAMEnXhAkTJhYQJumaMGHCxALCJF0TJkyYWECYpGvChAkTCwiTdE2YMGFiAWGSrgkTJkwsIEzSNWHChIkFxEx+uqbvowkTJkyUD1HoAzPSNWHChIkFhEm6JkyYMLGAMEnXhAkTJhYQJumaMGHCxAKi7IUp4/E4fX19RCKRc9EfE2WgoqKCtrY2bDbb+e6KCRMmSsRMC1PmfNjd3U1lZSX19fUIUXCCzsQ5hpSS8fFx/H4/q1atOt/dMWHCRCbmL3shEomYhHsBQAhBfX29OeIwYWKRYVaarkm4FwbM38GEicUHcyLNhAkTJhYQFwXp3nDDDUxNTRXd5tprr+XgwYM57x8+fJjnn3/+XHUtL+666y6OHTu2oPs0YcLEhYGysxcuJEgpkVLOiTQPHz7MwYMHueGGG+axZ8Xxne98Z8H2ZcKEiQsLiy7S7enpYePGjXzyk59kx44d9Pb2snLlSsbGxgD4whe+wIYNG3j3u9/Nbbfdxt///d8b3/3Xf/1Xdu3axbp163j55ZeJxWI8+OCDPPXUU2zfvp2nnnoqY1+PP/44733ve7nxxhtZtWoVjzzyCF/+8pe59NJL2b17NxMTE4BG3Lt372br1q28733vY3JykuPHj7Nr166Mfm/duhXIjLo9Hg//63/9L7Zt28bu3bsZHh4GoLOzk927d7Nz504efPBBPB7PuTupJkyYWDDMKdJ9+Nm3ODbgm6++AHBJSxWfu3FT0W1OnjzJ9773Pb7xjW9kvH/w4EGefvppXn/9dRKJBDt27OCyyy4zPk8kEuzfv5/nn3+ehx9+mF/96ld8/vOf5+DBgzzyyCN59/Xmm2/y+uuvE4lEWLNmDX/zN3/D66+/zl/8xV/wz//8z/z5n/85t99+O1//+td55zvfyYMPPsjDDz/MV7/6VWKxGF1dXbS3t/PUU0/xwQ9+MKf9YDDI7t27+eIXv8hnPvMZvv3tb/O///f/5r777uO+++7jtttu49FHH53FmTRhwsSFiEUX6QKsWLGC3bt357z/yiuvcPPNN+N0OqmsrOTGG2/M+PyWW24B4LLLLqOnp6ekfb3rXe+isrKSJUuWUF1dbbS5ZcsWenp68Hq9TE1N8c53vhOAO+64g5deegmAD37wg/z4xz8G4KmnnuLWW2/Nad9ut/Oe97wnp1+vvvoqH/jABwD44z/+45L6utCQUpJMJgmFQvh8PkKhEJFIhHg8TjKZZIYccBMm3paYU6Q7U0R6ruB2u/O+P9NN7nA4ALBYLCQSiZL2pX8HQFEU47WiKDO2ceutt/KBD3yAW265BSEEa9euzdnGZrMZqV/l9Ot8QifbRCKBlBJVVVFVlVgshpQyI5VNURQsFovxT1EUFEUx091MvG2xKCPdQrjqqqt49tlniUQiBAIBnnvuuRm/U1lZid/vn/U+q6urqa2t5eWXXwbgBz/4gRH1rl69GovFwhe+8IW8UW4x7N69m6effhqAJ598ctb9m09IKUkkEkSjUeLxOIBBoEIILBYLVqs1g2CllMTjccLhMG+99RbDw8N4vV58Ph/BYNCMjE287XBRke7OnTu56aab2LZtG7fccguXX3451dXVRb/zrne9i2PHjuWdSCsV3//+97n//vvZunUrhw8f5sEHHzQ+u/XWW/nhD3+YV88thq9+9at8+ctfZteuXQwODs54HOcSOnHqcgpQUrQqhDAiXavVSiwWM74LmsYeiUQIBoP4/X6TjE28LVC298Lx48fZuHHjuevRHBEIBPB4PIRCIa655hq+9a1vsWPHjvPdrbIRCoVwOp0IIXjyySd54okneOaZZ3K2O5e/hx7Z6pLHyZMnWbp0KXV
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "frequencies = 1/np.array([2,-2])\n",
    "velocities = np.array([1, 4])\n",
    "signal  = TravellingSine(0, 0, frequencies)\n",
    "\n",
    "zeros = np.zeros(500)\n",
    "x = np.linspace(-10, 10, 500)\n",
    "t = x\n",
    "\n",
    "x_mesh, t_mesh = np.meshgrid(x,t, sparse=True)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.suptitle(\"Travelling sine wave\")\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot(x, zeros, signal(x, zeros, velocities), label='right moving')\n",
    "ax.plot(x, t, signal(x, t, velocities), label='left moving')\n",
    "ax.plot(x, t, signal(x, t, 1), label='flat line')\n",
    "ax.set_xlabel(\"spatial distance\")\n",
    "ax.set_ylabel('temporal distance')\n",
    "ax.set_zlabel(\"amplitude\")\n",
    "ax.legend()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# animation\n",
    "fig = plt.figure()\n",
    "plt.suptitle(\"Travelling sine wave\")\n",
    "ax = fig.add_subplot(111,projection='3d')\n",
    "ax.set_xlabel(\"spatial distance\")\n",
    "ax.set_ylabel('temporal distance')\n",
    "ax.set_zlabel(\"amplitude\")\n",
    "l1 = ax.plot(x[0], t[0], signal(x[0], t[0], velocities), label='right moving')\n",
    "l2 = ax.plot(x[0], t[0], signal(x[0], t[0], velocities[::-1]), label='left moving')\n",
    "l3 = ax.plot(x[0], t[0], signal(x[0], t[0], 1), label='flat line')\n",
    "\n",
    "lines = [l1, l2, l3]\n",
    "\n",
    "def update(i, fig, ax):\n",
    "    ax.view_init(elev=20., azim=i)\n",
    "    \n",
    "    \n",
    "    ax.legend()\n",
    "    return fig, ax\n",
    " \n",
    "anim = FuncAnimation(fig, update, frames=np.arange(0, 360, 2), repeat=True, fargs=(fig, ax))\n",
    "anim.save('rgb_cube.gif', dpi=80, writer='imagemagick', fps=24)\n",
    "plt.show();"
   ]
  }
 ],
 "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
}