Dynamics of a particle.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6ca7277e",
   "metadata": {},
   "source": [
    "# Dynamics of a particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "063bc1d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import lru_cache\n",
    "\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c198870",
   "metadata": {},
   "source": [
    "## Death Star dynamics (Dynamics & Relativity Sheet 1 Question 4)\n",
    "\n",
    "The equation of motion\n",
    "$$ \\ddot{\\mathbf{r}} = k \\mathbf{r} \\wedge \\dot{\\mathbf{r}} $$\n",
    "could represent the dynamics of a charged particle moving under a magnetic field\n",
    "$$ \\mathbf{B} = \\frac{-km}{q} \\mathbf{r} $$\n",
    "with field lines pointing radially inwards. Hence a particle moving radially (in a straight line and parallel to a magnetic field line) will experience no force, and move with constant velocity.\n",
    "\n",
    "Because the magnetic force is always perpendicular to the direction of velocity, the force does no work on the particle and therefore the speed is constant. \n",
    "\n",
    "If the particle has a velocity component in an angular (non-radial) direction, then the magnetic force will cause the particle to execute a spiralling motion around a radial axis. The frequency of the spiralling and the angular momentum in fact increase with time, since, far from the origin, the magnetic field is very strong and points in a more uniform direction.\n",
    "\n",
    "Note that the high-frequency motion is not very well captured by the following numerical solutions: for example, [aliasing](https://en.wikipedia.org/wiki/Aliasing) is visible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aea415d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun(t, xu, k):\n",
    "    x, u = xu[:3], xu[3:]\n",
    "    return [*u, *(k * np.cross(x, u))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ad040b",
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(\n",
    "    k=widgets.FloatSlider(min=0, max=30, value=6, continuous_update=False),\n",
    "    r0=widgets.FloatSlider(min=0, max=4, value=1, continuous_update=False),\n",
    "    v0_rad=widgets.FloatSlider(min=-2, max=2, value=0.5, continuous_update=False),\n",
    "    v0_perp=widgets.FloatSlider(min=0, max=4, value=2, continuous_update=False),\n",
    ")\n",
    "def death_star_demo(k, r0, v0_rad, v0_perp):\n",
    "    sol = solve_ivp(\n",
    "        fun, (0, 10), \n",
    "        [r0, 0, 0, v0_rad, v0_perp, 0], \n",
    "        args=(k,),\n",
    "#         dense_output=True,\n",
    "        t_eval=np.arange(0, 10, 0.01)\n",
    "    )\n",
    "    assert sol.success, sol.message\n",
    "    xs, ys, zs = sol.y[0, :], sol.y[1, :], sol.y[2, :]\n",
    "    us, vs, ws = sol.y[3, :], sol.y[4, :], sol.y[5, :]\n",
    "\n",
    "    fig = plt.figure(figsize=(14, 8))\n",
    "    ax = fig.add_subplot(121, projection=\"3d\")\n",
    "    ax.plot(xs, ys, zs)\n",
    "    ax.set_xlim((-0.2, 3))\n",
    "    ax.set_ylim((-0.2, 3))\n",
    "    ax.set_zlim((-0.2, 3))\n",
    "\n",
    "    ax = fig.add_subplot(322)\n",
    "    ax.plot(sol.t, xs, 'r', sol.t, ys, 'g', sol.t, zs, 'b')\n",
    "    ax.plot(sol.t, (xs**2 + ys**2 + zs**2)**0.5, 'k')\n",
    "\n",
    "    ax.grid()\n",
    "    ax.legend(('$x$', '$y$', '$z$', r'$|\\mathbf{x}|$',))\n",
    "\n",
    "    ax = fig.add_subplot(324)\n",
    "    ax.plot(sol.t, us, 'r', sol.t, vs, 'g', sol.t, ws, 'b')\n",
    "    ax.plot(sol.t, (us**2 + vs**2 + ws**2)**0.5, 'k')\n",
    "\n",
    "    ax.grid()\n",
    "    ax.legend(('$u$', '$v$', '$w$', r'$|\\dot{\\mathbf{x}}|$',))\n",
    "    \n",
    "    ax = fig.add_subplot(326)\n",
    "    ang_mom = ((ys * ws - zs * vs)**2 + (zs * us - xs * ws)**2 + (xs * vs - ys * us)**2) ** 0.5\n",
    "    ax.plot(sol.t, ang_mom, 'k')\n",
    "#     ax.set_ylim((0, 6))\n",
    "    \n",
    "    ax.legend((r'$|\\mathbf{L}|$',))\n",
    "    ax.grid()\n",
    "#     ax.set_xscale('log')\n",
    "#     ax.set_yscale('log')\n",
    "    \n",
    "    ax.set_xlabel('$t$')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c37f05dd",
   "metadata": {},
   "source": [
    "## Planar polar dynamics\n",
    "\n",
    "In planar polar coordinates $(r, \\theta)$, the equations of motion are\n",
    "$$\n",
    "\\ddot{r} = \\frac{1}{m} f_r + r \\dot{\\theta}^2 , \\qquad\n",
    "\\ddot{\\theta} = \\frac{1}{r} \\left(\\frac{1}{m} f_\\theta - 2\\dot{r}\\dot{\\theta} \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d453f11",
   "metadata": {},
   "source": [
    "## Gravitational field and radiation pressure (Dynamics & Relativity Sheet 2 Question 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2587c4de",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun(t, y, mu):\n",
    "    r, th, r_dot, th_dot = y\n",
    "    \n",
    "    if r < 0.05:\n",
    "        return (0, 0, 0, 0)\n",
    "    \n",
    "    f_r = - 1 / r**2 - mu * r_dot\n",
    "    f_th = - mu * r * th_dot\n",
    "    \n",
    "    r_ddot = f_r + r * th_dot**2\n",
    "    th_ddot = (1/r) * (f_th - 2 * r_dot * th_dot)\n",
    "    return (r_dot, th_dot, r_ddot, th_ddot)\n",
    "\n",
    "\n",
    "tmax = 30\n",
    "\n",
    "@lru_cache()\n",
    "def solver(mu, r0, v0):\n",
    "    sol = solve_ivp(\n",
    "        fun, (0, tmax), [r0, 0, 0, v0],\n",
    "        args=(mu,),\n",
    "        t_eval=np.arange(0, tmax, 0.01),\n",
    "        method='Radau'\n",
    "    )\n",
    "    assert sol.success, sol.message\n",
    "    return sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64b34408",
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(\n",
    "    mu=widgets.FloatSlider(min=0, max=2, value=0.1),\n",
    "    r0=widgets.FloatSlider(min=0.1, max=2, value=1),\n",
    "    v0=widgets.FloatSlider(min=-2, max=2, value=1)\n",
    ")\n",
    "def radiation_pressure_demo(mu, r0, v0):\n",
    "    sol = solver(mu, r0, v0)\n",
    "    rs, ths, r_dots, th_dots = sol.y\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    ax.set_facecolor(\"black\")\n",
    "    ax.plot(\n",
    "        rs * np.cos(ths), rs * np.sin(ths), \n",
    "        'gray'\n",
    "    )\n",
    "    ax.scatter(\n",
    "        rs * np.cos(ths), rs * np.sin(ths), \n",
    "        s=((sol.t-min(sol.t))/(max(sol.t) - min(sol.t)))*30,\n",
    "    #     s=2,\n",
    "        c=sol.t, cmap='hot'\n",
    "    )\n",
    "    ax.set_aspect('equal')\n",
    "\n",
    "    xmin = ymin = -2\n",
    "    xmax = ymax = 2\n",
    "\n",
    "    ax.set_xlim((xmin, xmax))\n",
    "    ax.set_xticks(np.arange(xmin, xmax, 0.5))\n",
    "    ax.set_ylim((ymin, ymax))\n",
    "    ax.set_yticks(np.arange(ymin, ymax, 0.5))\n",
    "    ax.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3981811",
   "metadata": {},
   "source": [
    "## A $1/r^3$ attractive force (Dynamics & Relativity Sheet 2 Question 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b0f2a7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun(t, y, k):\n",
    "    r, th, r_dot, th_dot = y\n",
    "    \n",
    "    if r < 0.05:\n",
    "        return (0, 0, 0, 0)\n",
    "    \n",
    "    f_r = - k / r**3\n",
    "    f_th = 0\n",
    "    \n",
    "    r_ddot = f_r + r * th_dot**2\n",
    "    th_ddot = (1/r) * (f_th - 2 * r_dot * th_dot)\n",
    "    return (r_dot, th_dot, r_ddot, th_ddot)\n",
    "\n",
    "\n",
    "tmax = 30\n",
    "\n",
    "@lru_cache()\n",
    "def solver(r0, v0, k):\n",
    "    sol = solve_ivp(\n",
    "        fun, (0, tmax), [r0, 0, 0, v0],\n",
    "        args=(k,),\n",
    "        t_eval=np.arange(0, tmax, 0.01),\n",
    "        method='Radau'\n",
    "    )\n",
    "    assert sol.success, sol.message\n",
    "    return sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dcd8531",
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(\n",
    "    r0=widgets.FloatSlider(min=0.1, max=2, value=1),\n",
    "    v0=widgets.FloatSlider(min=0, max=2, value=0.975, step=0.025),\n",
    "    k=widgets.FloatSlider(min=0, max=2, value=1),\n",
    ")\n",
    "def inverse_cube_demo(r0, v0, k):\n",
    "    sol = solver(r0, v0, k)\n",
    "    rs, ths, r_dots, th_dots = sol.y\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    ax.set_facecolor(\"black\")\n",
    "    ax.plot(\n",
    "        rs * np.cos(ths), rs * np.sin(ths), \n",
    "        'gray'\n",
    "    )\n",
    "    ax.scatter(\n",
    "        rs * np.cos(ths), rs * np.sin(ths), \n",
    "        s=((sol.t-min(sol.t))/(max(sol.t) - min(sol.t)))*30,\n",
    "    #     s=2,\n",
    "        c=sol.t, cmap='hot'\n",
    "    )\n",
    "    ax.set_aspect('equal')\n",
    "\n",
    "    xmin = ymin = -2\n",
    "    xmax = ymax = 2\n",
    "\n",
    "    ax.set_xlim((xmin, xmax))\n",
    "    ax.set_xticks(np.arange(xmin, xmax, 0.5))\n",
    "    ax.set_ylim((ymin, ymax))\n",
    "    ax.set_yticks(np.arange(ymin, ymax, 0.5))\n",
    "    ax.grid()\n",
    "    plt.show()"
   ]
  }
 ],
 "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}