Improving instructions + added comparison of ptl integration methods.

284142e4 · Jan Benáček · b167502a · 284142e4 · 284142e4
Commit 284142e4 authored Jul 30, 2024 by Jan Benáček
--- a/README.md
+++ b/README.md
@@ -7,6 +7,11 @@
 - fields_mhd.h5: Input datafile from a MHD simulation
 - environment.yml: Downloads and loads the required python packages for the code when Binder stars
+## How to open Jupyter Hub
+- Go to [https://notebooks.mpcdf.mpg.de/isss](https://notebooks.mpcdf.mpg.de/isss)
+- *Sign up first* using your ISSS username and password
+- Login using the same username and password
 ## How to run the code
 - In the left panel, click "particle-basic-motion.ipynb" or "magnetic-reconnection.ipynb"
 - Click Menu "Run", then "Run All Cells" (or click the "fast forward" symbol under the notebook tab)
@@ -32,8 +37,8 @@
  - Use protons instead of electrons
 ## Notes
- Plots can be zoomed in (using the square symbol on the left of each plot)
+- Plots can be zoomed in (using the square symbol on the left of each plot), if you use '%matplotlib widget' in the header of the notebook
- Please avoid to use more than 2 Gb of memory
+- Please avoid to use more than 4 GB of memory
 - If the memory consumption is too high, you can restart everything by clicking the menu "Kernel", then "Restart Kernel and Run All Cells..."
--- a/particle-basic-motion.ipynb
+++ b/particle-basic-motion.ipynb
@@ -254,8 +254,8 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "N = 100 # Number of time steps\n",
+    "N = 1000 # Number of time steps\n",
-    "dt = 0.1 # Time step size in relative units"
+    "dt = 0.01 # Time step size in relative units"
   ]
  },
  {
@@ -363,6 +363,7 @@
   "metadata": {},
   "outputs": [],
   "source": [
+    "plt.figure()\n",
    "plt.plot(trajectory[:,0], trajectory[:,3])\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"v$_x$\")"
@@ -404,8 +405,8 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "N = 200\n",
+    "N = 10000\n",
-    "dt = 0.1\n",
+    "dt = 0.01\n",
    "trajectory = np.zeros((N,6))\n",
    "# Nonzero velocity in x\n",
    "trajectory[0,3] = 1.0"
@@ -431,9 +432,11 @@
   "metadata": {},
   "outputs": [],
   "source": [
+    "plt.figure()\n",
    "plt.plot(trajectory[:,0], trajectory[:,1])\n",
    "plt.xlabel(\"x\")\n",
-    "plt.ylabel(\"x\")"
+    "plt.ylabel(\"x\")\n",
+    "plt.show()"
   ]
  },
  {
@@ -472,7 +475,7 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "N = 200\n",
+    "N = 2000\n",
    "trajectory = np.zeros((N,6))\n",
    "# Nonzero velocity in x\n",
    "trajectory[0,3] = 1.0"
@@ -498,9 +501,11 @@
   "metadata": {},
   "outputs": [],
   "source": [
+    "plt.figure()\n",
    "plt.plot(trajectory[:,0], trajectory[:,1])\n",
    "plt.xlabel(\"x\")\n",
-    "plt.ylabel(\"y\")"
+    "plt.ylabel(\"y\")\n",
+    "plt.show()"
   ]
  },
  {
@@ -539,7 +544,7 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "N = 200\n",
+    "N = 2000\n",
    "trajectory = np.zeros((N,6))\n",
    "# Nonzero velocity in x\n",
    "trajectory[0,3] = 1.0"
@@ -565,9 +570,11 @@
   "metadata": {},
   "outputs": [],
   "source": [
+    "plt.figure()\n",
    "plt.plot(trajectory[:,0], trajectory[:,1])\n",
    "plt.xlabel(\"x\")\n",
-    "plt.ylabel(\"y\")"
+    "plt.ylabel(\"y\")\n",
+    "plt.show()"
   ]
  },
  {
@@ -586,6 +593,86 @@
   "outputs": [],
   "source": []
  },
+  {
+   "cell_type": "markdown",
+   "id": "35af2c11-d661-4445-8823-ff3bad9b12dd",
+   "metadata": {},
+   "source": [
+    "### Various solvers for cyclotron motion"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8710f371-3773-48da-adf3-5eba3a284753",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def E(x,y,z):\n",
+    "    return np.array((0,0,0))\n",
+    "def B(x,y,z):\n",
+    "    return np.array((0,0,1))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3ce66e2-1156-440c-a9b3-7a3e86fcc054",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 200\n",
+    "dt = 0.1\n",
+    "\n",
+    "trajectory1 = np.zeros((N,6))\n",
+    "trajectory2 = np.zeros((N,6))\n",
+    "trajectory3 = np.zeros((N,6))\n",
+    "trajectory4 = np.zeros((N,6))\n",
+    "\n",
+    "# Nonzero velocity in x\n",
+    "trajectory1[0,3] = 1.0\n",
+    "trajectory2[0,3] = 1.0\n",
+    "trajectory3[0,3] = 1.0\n",
+    "trajectory4[0,3] = 1.0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7693d329-25b4-4aa9-b563-97da1e105b5d",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "integrate_particle(Euler_push, dt, trajectory1,\n",
+    "                   q_m_ratio=1, E=E, B=B)\n",
+    "integrate_particle(Runge_Kutta4_push, dt, trajectory2,\n",
+    "                   q_m_ratio=1, E=E, B=B)\n",
+    "integrate_particle(Boris_push, dt, trajectory3,\n",
+    "                   q_m_ratio=1, E=E, B=B)\n",
+    "integrate_particle(Vay_push, dt, trajectory4,\n",
+    "                   q_m_ratio=1, E=E, B=B)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "296550c7-aba3-405e-9018-44fcc2a0211d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(trajectory1[:,0], trajectory1[:,1], \"k\", label=\"Euler\")\n",
+    "plt.plot(trajectory2[:,0], trajectory2[:,1], \"b\", label=\"RK4\")\n",
+    "plt.plot(trajectory3[:,0], trajectory3[:,1], \"g.\", label=\"Boris\")\n",
+    "plt.plot(trajectory4[:,0], trajectory4[:,1], \"r.\", label=\"Vay\")\n",
+    "plt.legend()\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.show()"
+   ]
+  },
  {
   "cell_type": "code",
   "execution_count": null,
@@ -602,6 +689,112 @@
   "outputs": [],
   "source": []
  },
+  {
+   "cell_type": "markdown",
+   "id": "2cb22ef8-8457-435c-b1a4-39f50921ac05",
+   "metadata": {},
+   "source": [
+    "### Various solvers for E x B drift"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3e80ba03-a827-438d-b7c8-9a18576b2351",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def E(x,y,z):\n",
+    "    return np.array((1,0,0))\n",
+    "def B(x,y,z):\n",
+    "    return np.array((0,0,1))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e000d6c6-d8e4-4f18-861a-67f4bb4f0ddb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 200\n",
+    "dt = 0.1\n",
+    "\n",
+    "trajectory1 = np.zeros((N,6))\n",
+    "trajectory2 = np.zeros((N,6))\n",
+    "trajectory3 = np.zeros((N,6))\n",
+    "trajectory4 = np.zeros((N,6))\n",
+    "\n",
+    "# Nonzero velocity in x\n",
+    "trajectory1[0,3] = 1.0\n",
+    "trajectory2[0,3] = 1.0\n",
+    "trajectory3[0,3] = 1.0\n",
+    "trajectory4[0,3] = 1.0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "09198135-b6d4-4b61-b423-c7c5637e644f",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "integrate_particle(Euler_push, dt, trajectory1,\n",
+    "                   q_m_ratio=1, E=E, B=B)\n",
+    "integrate_particle(Runge_Kutta4_push, dt, trajectory2,\n",
+    "                   q_m_ratio=1, E=E, B=B)\n",
+    "integrate_particle(Boris_push, dt, trajectory3,\n",
+    "                   q_m_ratio=1, E=E, B=B)\n",
+    "integrate_particle(Vay_push, dt, trajectory4,\n",
+    "                   q_m_ratio=1, E=E, B=B)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9024774c-d074-4cb7-a841-5933bd6deb44",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(trajectory1[:,0], trajectory1[:,1], \"k\", label=\"Euler\")\n",
+    "plt.plot(trajectory2[:,0], trajectory2[:,1], \"b\", label=\"RK4\")\n",
+    "plt.plot(trajectory3[:,0], trajectory3[:,1], \"g\", label=\"Boris\")\n",
+    "plt.plot(trajectory4[:,0], trajectory4[:,1], \"r--\", label=\"Vay\")\n",
+    "plt.legend()\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6071c66e-df80-4350-af9f-54ac2b74f5d2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(trajectory1[:,0], trajectory1[:,4], \"k\", label=\"Euler\")\n",
+    "plt.plot(trajectory2[:,0], trajectory2[:,4], \"b\", label=\"RK4\")\n",
+    "plt.plot(trajectory3[:,0], trajectory3[:,4], \"g\", label=\"Boris\")\n",
+    "plt.plot(trajectory4[:,0], trajectory4[:,4], \"r--\", label=\"Vay\")\n",
+    "plt.legend()\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7c9eb040-3686-476f-be1e-18489c154bfd",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
  {
   "cell_type": "code",
   "execution_count": null,
@@ -609,6 +802,24 @@
   "metadata": {},
   "outputs": [],
   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "682cf5f3-7f71-4770-a1c1-71af4826984f",
+   "metadata": {},
+   "source": [
+    "## Difference between Boris and Vay pushers\n",
+    "<img src=\"Belyaev+2015-drift.jpg\" width=\"800\">\n",
+    "Belyaev+ (2015)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65cc9503-dc2a-4623-aa81-5b9b7b81dc44",
+   "metadata": {},
+   "outputs": [],
+   "source": []
  }
 ],
 "metadata": {


 %% Cell type:code id:5203c460-72d1-421a-94bb-470659b07b69 tags:
 ``` python
 %matplotlib inline
 #%matplotlib widget
 ```
 %% Cell type:code id:4eb28d9a-fac9-489c-94d5-a3030dbc5385 tags:
 ``` python
 import matplotlib.pyplot as plt
 import numpy as np
 ```
 %% Cell type:code id:1eada8ff-f9eb-4bda-ad62-6ef86ebe2069 tags:
 ``` python
 ```
 %% Cell type:markdown id:14a22f87-8620-4344-bdd2-5019067f7bd4 tags:
 # Basic particle motion
 - Observe particle trajectories and visualization of their motion.
 - Investigate how various pushers impact the precision of the Larmor motion in the magnetic field.
 - Investigate various particle drifts.
 - How the size of time step dt influences the precision of particle position?
 - Which of the particle pushers is the most precise?
 %% Cell type:code id:75f151a0-f464-4824-b3e3-196950d0ce22 tags:
 ``` python
 ```
 %% Cell type:markdown id:f34a8244-4211-400c-bde9-0e0faaf7c52b tags:
 ## Define various particle pushers
 - We want to solve the Lorentz force on a charged particle.
 - Note: All pushers are non-relativistic
 %% Cell type:code id:00ef1b41-99d5-489f-8c2e-a4d5456725d7 tags:
 ``` python
 # Euler method
 def Euler_push(particle, dt, E, B, q_m_ratio):
    velocity = particle[3:]
    E_local = E(particle[0], particle[1], particle[2])
    B_local = B(particle[0], particle[1], particle[2])
    #print("E,B", E_local, B_local)
    vx = velocity[0] + dt * q_m_ratio * (E_local[0] + velocity[1]*B_local[2] - velocity[2]*B_local[1])
    vy = velocity[1] + dt * q_m_ratio * (E_local[1] + velocity[2]*B_local[0] - velocity[0]*B_local[2])
    vz = velocity[2] + dt * q_m_ratio * (E_local[2] + velocity[0]*B_local[1] - velocity[1]*B_local[0])
    #print("New v", vx,vy,vz)
    x = particle[0] + dt * vx
    y = particle[1] + dt * vy
    z = particle[2] + dt * vz
    return np.array((x,y,z, vx,vy,vz))
 ```
 %% Cell type:code id:ad995fee-b938-439d-b9e0-f20bbc32bb77 tags:
 ``` python
 # Kunge-Kutta fourth order method
 def df(particle, dt, E, B, q_m_ratio):
    position = particle[:3]
    velocity = particle[3:]
    # Time dt is useless here because electromagnetic fields do not evolve in time
    E_local = E(position[0], position[1], position[2])
    B_local = B(position[0], position[1], position[2])
    #print("E,B", E_local, B_local)
    dvx = q_m_ratio * (E_local[0] + velocity[1]*B_local[2] - velocity[2]*B_local[1])
    dvy = q_m_ratio * (E_local[1] + velocity[2]*B_local[0] - velocity[0]*B_local[2])
    dvz = q_m_ratio * (E_local[2] + velocity[0]*B_local[1] - velocity[1]*B_local[0])
    return np.array((velocity[0], velocity[1], velocity[2], dvx,dvy,dvz))
 def Runge_Kutta4_push(particle, dt, E, B, q_m_ratio):
    k1 = dt * df(particle, 0, E, B, q_m_ratio)
    k2 = dt * df(particle+k1/2, dt/2, E, B, q_m_ratio)
    k3 = dt * df(particle+k2/2, dt/2, E, B, q_m_ratio)
    k4 = dt * df(particle+k3, dt, E, B, q_m_ratio)
    particle += (k1+2.*k2+2.*k3+k4)/6.
    return particle
 ```
 %% Cell type:code id:151f5867-7828-4424-b1e8-e355a60686f3 tags:
 ``` python
 # Boris method
 def Boris_push(particle, dt, E, B, q_m_ratio):
    E_local = E(particle[0], particle[1], particle[2])
    B_local = B(particle[0], particle[1], particle[2])
    velocity = particle[3:]
    uminus = velocity + q_m_ratio * E_local * dt/2.0
    h1 = q_m_ratio * B_local * dt / 2.0
    uprime = uminus + np.cross(uminus, h1)
    h2 = 2.0 * h1 / (1. + np.dot(h1, h1));
    uplus = uminus + np.cross(uprime, h2);
    u_new = uplus + q_m_ratio * E_local * dt / 2.0
    x = particle[0] + dt * u_new[0]
    y = particle[1] + dt * u_new[1]
    z = particle[2] + dt * u_new[2]
    return np.array((x,y,z, u_new[0],u_new[1],u_new[2]))
 ```
 %% Cell type:code id:92e0171d-e5c4-4219-91af-c26a44246bc4 tags:
 ``` python
 # Vay method
 def Vay_push(particle, dt, E, B, q_m_ratio):
    E_local = E(particle[0], particle[1], particle[2])
    B_local = B(particle[0], particle[1], particle[2])
    velocity = particle[3:]
    vstart = velocity;
    un = vstart + q_m_ratio * dt / 2.0 * (E_local + np.cross(vstart, B_local))
    uprime = un + q_m_ratio * dt / 2.0 * E_local
    tau = q_m_ratio * dt / 2.0 * B_local
    sigma = 1 - np.dot(tau,tau);
    ustar = np.dot(uprime, tau);
    gammanew = np.sqrt(0.5 * (sigma + np.sqrt(sigma*sigma + 4.*(np.dot(tau, tau) + ustar*ustar))))
    t = tau / gammanew
    s = 1./(1. + np.dot(t, t))
    u_new = s * (uprime + np.dot(uprime, t)*t + np.cross(uprime, t));
    x = particle[0] + dt * u_new[0]
    y = particle[1] + dt * u_new[1]
    z = particle[2] + dt * u_new[2]
    return np.array((x,y,z, u_new[0],u_new[1],u_new[2]))
 ```
 %% Cell type:code id:5a009b6f-7dbd-4606-8d21-f2512b68d054 tags:
 ``` python
 ```
 %% Cell type:markdown id:33d2ab58-33d5-4f2f-8e44-31cde05d78c6 tags:
 ### Select your own pusher for the integrations below
 %% Cell type:code id:cb24b7d5-d8ce-45e3-97cd-2144d2ba23df tags:
 ``` python
 my_push = Euler_push
 # my_push = Runge_Kutta4_push
 # my_push = Boris_push
 # my_push = Vay_push
 ```
 %% Cell type:code id:2a4ba614-4ac4-4ca1-9c71-22655b113201 tags:
 ``` python
 ```
 %% Cell type:code id:83211e9a-72df-4792-a0f6-a780a922e17c tags:
 ``` python
 ```
 %% Cell type:markdown id:27d53681-9f73-45a7-8357-1a14c7396429 tags:
 ### Particle path integration
 %% Cell type:code id:89894aaa-f44c-45e7-b5c7-137f3f6e6a42 tags:
 ``` python
 def integrate_particle(pusher, dt,
                       trajectory,
                       q_m_ratio, E, B):
    # Itegrates particle over a path of N time steps
    for i in range(1, trajectory.shape[0]):
        trajectory[i,:] = pusher(trajectory[i-1,:], dt, E, B, q_m_ratio)
    #return trajectory
 ```
 %% Cell type:code id:b3705d0b-f6ee-4973-a048-d2a05c888707 tags:
 ``` python
-N = 100 # Number of time steps
+N = 1000 # Number of time steps
-dt = 0.1 # Time step size in relative units
+dt = 0.01 # Time step size in relative units
 ```
 %% Cell type:code id:79369534-58fe-4a81-b243-0d3b7cc98af9 tags:
 ``` python
 # particle array
 # index 0 - position x
 # index 1 - position y
 # index 2 - position z
 # index 3 - position v_x
 # index 4 - position v_y
 # index 5 - position v_z
 trajectory = np.zeros((N,6))
 ```
 %% Cell type:code id:53a6e5d1-76c1-4dbb-9432-4c2c5de87623 tags:
 ``` python
 ```
 %% Cell type:markdown id:07708ccb-6352-49a6-b1b5-fcee4ece2576 tags:
 # Various particle motions
 %% Cell type:markdown id:14b63166-d55e-4cc4-ad75-2a099cece407 tags:
 ### Acceleration in electric field, no magnetic field
 %% Cell type:code id:50409212-0e86-4ac7-ab54-905c29964da7 tags:
 ``` python
 def E(x,y,z):
    return np.array((1,0,0))
 def B(x,y,z):
    return np.array((0,0,0))
 ```
 %% Cell type:code id:09912138-3dce-48a2-b2b8-0fb89e16bd3c tags:
 ``` python
 N = 100
 trajectory = np.zeros((N,6))
 ```
 %% Cell type:code id:df621d50-da81-45b4-a6ef-c51b1957d6d8 tags:
 ``` python
 integrate_particle(my_push, dt, trajectory,
                   q_m_ratio=1, E=E, B=B)
 ```
 %% Cell type:code id:42414a57-d4b4-4b44-ba8b-c5eb6dabe06b tags:
 ``` python
 trajectory
 ```
 %% Cell type:code id:27da5796-1521-4a55-9b94-6844f0701c0f tags:
 ``` python
 ```
 %% Cell type:code id:e22d552f-a897-45ec-8901-f353b5f3c189 tags:
 ``` python
+plt.figure()
 plt.plot(trajectory[:,0], trajectory[:,3])
 plt.xlabel("x")
 plt.ylabel("v$_x$")
 ```
 %% Cell type:code id:4727f6ea-67b0-4b92-8a52-921a67d55aa4 tags:
 ``` python
 ```
 %% Cell type:markdown id:9adf43e5-24d2-4abc-86a6-fd3c3f9f086f tags:
 ### Cyclotron motion
 %% Cell type:code id:1773f845-661f-476f-b128-c9ee67aeb17d tags:
 ``` python
 def E(x,y,z):
    return np.array((0,0,0))
 def B(x,y,z):
    return np.array((0,0,1))
 ```
 %% Cell type:code id:80755d97-69b0-4253-8f43-09ac3056f2c5 tags:
 ``` python
-N = 200
+N = 10000
-dt = 0.1
+dt = 0.01
 trajectory = np.zeros((N,6))
 # Nonzero velocity in x
 trajectory[0,3] = 1.0
 ```
 %% Cell type:code id:4ea44616-c5fd-47bb-a6e2-91d7194eeb3b tags:
 ``` python
 integrate_particle(my_push, dt, trajectory,
                   q_m_ratio=1, E=E, B=B)
 ```
 %% Cell type:code id:15e2dbc2-039c-4096-b2e7-163e8a1e3939 tags:
 ``` python
+plt.figure()
 plt.plot(trajectory[:,0], trajectory[:,1])
 plt.xlabel("x")
 plt.ylabel("x")
+plt.show()
 ```
 %% Cell type:code id:6ca9d1e8-d607-4483-b6fe-02fd9323e956 tags:
 ``` python
 ```
 %% Cell type:markdown id:56dd66db-80ce-489a-b1bd-f629ff6ddbb1 tags:
 ### E x B drift
 %% Cell type:code id:9974b9ce-d7b2-4833-b257-a88b20e63b11 tags:
 ``` python
 def E(x,y,z):
    return np.array((1,0,0))
 def B(x,y,z):
    return np.array((0,0,1))
 ```
 %% Cell type:code id:f16ceac0-367f-466b-a7e3-475bb9083fc9 tags:
 ``` python
-N = 200
+N = 2000
 trajectory = np.zeros((N,6))
 # Nonzero velocity in x
 trajectory[0,3] = 1.0
 ```
 %% Cell type:code id:e0da1725-c804-40b2-a5e2-739d34bda6e2 tags:
 ``` python
 integrate_particle(my_push, dt, trajectory,
                   q_m_ratio=1, E=E, B=B)
 ```
 %% Cell type:code id:f7aaca0f-db74-4eee-bada-3a74c742ecda tags:
 ``` python
+plt.figure()
 plt.plot(trajectory[:,0], trajectory[:,1])
 plt.xlabel("x")
 plt.ylabel("y")
+plt.show()
 ```
 %% Cell type:code id:5a490875-809b-48be-a007-a39271f44b56 tags:
 ``` python
 ```
 %% Cell type:markdown id:6d8c0328-baf6-4c33-82bc-d1c7af79eb27 tags:
 ### grad B drift
 %% Cell type:code id:f37064f8-d9c5-4c53-8299-5969036006c8 tags:
 ``` python
 def E(x,y,z):
    return np.array((0,0,0))
 def B(x,y,z):
    return np.array((0,0,1+0.1*x))
 ```
 %% Cell type:code id:eaa6da67-08ec-4a6e-a5ef-d53fb8be3efe tags:
 ``` python
-N = 200
+N = 2000
 trajectory = np.zeros((N,6))
 # Nonzero velocity in x
 trajectory[0,3] = 1.0
 ```
 %% Cell type:code id:7f4ab5d0-ec4d-438d-b7c2-c553c44ef2f0 tags:
 ``` python
 integrate_particle(my_push, dt, trajectory,
                   q_m_ratio=1, E=E, B=B)
 ```
 %% Cell type:code id:33462467-76c7-442d-9b08-12a752118df2 tags:
 ``` python
+plt.figure()
 plt.plot(trajectory[:,0], trajectory[:,1])
 plt.xlabel("x")
 plt.ylabel("y")
+plt.show()
 ```
 %% Cell type:code id:fb0e1fd2-3888-4708-8ae9-10660d062d39 tags:
 ``` python
 ```
 %% Cell type:code id:6ebedbc2-ec0a-4d59-b185-3c37c50400f9 tags:
 ``` python
 ```
+%% Cell type:markdown id:35af2c11-d661-4445-8823-ff3bad9b12dd tags:
+### Various solvers for cyclotron motion
+%% Cell type:code id:8710f371-3773-48da-adf3-5eba3a284753 tags:
+``` python
+def E(x,y,z):
+    return np.array((0,0,0))
+def B(x,y,z):
+    return np.array((0,0,1))
+```
+%% Cell type:code id:a3ce66e2-1156-440c-a9b3-7a3e86fcc054 tags:
+``` python
+N = 200
+dt = 0.1
+trajectory1 = np.zeros((N,6))
+trajectory2 = np.zeros((N,6))
+trajectory3 = np.zeros((N,6))
+trajectory4 = np.zeros((N,6))
+# Nonzero velocity in x
+trajectory1[0,3] = 1.0
+trajectory2[0,3] = 1.0
+trajectory3[0,3] = 1.0
+trajectory4[0,3] = 1.0
+```
+%% Cell type:code id:7693d329-25b4-4aa9-b563-97da1e105b5d tags:
+``` python
+integrate_particle(Euler_push, dt, trajectory1,
+                   q_m_ratio=1, E=E, B=B)
+integrate_particle(Runge_Kutta4_push, dt, trajectory2,
+                   q_m_ratio=1, E=E, B=B)
+integrate_particle(Boris_push, dt, trajectory3,
+                   q_m_ratio=1, E=E, B=B)
+integrate_particle(Vay_push, dt, trajectory4,
+                   q_m_ratio=1, E=E, B=B)
+```
+%% Cell type:code id:296550c7-aba3-405e-9018-44fcc2a0211d tags:
+``` python
+plt.figure()
+plt.plot(trajectory1[:,0], trajectory1[:,1], "k", label="Euler")
+plt.plot(trajectory2[:,0], trajectory2[:,1], "b", label="RK4")
+plt.plot(trajectory3[:,0], trajectory3[:,1], "g.", label="Boris")
+plt.plot(trajectory4[:,0], trajectory4[:,1], "r.", label="Vay")
+plt.legend()
+plt.xlabel("x")
+plt.ylabel("y")
+plt.show()
+```
 %% Cell type:code id:88883e35-08e1-4da2-9962-503bc1c2c176 tags:
 ``` python
 ```
 %% Cell type:code id:8295aba9-c108-43fe-87ca-b2e89a011f7c tags:
 ``` python
 ```
+%% Cell type:markdown id:2cb22ef8-8457-435c-b1a4-39f50921ac05 tags:
+### Various solvers for E x B drift
+%% Cell type:code id:3e80ba03-a827-438d-b7c8-9a18576b2351 tags:
+``` python
+def E(x,y,z):
+    return np.array((1,0,0))
+def B(x,y,z):
+    return np.array((0,0,1))
+```
+%% Cell type:code id:e000d6c6-d8e4-4f18-861a-67f4bb4f0ddb tags:
+``` python
+N = 200
+dt = 0.1
+trajectory1 = np.zeros((N,6))
+trajectory2 = np.zeros((N,6))
+trajectory3 = np.zeros((N,6))
+trajectory4 = np.zeros((N,6))
+# Nonzero velocity in x
+trajectory1[0,3] = 1.0
+trajectory2[0,3] = 1.0
+trajectory3[0,3] = 1.0
+trajectory4[0,3] = 1.0
+```
+%% Cell type:code id:09198135-b6d4-4b61-b423-c7c5637e644f tags:
+``` python
+integrate_particle(Euler_push, dt, trajectory1,
+                   q_m_ratio=1, E=E, B=B)
+integrate_particle(Runge_Kutta4_push, dt, trajectory2,
+                   q_m_ratio=1, E=E, B=B)
+integrate_particle(Boris_push, dt, trajectory3,
+                   q_m_ratio=1, E=E, B=B)
+integrate_particle(Vay_push, dt, trajectory4,
+                   q_m_ratio=1, E=E, B=B)
+```
+%% Cell type:code id:9024774c-d074-4cb7-a841-5933bd6deb44 tags:
+``` python
+plt.figure()
+plt.plot(trajectory1[:,0], trajectory1[:,1], "k", label="Euler")
+plt.plot(trajectory2[:,0], trajectory2[:,1], "b", label="RK4")
+plt.plot(trajectory3[:,0], trajectory3[:,1], "g", label="Boris")
+plt.plot(trajectory4[:,0], trajectory4[:,1], "r--", label="Vay")
+plt.legend()
+plt.xlabel("x")
+plt.ylabel("y")
+plt.show()
+```
+%% Cell type:code id:6071c66e-df80-4350-af9f-54ac2b74f5d2 tags:
+``` python
+plt.figure()
+plt.plot(trajectory1[:,0], trajectory1[:,4], "k", label="Euler")
+plt.plot(trajectory2[:,0], trajectory2[:,4], "b", label="RK4")
+plt.plot(trajectory3[:,0], trajectory3[:,4], "g", label="Boris")
+plt.plot(trajectory4[:,0], trajectory4[:,4], "r--", label="Vay")
+plt.legend()
+plt.xlabel("x")
+plt.ylabel("y")
+plt.show()
+```
+%% Cell type:code id:7c9eb040-3686-476f-be1e-18489c154bfd tags:
+``` python
+```
 %% Cell type:code id:42b526e9-ab37-4de9-ab41-c1650071b532 tags:
 ``` python
 ```
+%% Cell type:markdown id:682cf5f3-7f71-4770-a1c1-71af4826984f tags:
+## Difference between Boris and Vay pushers
+<img src="Belyaev+2015-drift.jpg" width="800">
+Belyaev+ (2015)
+%% Cell type:code id:65cc9503-dc2a-4623-aa81-5b9b7b81dc44 tags:
+``` python
+```