{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from astropy import units as u\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from plasmapy import formulary, particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physics of the ExB drift"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a single particle of mass $m$ and charge $q$ in a constant, uniform magnetic field $\\mathbf{B}=B\\ \\hat{\\mathbf{z}}$. In the absence of external forces, it travels with velocity $\\mathbf{v}$ governed by the equation of motion\n",
    "\n",
    "$$m\\frac{d\\mathbf{v}}{dt} = q\\mathbf{v}\\times\\mathbf{B}$$\n",
    "\n",
    "which simply equates the net force on the particle to the corresponding Lorentz force. Assuming the particle initially (at time $t=0$) has $\\mathbf{v}$ in the $x,z$ plane (with $v_y=0$), solving reveals \n",
    "\n",
    "$$v_x = v_\\perp\\cos\\omega_c t \\quad\\mathrm{;}\\quad v_y = -\\frac{q}{\\lvert q \\rvert}v_\\perp\\sin \\omega_c t$$\n",
    "\n",
    "while the parallel velocity $v_z$ is constant. This indicates that the particle gyrates in a circular orbit in the $x,y$ plane with constant speed $v_\\perp$, angular frequency $\\omega_c = \\frac{\\lvert q\\rvert B}{m}$, and Larmor radius $r_L=\\frac{v_\\perp}{\\omega_c}$.\n",
    "\n",
    "As an example, take one proton `p+` moving with velocity $1\\ m/s$ in the $x$-direction at $t=0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup proton in uniform B field\n",
    "B = 5 * u.T\n",
    "proton = particles.Particle(\"p+\")\n",
    "omega_c = formulary.parameters.gyrofrequency(B, proton)\n",
    "v_perp = 1 * u.m / u.s\n",
    "r_L = formulary.parameters.gyroradius(B, proton, Vperp=v_perp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can define a function that evolves the particle's position according to the relations above describing $v_x,v_y$, and $v_z$. The option to add a constant drift velocity $v_d$ to the solution is included as an argument, though this drift velocity is zero by default:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def single_particle_trajectory(v_d=np.array([0, 0, 0])):\n",
    "    # Set time resolution & velocity such that proton goes 1 meter along B per rotation\n",
    "    T = 2 * math.pi / omega_c.value  # rotation period\n",
    "    v_parallel = 1 / T * u.m / u.s\n",
    "    dt = T / 1e2 * u.s\n",
    "\n",
    "    # Set initial particle position\n",
    "    x = []\n",
    "    y = []\n",
    "    xt = 0 * u.m\n",
    "    yt = -r_L\n",
    "\n",
    "    # Evolve motion\n",
    "    timesteps = np.arange(0, 10 * T, dt.value)\n",
    "    for t in list(timesteps):\n",
    "        v_x = v_perp * math.cos(omega_c.value * t) + v_d[0]\n",
    "        v_y = v_perp * math.sin(omega_c.value * t) + v_d[1]\n",
    "        xt += +v_x * dt\n",
    "        yt += +v_y * dt\n",
    "        x.append(xt.value)\n",
    "        y.append(yt.value)\n",
    "    x = np.array(x)\n",
    "    y = np.array(y)\n",
    "    z = v_parallel.value * timesteps\n",
    "\n",
    "    return x, y, z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Executing with the default argument and plotting the particle trajectory gives the expected helical motion, with a radius equal to the Larmor radius:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = single_particle_trajectory()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111, projection=\"3d\")\n",
    "ax.plot(x, y, z, label=\"$\\mathbf{F}=0$\")\n",
    "ax.legend()\n",
    "bound = 3 * r_L.value\n",
    "ax.set_xlim([-bound, bound])\n",
    "ax.set_ylim([-bound, bound])\n",
    "ax.set_zlim([0, 10])\n",
    "ax.set_xlabel(\"x [m]\")\n",
    "ax.set_ylabel(\"y [m]\")\n",
    "ax.set_zlabel(\"z [m]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"r_L = {r_L.value:.2e} m\")\n",
    "print(f\"omega_c = {omega_c.value:.2e} rads/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does this motion change when a constant external force $\\mathbf{F}$ is added? The new equation of motion is\n",
    "\n",
    "$$m\\frac{d\\mathbf{v}}{dt} = q\\mathbf{v}\\times\\mathbf{B} + \\mathbf{F}$$\n",
    "\n",
    "and we can find a solution by considering velocities of the form $\\mathbf{v}=\\mathbf{v}_\\parallel + \\mathbf{v}_L + \\mathbf{v}_d$. Here, $\\mathbf{v}_\\parallel$ is the velocity parallel to the magnetic field, $\\mathbf{v}_L$ is the Larmor gyration velocity in the absence of $\\mathbf{F}$ found previously, and $\\mathbf{v}_d$ is some constant drift velocity perpendicular to the magnetic field. Then, we find that\n",
    "\n",
    "$$F_\\parallel = m\\frac{dv_\\parallel}{dt} \\quad\\mathrm{and}\\quad \\mathbf{F}_\\perp = q\\mathbf{B}\\times \\mathbf{v}_d$$\n",
    "\n",
    "and applying the vector triple product yields\n",
    "\n",
    "\n",
    "$$\\mathbf{v}_d = \\frac{1}{q}\\frac{\\mathbf{F}_\\perp\\times\\mathbf{B}}{B^2}$$\n",
    "\n",
    "In the case where the external force $\\mathbf{F} = q\\mathbf{E}$ is due to a constant electric field, this is the constant $\\mathbf{E}\\times\\mathbf{B}$ drift velocity:\n",
    "\n",
    "$$\\boxed{\n",
    "    \\mathbf{v}_d = \\frac{\\mathbf{E}\\times\\mathbf{B}}{B^2}\n",
    "    }$$\n",
    "\n",
    "Built in drift functions allow you to account for the new force added to the system in two different ways:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = 0.2 * u.V / u.m  # E-field magnitude\n",
    "ey = np.array([0, 1, 0])\n",
    "ez = np.array([0, 0, 1])\n",
    "F = proton.charge * E  # force due to E-field\n",
    "\n",
    "v_d = formulary.drifts.force_drift(F * ey, B * ez, proton.charge)\n",
    "print(\"F drift velocity: \", v_d)\n",
    "v_d = formulary.drifts.ExB_drift(E * ey, B * ez)\n",
    "print(\"ExB drift velocity: \", v_d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting particle trajectory can be compared to the case without drifts by calling our previously defined function with the drift velocity now as an argument. As expected, there is a constant drift in the direction of $\\mathbf{E}\\times\\mathbf{B}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_d, y_d, z_d = single_particle_trajectory(v_d=v_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "output-index": 0
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111, projection=\"3d\")\n",
    "ax.plot(x, y, z, label=\"$\\mathbf{F}=0$\")\n",
    "ax.plot(x_d, y_d, z_d, label=\"$\\mathbf{F}=q\\mathbf{E}$\")\n",
    "\n",
    "bound = 3 * r_L.value\n",
    "ax.set_xlim([-bound, bound])\n",
    "ax.set_ylim([-bound, bound])\n",
    "ax.set_zlim([0, 10])\n",
    "ax.set_xlabel(\"x [m]\")\n",
    "ax.set_ylabel(\"y [m]\")\n",
    "ax.set_zlabel(\"z [m]\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"r_L = {r_L.value:.2e} m\")\n",
    "print(f\"omega_c = {omega_c.value:.2e} rads/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, the implementation in our `single_particle_trajectory()` function requires the analytical solution for the velocity $\\mathbf{v}_d$. This solution can be compared with that implemented in the [particle stepper notebook](simulation/particle_stepper.ipynb). It uses the Boris algorithm to evolve the particle along its trajectory in prescribed $\\mathbf{E}$ and $\\mathbf{B}$ fields, and thus does not require the analytical solution."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}