# New Algorithms for Simulating Dynamical Friction

## Dan T. Abell, David L. Bruhwiler, Stephen D. Webb — *RadiaSoft, LLC*

This notebook describes—and documents in code—algorithms for simulating  
the dynamical friction experienced by ions in the presence of magnetized electrons.

The $\LaTeX$ preamble is *here*.
$$
%% math text
\newcommand{\hmhsp}{\mspace{1mu}}% math hair space
\newcommand{\mhsp}{\mspace{2mu}}% math hair space
\newcommand{\ud}{\mathop{}\!\mathrm{d}}% upright d for differential
\newcommand{\ui}{\mathrm{i}}% upright i for imaginary unit
\newcommand{\ue}{\mathrm{e}}% upright e for Euler number
%%
\newcommand{\Mion}{m_\text{ion}}
\newcommand{\Me}{m_\text{e}}
%%
\newcommand{\vQion}{\vec{q}_\text{ion}}
\newcommand{\vPion}{\vec{p}_\text{ion}}
\newcommand{\Qion}[1]{#1_\text{ion}}
\newcommand{\Pion}[1]{p_{\text{ion},\hmhsp#1}}
%%
\newcommand{\vQe}{\vec{q}_\text{e}}
\newcommand{\vPe}{\vec{p}_\text{e}}
\newcommand{\Qe}[1]{#1_\text{e}}
\newcommand{\Pe}[1]{p_{\text{e},\hmhsp#1}}
%%
\newcommand{\Map}[2][]{\mathcal{#2}^{#1}}
%%
\newcommand{\pgc}{p_\text{gc}}
\newcommand{\xgc}{x_\text{gc}}
\newcommand{\ygc}{y_\text{gc}}
$$


```python
""" Python preamble """
%matplotlib inline
```


```python
""" Python preamble (cont.) """
from __future__ import division

import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib as mpl

from scipy.constants import pi
from scipy.constants import speed_of_light as clight
from scipy.constants import epsilon_0 as eps0
from scipy.constants import mu_0 as mu0
from scipy.constants import elementary_charge as qe
from scipy.constants import electron_mass as me
from scipy.constants import proton_mass as mp
from scipy.constants import Boltzmann as kB
fourPiEps0 = 4 * pi * eps0
invFourPiEps0 = 1 / fourPiEps0

""" reset some default options """
np.set_printoptions(linewidth=96)

""" indexing """
(Ix, Ipx, Iy, Ipy, Iz, Ipz) = range(6)

""" prefixes """
(femto, pico, nano, micro, milli, one, kilo, mega, giga, tera, peta) = \
    10. ** np.asarray(range(-15, 15+1, 3))
```

We define the ion charge and mass here as global parameters.  
We do the same for the magnetic field strength $B$ and the  
thermal velocity $v_\text{th}$.  

Then we compute various related derived quantities.


```python
"""
angular frequency of Larmor rotations
NB: This is a *signed* quantity, which means
that for electrons, say, you must set Z = -1.
"""
def omega_Larmor(mass, B, Z = 1):
    return Z * qe * B / mass

Z_ion = 1
M_ion = mp
B_mag = 1. # Tesla
e_temp = 300. # Kelvin

N_gyro = 100 # a somewhat arbitrary choice, range [100, 160] 

""" derived quantities """
V_th = math.sqrt(2 * kB * e_temp / me)
rho_gc = me * V_th / (qe * B_mag)
Omega_e = omega_Larmor(me, B_mag, Z = -1)
T_e = (2 * pi) / abs(Omega_e)
T_intxn = N_gyro * T_e

print "V_th = ", V_th
print "rho_gc / µm = ", rho_gc / micro
print "Omega_e / s^(-1) = ", Omega_e
print "frequency / GHz = ", Omega_e / (2 * pi) / giga
print "T_e / ns = ", T_e / nano
print "T_intxn / ns = ", T_intxn / nano
```

    V_th =  95361.4171888
    rho_gc / µm =  0.542189740332
    Omega_e / s^(-1) =  -1.7588200236e+11
    frequency / GHz =  -27.9924900765
    T_e / ns =  0.0357238672682
    T_intxn / ns =  3.57238672682


## Two-body Magnetized Collisions

The Hamiltonian for a two-body interaction between an ion and a magnetized electron is
$$
\vphantom{\Big]}
  H(\vQion, \vPion, \vQe, \vPe)
    = H_0(\vPion, \Qe{y}, \vPe)
    + H_\text{C}(\vQion, \vQe)
$$
where
$$
\begin{align}
  H_0(\vPion, \Qe{y}, \vPe)
    &= \frac{1}{2\Mion}\bigl(\Pion{x}^2 + \Pion{y}^2 + \Pion{z}^2\bigr)
     + \frac{1}{2\Me}\bigl((\Pe{x} + e B \Qe{y})^2 + \Pe{y}^2 + \Pe{z}^2\bigr),\\[1ex]
  H_\text{C}(\vQion, \vQe)
    &= -\frac{Ze^2}{4\pi\varepsilon_0} \big/
                {\sqrt{(\Qion{x}-\Qe{x})^2 + (\Qion{y}-\Qe{y})^2 + (\Qion{z}-\Qe{z})^2}},
\end{align}\\[1ex]
$$
and $e$ denotes the elementary quantum of charge.

The simplest second-order scheme for integrating this system uses  
a *split-operator* approach: We approximate the total map $\Map{M}$ for a  
time step of size $h$ by the symmetric form
$$
\vphantom{\Big]}
  \Map{M}(h) \approx \Map{M}_0(h/2) \Map{M}_C(h) \Map{M}_0(h/2)
$$
where $\Map{M}_0$ and $\Map{M}_C$ are the exact maps for the Hamiltonians $H_0$  
and $H_C$ respectively. The map $\Map{M}_0$ is a simple linear map. The map  
$\Map{M}_C$ generates a nonlinear kick of both ion and electron momenta.

### Hamiltonians for Two-body Magnetized Collisions


```python
"""
Hamiltonian for free ion and electron in a magnetic field, under
the assuption that the ion is unaffected by that magnetic field.

Arguments:
  z_i (ndArray): 6 x N array of canonical coördinates
                 and conjugate momenta for the ions
  z_e (ndArray): 6 x N array of canonical coördinates
                 and conjugate momenta for the electrons
  In both of the above arrays, the six phase-space variables
  are given in the order(x, px, y, py, z, pz)

Return:
  the total 'free' energy of each ion-electron pair
  
"""
def H_twobody_0(z_i, z_e):
    ham_i = ((z_i[Ipx,:] ** 2 + z_i[Ipy,:] ** 2 + z_i[Ipz,:] ** 2)
             / (2 * M_ion))
    ham_e = ((z_e[Ipx,:] + (-qe) * B_mag * z_e[Iy,:]) ** 2
             + z_e[Ipy,:] ** 2 + z_e[Ipz,:] ** 2) / (2 * me)
    return ham_i + ham_e

"""
Hamiltonian for the interaction of each ion-electron pair.
"""
def H_twobody_C(z_i, z_e):
    g_ie = -(Z_ion * qe ** 2) / (4 * pi * eps0)
    intxn = g_ie / np.sqrt(
        + (z_i[Ix,:] - z_e[Ix,:]) ** 2
        + (z_i[Iy,:] - z_e[Iy,:]) ** 2
        + (z_i[Iz,:] - z_e[Iz,:]) ** 2)
    return intxn

"""
Total Hamiltonian for each ion-electron pair.
"""
def H_twobody(z_i, z_e):
    ham_0 = H_twobody_0(z_i, z_e)
    ham_C = H_twobody_C(z_i, z_e)
    return ham_0 + ham_C
```

### Maps for Two-body Magnetized Collisions


```python
"""
define transfer maps for ions and electrons
There are three maps to define here: one each
for ions and electrons under H_0, and another
"""

""" matrix for a linear drift """
def MatD(mass, h):
    Mdrift = np.identity(6)
    for i in (Ix, Iy, Iz):
        Mdrift[i, i + 1] = h / mass
    return Mdrift

""" matrix for linear electron dynamics in a solenoidal field """
def MatK0_e(h):
    mw = me * Omega_e
    wh = Omega_e * h
    cwh = math.cos(wh)
    swh = math.sin(wh)
    cwh1m = 2 * math.sin(wh / 2) ** 2  # 1 - cos(a) = 2 sin^2(a / 2)
    MK0 = np.identity(6)
    MK0[Iy,  Iy ] = cwh
    MK0[Ipy, Ipy] = cwh
    MK0[Iy,  Ipy] = swh / mw
    MK0[Ipy, Iy ] = -mw * swh
    MK0[Iz,  Ipz] = h / me
    MK0[Ix,  Ipx] = swh / mw
    MK0[Ix,  Iy ] = swh
    MK0[Ix,  Ipy] = cwh1m / mw
    MK0[Iy,  Ipx] = -cwh1m / mw
    MK0[Ipy, Ipx] = -swh
    return MK0

"""
map phase-space coördinates forward in time by amount h
based on the Hamiltonian H_0, which describes the free
motion of ions and the motion of electrons in a solenoidal
magnetic field
"""
def MapZ_0(h, z_i, z_e):
    mat = MatD(M_ion, h)
    zf_i = mat.dot(z_i)
    mat = MatK0_e(h)
    zf_e = mat.dot(z_e)
    return zf_i, zf_e

"""
map phase-space coördinates forward in time by amount h
based on the Hamiltonian H_C, which describes the collision
between a single ion-electron pair
"""
def MapZ_C(h, z_i, z_e):
    g = h * Z_ion * qe ** 2 / (4 * pi * eps0)
    dz = z_i - z_e
    denom = (dz[Ix,:] ** 2 + dz[Iy,:] ** 2 + dz[Iz,:] ** 2) ** (3/2)
    zf_i = z_i.copy()
    zf_e = z_e.copy()
    for ip in (Ipx, Ipy, Ipz):
        zf_i[ip,:] = z_i[ip,:] - g * dz[ip - 1] / denom
        zf_e[ip,:] = z_e[ip,:] + g * dz[ip - 1] / denom
    return zf_i, zf_e

def apply_MapZ_0(h, n, z_i, z_e):
    mat_i = MatD(M_ion, h)
    mat_e = MatK0_e(h)
    zf_i = [z_i]
    zf_e = [z_e]
    for i in range(n):
        z_i = mat_i.dot(z_i)
        z_e = mat_e.dot(z_e)
        zf_i.append(z_i)
        zf_e.append(z_e)
    return np.asarray(zf_i), np.asarray(zf_e)

""" second-order split-operator integration for the total Hamiltonian """
def apply_MapZ(h, n, z_i, z_e):
    hh = 0.5 * h
    mat_i = MatD(M_ion, hh)
    mat_e = MatK0_e(hh)
    zf_i = [z_i]
    zf_e = [z_e]
    for i in range(n):
        z_i = mat_i.dot(z_i)
        z_e = mat_e.dot(z_e)
        z_i, z_e = MapZ_C(h, z_i, z_e)
        z_e = mat_e.dot(z_e)
        z_i = mat_i.dot(z_i)
        zf_i.append(z_i)
        zf_e.append(z_e)
    return np.asarray(zf_i), np.asarray(zf_e)
```

## Guiding-center Coördinates and $\Theta$-J Coördinates
### Transformations To and From Guiding-center Coördinates and $\Theta$-J Coördinates

We transform the electron's transverse phase-space coördinates  
using the type-1 generating function
$$
  F_1(x,y;\, \phi,\ygc)
    = m\Omega\Bigl[\frac{1}{2}(y - \ygc)^2\cot\phi - y \ygc\Bigr].
$$
This yields the following transformation rules:  
*to* guiding-center coördinates
$$
\begin{align}
  m\Omega &= qB_0,
               \quad\text{(this is a signed quantity)}\\[1ex]
  \phi    &= \arctan\Bigl(\frac{p_x + e B y}{p_y}\Bigr),\\[1ex]
  p_\phi  &= \frac{1}{2m\Omega}\bigl[(p_x + m\Omega y)^2 + p_y^2\bigr],\\[1ex]
  \ygc    &= -\frac{p_x}{m\Omega},\\[1ex]
  \pgc    &= p_y + m\Omega x.
\end{align}
$$
*from* guiding-center coördinates
$$
\begin{align}
  r_L &= \frac{1}{m\Omega}\sqrt{2m\Omega\,p_\phi},
           \quad\text{(this is a signed quantity)}\\[1ex]
  x   &= \frac{\pgc}{m\Omega} - r_L\cos\phi,\\[1ex]
  p_x &= -m\Omega\,\ygc,\\[1ex]
  y   &= \ygc + r_L\sin\phi,\\[1ex]
  p_y &= m\Omega\,r_L\cos\phi.
\end{align}
$$

We also require the transformation to and from the coördinates $\Theta$-J:
$$
\begin{align}
  \Theta &= \dotsb, \\
  J &= p_\phi + \frac{Ze^2}{4\pi\varepsilon_0} \frac{r_L}{\Omega}
      \frac{(\Qion{x}-\xgc)\cos\phi - (\Qion{y}-\ygc)\sin\phi}{%
        \bigl[(\Qion{x}-\Qe{x})^2 + (\Qion{y}-\Qe{y})^2 + (\Qion{z}-\Qe{z})^2
              + r_L^2\bigr]^{3/2}}.
\end{align}
$$

$$
\begin{align}
  \phi &= \dotsb, \\
  p_\phi &=  \dotsb.
\end{align}
$$


```python
""" convert to guiding-center coordinates """
def toGuidingCenter(z_e):
    mOmega = me * Omega_e
    zgc = z_e.copy()
    zgc[Ix,:]  = np.arctan2(z_e[Ipx,:] + mOmega * z_e[Iy,:], z_e[Ipy,:])
    zgc[Ipx,:] = (((z_e[Ipx,:] + mOmega * z_e[Iy,:]) ** 2
                   + z_e[Ipy,:] ** 2) / (2 * mOmega))
    zgc[Iy,:]  = - z_e[Ipx,:] / mOmega
    zgc[Ipy,:] = z_e[Ipy,:] + mOmega * z_e[Ix,:]
    return zgc

""" convert from guiding-center coordinates """
def fromGuidingCenter(zgc):
    mOmega = me * Omega_e
    rhoL = np.sqrt(2 * mOmega * zgc[Ipx,:]) / mOmega
    z_e = zgc.copy()
    z_e[Ix,:]  = zgc[Ipy,:] / mOmega - rhoL * np.cos(zgc[Ix,:])
    z_e[Ipx,:] = - mOmega * zgc[Iy,:]
    z_e[Iy,:]  = zgc[Iy,:] + rhoL * np.sin(zgc[Ix,:])
    z_e[Ipy,:] = mOmega * rhoL * np.cos(zgc[Ix,:])
    return z_e

""" return J(z_gc) coordinates using the (probably correct) minus sign """
def actionJ(z_i, zgc):
    g = Z_ion * qe ** 2 / (4 * pi * eps0)
    mOmega = me * Omega_e
    rhoL = np.sqrt(2 * zgc[Ipx,:] / mOmega)
    num = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:])
        -  (z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:]))
    den = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
        +  (z_i[Iy,:] - zgc[Iy,:]) ** 2
        +  (z_i[Iz,:] - zgc[Iz,:]) ** 2 + rhoL ** 2) ** (3/2)
    return zgc[Ipx,:] + g * (rhoL / Omega_e) * num / den
```


```python
""" return the Larmor radius """
def rLarmor(z_e):
    mOmega = me * Omega_e
    return np.sqrt((z_e[Ipx,:] + mOmega * z_e[Iy,:]) ** 2 + z_e[Ipy,:] ** 2) / mOmega

""" return the Larmor radius """
def rLarmor_gc(zgc):
    mOmega = me * Omega_e
    return np.sqrt(2 * mOmega * zgc[Ipx,:]) / mOmega

# """ return the quantity 2*r_L[dX cos\phi - dY sin\phi] """
# def minus2RLsq(z_i, zgc):
#     mOmega = me * Omega_e
#     rhoL = np.sqrt(2 * zgc[Ipx,:] / mOmega)
#     res = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:])
#         -  (z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:]))
#     return 2 * rhoL * res
#
# """ return the quantity 2*r_L[dX cos\phi + dY sin\phi] """
# def plus2RLsq(z_i, zgc):
#     mOmega = me * Omega_e
#     rhoL = np.sqrt(2 * zgc[Ipx,:] / mOmega)
#     res = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:])
#         +  (z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:]))
#     return 2 * rhoL * res
#
# """ return J(z_gc) coordinates using the probably wrong plus sign """
# def actionJp(z_i, zgc):
#     g = Z_ion * qe ** 2 / (4 * pi * eps0)
#     mOmega = me * Omega_e
#     rhoL = np.sqrt(2 * zgc[Ipx,:] / mOmega)
#     num = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:])
#         -  (z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:]))
#     den = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
#         +  (z_i[Iy,:] - zgc[Iy,:]) ** 2
#         +  (z_i[Iz,:] - zgc[Iz,:]) ** 2 + rhoL ** 2) ** (3/2)
#     return zgc[Ipx,:] + g * (rhoL / Omega_e) * num / den
#
""" return the perturbation ratio (uses the minus sign) """
def pertubationRatio(z_i, zgc):
    mOmega = me * Omega_e
    rhoL = np.sqrt(2 * mOmega * zgc[Ipx,:]) / mOmega
    num = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:])
           - (z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:]))
    den = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
        +  (z_i[Iy,:] - zgc[Iy,:]) ** 2
        +  (z_i[Iz,:] - zgc[Iz,:]) ** 2 + rhoL ** 2)
    return 2 * rhoL * num / den

# """ return the perturbation ratio using the absolute value """
# def pertubationRatioA(z_i, zgc):
#     mOmega = me * Omega_e
#     rhoL = np.sqrt(2 * mOmega * zgc[Ipx,:]) / mOmega
#     num = (np.abs((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:]))
#            + np.abs((z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:])))
#     den = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
#         +  (z_i[Iy,:] - zgc[Iy,:]) ** 2
#         +  (z_i[Iz,:] - zgc[Iz,:]) ** 2 + rhoL ** 2)
#     return 2 * rhoL * num / den

# """ return the perturbation ratio using the plus sign """
# def pertubationRatioP(z_i, zgc):
#     mOmega = me * Omega_e
#     rhoL = np.sqrt(2 * mOmega * zgc[Ipx,:]) / mOmega
#     num = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * np.cos(zgc[Ix,:])
#            - (z_i[Iy,:] - zgc[Iy,:]) * np.sin(zgc[Ix,:]))
#     den = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
#         +  (z_i[Iy,:] - zgc[Iy,:]) ** 2
#         +  (z_i[Iz,:] - zgc[Iz,:]) ** 2 + rhoL ** 2)
#     return 2 * rhoL * num / den

""" return the ratio (impact parameter) / (Larmor radius) """
def impactParamOverRL(z_i, zgc):
    mOmega = me * Omega_e
    rhoL = np.sqrt(2 * mOmega * zgc[Ipx,:]) / mOmega
    b = np.sqrt((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
               + (z_i[Iy,:] - zgc[Iy,:]) ** 2)
    return b / rhoL
```

### Hamiltonians using Guiding-center Coördinates


```python
"""
Hamiltonian for free ion and electron in a magnetic field, under
the assuption that the ion is unaffected by that magnetic field.
"""
def H_gc_0(z_i, zgc):
    ham_i = ((z_i[Ipx,:] ** 2 + z_i[Ipy,:] ** 2 + z_i[Ipz,:] ** 2)
             / (2 * M_ion))
    ham_e = Omega_e * actionJ(z_i, zgc) + zgc[Ipz,:] ** 2 / (2 * me)
    return ham_i + ham_e

"""
Hamiltonian for the ion-electron interaction
in guiding-center (Θ,J) coördinates.
"""
def H_gc_C(z_i, zgc):
    g_ie = (Z_ion * qe ** 2) / (4 * pi * eps0)
    mOmega = me * Omega_e
    intxn = -g_ie / np.sqrt(
        + (z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
        + (z_i[Iy,:] - zgc[Iy,:]) ** 2
        + (z_i[Iz,:] - zgc[Iz,:]) ** 2
        + 2 * actionJ(z_i, zgc) / mOmega)
    return intxn

"""
total Hamiltonian for the ion-electron system in GC coördinates
"""
def H_gc(z_i, zgc):
    ham_0 = H_gc_0(z_i, zgc)
    ham_C = H_gc_C(z_i, zgc)
    return ham_0 + ham_C

"""
Hamiltonian for the ion-electron interaction
in guiding-center (Θ,J) coördinates.
"""
def H_gc_Cp(z_i, zgc):
    g_ie = (Z_ion * qe ** 2) / (4 * pi * eps0)
    mOmega = me * Omega_e
    intxn = -g_ie / np.sqrt(
        + (z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
        + (z_i[Iy,:] - zgc[Iy,:]) ** 2
        + (z_i[Iz,:] - zgc[Iz,:]) ** 2
        + 2 * actionJp(z_i, zgc) / mOmega)
    return intxn

"""
total Hamiltonian for the ion-electron system in GC coördinates
"""
def H_gc_m(z_i, zgc):
    ham_0 = H_gc_0(z_i, zgc)
    ham_C = H_gc_Cm(z_i, zgc)
    return ham_0 + ham_C
```

### Maps using Guiding-center Coördinates


```python
"""
define transfer maps for ions and electrons
There are three maps to define here: one each
for ions and electrons under H_0, and another
for the interaction under H_c
"""

"""
Map phase-space coördinates forward in time by amount h.
This map is based on the Hamiltonian H_gc_0, which describes
the free motion of ions and the averaged motion of electrons
in a solenoidal magnetic field.
NB: We do not update the \Theta coördinate, as it does not
contribute to the dynamics of any other variables.
"""
def MapZgc_0(h, z_i, zgc):
    matD = MatD(M_ion, h)
    zf_i = matD.dot(z_i)
    zf_e = zgc.copy()
    zf_e[Iz,:] += (h / me) * zgc[Ipz,:]
    return zf_i, zf_e

"""
Map phase-space coördinates forward in time by amount h.
This map is based on the Hamiltonian H_gc_C, which describes
the collision between a single ion-electron pair in guiding-center
coördinates.
NB: We do not update the \Theta coördinate, as it does not
contribute to the dynamics of any other variables.
"""
def MapZgc_C(h, z_i, zgc):
    g_ie = Z_ion * qe ** 2 / (4 * pi * eps0)
    mOmega = me * Omega_e
    dr3 = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) ** 2
           + (z_i[Iy,:] - zgc[Iy,:]) ** 2
           + (z_i[Iz,:] - zgc[Iz,:]) ** 2
           + (2 / mOmega) * actionJ(z_i, zgc)) ** (3/2)
    Omega_gc = (g_ie / mOmega) / dr3
    S = np.sin(Omega_gc * h)
    C1 = 2 * np.sin(Omega_gc * (h / 2)) ** 2
    zf_i = z_i.copy()
    zf_e = zgc.copy()
    Dxgc = ((z_i[Ix,:] - zgc[Ipy,:] / mOmega) * C1
            + (z_i[Iy,:] - zgc[Iy,:]) * S)
    Dygc = ((z_i[Iy,:] - zgc[Iy,:]) * C1
            - (z_i[Ix,:]- zgc[Ipy,:] / mOmega) * S)
    Dpz = (Omega_gc * h) * mOmega * (z_i[Iz,:] - zgc[Iz,:])
    zf_i[Ipx,:] += mOmega * Dygc
    zf_i[Ipy,:] -= mOmega * Dxgc
    zf_i[Ipz,:] -= Dpz
    zf_e[Iy,:]  += Dygc
    zf_e[Ipy,:] += mOmega * Dxgc
    zf_e[Ipz,:] += Dpz
    return zf_i, zf_e

def apply_MapZgc_0(h, n, z_i, zgc):
    mat_i = MatD(M_ion, h)
    mat_e = np.identity(6)
    mat_e[Iz, Ipz] = h / me
    zf_i = [z_i]
    zf_e = [zgc]
    for i in range(n):
        z_i = mat_i.dot(z_i)
        zgc = mat_e.dot(zgc)
        zf_i.append(z_i)
        zf_e.append(zgc)
    return np.asarray(zf_i), np.asarray(zf_e)

def apply_MapZgc(h, n, z_i, zgc):
    hh = 0.5 * h
    mat_i = MatD(M_ion, hh)
    mat_e = np.identity(6)
    mat_e[Iz, Ipz] = hh / me
    zf_i = [z_i]
    zf_e = [zgc]
    for i in range(n):
        z_i = mat_i.dot(z_i)
        zgc = mat_e.dot(zgc)
        z_i, zgc = MapZgc_C(h, z_i, zgc)
        zgc = mat_e.dot(zgc)
        z_i = mat_i.dot(z_i)
        zf_i.append(z_i)
        zf_e.append(zgc) 
    return np.asarray(zf_i), np.asarray(zf_e)
```

### Maps using the Magnus Expansion


```python
"""
compute $\Delta P_\text{ion}$ using the Magnus expansion
"""
def magnus_deltaPIon(h, z_i, zgc):
    g_ie = (Z_ion * qe ** 2) / (4 * pi * eps0)
    mOmega = me * Omega_e
    xgc = zgc[Ipy,:] / mOmega
    C1 = ((z_i[Ix,:] - xgc) ** 2
        + (z_i[Iy,:] - zgc[Iy,:]) ** 2
        + (z_i[Iz,:] - zgc[Iz,:]) ** 2
        + (2 / mOmega) * actionJ(z_i, zgc))
    C2 = (2 * ((z_i[Ix,:] - xgc) * z_i[Ipx,:] / M_ion
             + (z_i[Iy,:] - zgc[Iy,:]) * z_i[Ipy,:] / M_ion
             + (z_i[Iz,:] - zgc[Iz,:]) * (z_i[Ipz,:] / M_ion - zgc[Ipz,:] / me)))
    C3 = ((z_i[Ipx,:] / M_ion) ** 2
        + (z_i[Ipy,:] / M_ion) ** 2
        + (z_i[Ipz,:] / M_ion - zgc[Ipz,:] / me) ** 2)
    B = np.sqrt(C1 + (C2 + C3 * h) * h)
    Delta = 4 * C1 * C3 - C2 ** 2
    D1 = (2 * C3 * h + C2) / B - C2 / np.sqrt(C1)
    D2 = (C2 * h + 2 * C1) / B - 2 * np.sqrt(C1)
    dPx = - ((2 * g_ie / Delta) *
             ((z_i[Ix,:] - xgc) * D1 - (z_i[Ipx,:] / M_ion) * D2))
    dPy = - ((2 * g_ie / Delta) *
             ((z_i[Iy,:] - zgc[Iy,:]) * D1 - (z_i[Ipy,:] / M_ion) * D2))
    dPz = - ((2 * g_ie / Delta) *
             ((z_i[Iz,:] - zgc[Iz,:]) * D1 - (z_i[Ipz,:] / M_ion - zgc[Ipz,:] / me) * D2))
    return np.asarray(( dPx, dPy, dPz)).T
```

### Conversion from (*Q*,*V*) to (*Q*,*P*)


```python
"""
define matrix that transforms ion coördinate-velocity data
to canonically conjugate phase-space variables
"""
QVtoZion = np.identity(6)
for i in (Ipx, Ipy, Ipz):
    QVtoZion[i, i] = M_ion

"""
define matrix that transforms electron coördinate-velocity data
to canonically conjugate phase-space variables
NB: This depends on the local magnetic vector potential,
hence on the local magnetic field.
"""
QVtoZe = np.identity(6)
for i in (Ipx, Ipy, Ipz):
    QVtoZe[i, i] = me
QVtoZe[Ipx, Iy] = (-qe) * (-B_mag)
```

## Test the Algorithms

### Initial conditions


```python
Z_ion = 1
M_ion = mp
e_temp = 300. # Kelvin
B_mag = 1. # Tesla

bOverLarmorR = 20.0 # 20 ==> max perturbation ratio of about 0.1
EkinOverVmax = 10.0 # 10 ==> V_closest_approach / E_kinetic = 0.1
N_gyro = 100    # enforce adequate averaging
tan_alpha = 5.0 # and an adequate opening angle

impactParameterB = micro * np.asarray((1, 2, 5, 10, 22, 50, 100))
# given B = 1 T, those first two values of b are just actually too small

""" derived quantities """
Omega_e = omega_Larmor(me, B_mag, Z = -1)
T_e = (2 * pi) / abs(Omega_e)
""" (the following depend on the impact parameter) """
LarmorR = impactParameterB / bOverLarmorR  # (this version is defined positive)
V_perp = - LarmorR * Omega_e
L_intxn = tan_alpha * impactParameterB #[-1] * np.ones_like(impactParameterB)
V_z = 2 * L_intxn / (N_gyro * T_e)
T_intxn = 2 * L_intxn / V_z
```


```python
print "Omega_e / GHz =", Omega_e / giga
print "T_e / ps     =", T_e / pico
print "T_intxn / ns =", T_intxn / nano
print "L_intxn / um ="
print "  ", L_intxn / micro
print "LarmorR / um ="
print "  ", LarmorR / micro
print "impactParameterB / um ="
print "  ", impactParameterB / micro
print "V_perp / (km/s) ="
print "  ", V_perp / kilo
```

    Omega_e / GHz = -175.88200236
    T_e / ps     = 35.7238672682
    T_intxn / ns = [ 3.57238673  3.57238673  3.57238673  3.57238673  3.57238673  3.57238673  3.57238673]
    L_intxn / um =
       [   5.   10.   25.   50.  105.  250.  500.]
    LarmorR / um =
       [ 0.05  0.1   0.25  0.5   1.05  2.5   5.  ]
    impactParameterB / um =
       [   1.    2.    5.   10.   21.   50.  100.]
    V_perp / (km/s) =
       [   8.79410012   17.58820024   43.97050059   87.94100118  184.67610248  439.7050059   879.4100118 ]



```python
"""
define matrix that transforms ion coördinate-velocity data
to canonically conjugate phase-space variables
"""
QVtoZion = np.identity(6)
for i in (Ipx, Ipy, Ipz):
    QVtoZion[i, i] = M_ion

"""
define matrix that transforms electron coördinate-velocity data
to canonically conjugate phase-space variables
NB: This depends on the local magnetic vector potential,
hence on the local magnetic field.
"""
QVtoZe = np.identity(6)
for i in (Ipx, Ipy, Ipz):
    QVtoZe[i, i] = me
QVtoZe[Ipx, Iy] = (-qe) * (-B_mag)
```


```python
""" initial condition for the ion --- just one, for now """
QVion = np.array([
          ( 0.0 * rho_gc,  0.000 * V_th, 0.0 * rho_gc,  0.000 * V_th, 0.0,  0.000 * V_th)
        ]).transpose()
Zion = QVtoZion.dot(QVion)

""" initial conditions for the electrons """
QVelec = np.asarray([
           np.zeros(impactParameterB.shape),
           V_perp,
           impactParameterB - LarmorR,
           np.zeros(impactParameterB.shape),
           - L_intxn,
           L_intxn * abs(Omega_e) / (pi * N_gyro)
        ])
Zelec = QVtoZe.dot(QVelec)

num_elec0 = Zelec.shape[1]
num_ion0 = Zion.shape[1]
```


```python
"""
=== create arrays of electrons and ions ===
Here we arrange them so that we can pair each ion with each
electron, and compute the \Delta{p} for each interaction.
"""
ZI_elec = np.hstack([Zelec for e in range(num_ion0)])
ZI_ion = Zion[:, np.arange(num_ion0 * num_elec0) // num_elec0]

num_elec = ZI_elec.shape[1]
num_ion = ZI_ion.shape[1]
```




```python
print ZI_elec
print ZI_ion
```

    [[  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
        0.00000000e+00   0.00000000e+00]
     [  1.60217662e-25   3.20435324e-25   8.01088310e-25   1.60217662e-24   3.36457090e-24
        8.01088310e-24   1.60217662e-23]
     [  9.50000000e-07   1.90000000e-06   4.75000000e-06   9.50000000e-06   1.99500000e-05
        4.75000000e-05   9.50000000e-05]
     [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
        0.00000000e+00   0.00000000e+00]
     [ -5.00000000e-06  -1.00000000e-05  -2.50000000e-05  -5.00000000e-05  -1.05000000e-04
       -2.50000000e-04  -5.00000000e-04]
     [  2.54994329e-27   5.09988658e-27   1.27497164e-26   2.54994329e-26   5.35488091e-26
        1.27497164e-25   2.54994329e-25]]
    [[ 0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.]]



```python
num_steps_per_gyro = 40
delta_Te = T_e / num_steps_per_gyro
num_steps = int(np.floor(T_intxn[0] / delta_Te))
print "num_steps =", num_steps
```

    num_steps = 4000



```python
ZF0_i, ZF0_e = apply_MapZ_0(delta_Te, num_steps, ZI_ion, ZI_elec)
```


```python
ZF_i, ZF_e = apply_MapZ(delta_Te, num_steps, ZI_ion, ZI_elec)
```


```python
nsteps = ZF0_e.shape[0]
print nsteps
print num_elec
```

    4001
    7



```python
plt.figure()
plt.xlim(-6, 6)
plt.ylim(-5, 115)
for ptcl in range(num_elec):
    plt.plot(ZF0_e[:, Ix, ptcl] / micro,
             ZF0_e[:, Iy, ptcl] / micro)
    plt.plot((ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e)) / micro,
             (- ZI_elec[Ipx, ptcl] / (me * Omega_e)) / micro,
             'or')
```


```python
plt.figure()
plt.xlim(-6, 6)
plt.ylim(-5, 115)
for ptcl in range(num_elec):
    plt.plot(ZF_e[:, Ix, ptcl] / micro,
             ZF_e[:, Iy, ptcl] / micro)
    plt.plot((ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e)) / micro,
             (- ZI_elec[Ipx, ptcl] / (me * Omega_e)) / micro,
             'or')
```


```python
plt.figure()
for ptcl in range(3):
    plt.plot((ZF_e[:, Ix, ptcl]
              - (ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e))) / micro,
             (ZF_e[:, Iy, ptcl]
              - (- ZI_elec[Ipx, ptcl] / (me * Omega_e))) / micro)
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot((ZF_e[:, Ix, ptcl]
              - (ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e))) / micro,
             (ZF_e[:, Iy, ptcl]
              - (- ZI_elec[Ipx, ptcl] / (me * Omega_e))) / micro)
```


```python
plt.figure()
for ptcl in range(2):
    plt.plot(ZF_e[:, Ix, ptcl] / micro,
             ZF_e[:, Iy, ptcl] / micro)
    plt.plot(ZF0_e[:, Ix, ptcl] / micro,
             ZF0_e[:, Iy, ptcl] / micro)
    plt.plot((ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e)) / micro,
             (- ZI_elec[Ipx, ptcl] / (me * Omega_e)) / micro,
             'ro')
    #plt.plot(0.0, 0.0, 'ro')
```


```python
plt.figure()
for ptcl in range(2,4):
    plt.plot(ZF_e[:, Ix, ptcl] / micro,
             ZF_e[:, Iy, ptcl] / micro)
    plt.plot(ZF0_e[:, Ix, ptcl] / micro,
             ZF0_e[:, Iy, ptcl] / micro)
    plt.plot((ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e)) / micro,
             (- ZI_elec[Ipx, ptcl] / (me * Omega_e)) / micro,
             'ro')
    #plt.plot(0.0, 0.0, 'ro')
```


```python
plt.figure()
for ptcl in range(4,6):
    plt.plot(ZF_e[:, Ix, ptcl] / micro,
             ZF_e[:, Iy, ptcl] / micro)
    plt.plot(ZF0_e[:, Ix, ptcl] / micro,
             ZF0_e[:, Iy, ptcl] / micro)
    plt.plot((ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e)) / micro,
             (- ZI_elec[Ipx, ptcl] / (me * Omega_e)) / micro,
             'ro')
    #plt.plot(0.0, 0.0, 'ro')
```


```python
plt.figure()
for ptcl in range(num_elec - 1, num_elec):
    plt.plot(ZF_e[:, Ix, ptcl] / micro,
             ZF_e[:, Iy, ptcl] / micro)
    plt.plot(ZF0_e[:, Ix, ptcl] / micro,
             ZF0_e[:, Iy, ptcl] / micro)
    plt.plot((ZI_elec[Ix, ptcl] + ZI_elec[Ipy, ptcl] / (me * Omega_e)) / micro,
             (- ZI_elec[Ipx, ptcl] / (me * Omega_e)) / micro,
             'ro')
    #plt.plot(0.0, 0.0, 'ro')
```


```python
plt.figure()
for ptcl in range(2):
    plt.plot((ZF_e - ZF0_e)[:, Ix, ptcl] / nano,
             (ZF_e - ZF0_e)[:, Iy, ptcl] / nano)
```


```python
plt.figure()
for ptcl in range(2,4):
    plt.plot((ZF_e - ZF0_e)[:, Ix, ptcl] / nano,
             (ZF_e - ZF0_e)[:, Iy, ptcl] / nano)
```


```python
plt.figure()
for ptcl in range(4,6):
    plt.plot((ZF_e - ZF0_e)[:, Ix, ptcl] / nano,
             (ZF_e - ZF0_e)[:, Iy, ptcl] / nano)
```


```python
plt.figure()
for ptcl in range(6,7):
    plt.plot((ZF_e - ZF0_e)[:, Ix, ptcl] / nano,
             (ZF_e - ZF0_e)[:, Iy, ptcl] / nano)
```


```python
plt.figure()
for ptcl in range(3):
    plt.plot( ZF_e[:, Ix,  ptcl] / micro,
             (ZF_e[:, Ipx, ptcl] - (-qe) * (-B_mag * ZF_e[:, Iy, ptcl])) / me / kilo)
```


```python
plt.figure()
for ptcl in range(3, num_elec):
    plt.plot( ZF_e[:, Ix,  ptcl] / micro,
             (ZF_e[:, Ipx, ptcl] - (-qe) * (-B_mag * ZF_e[:, Iy, ptcl])) / me / kilo)
```


```python
plt.figure()
for ptcl in range(3):
    plt.plot((ZF_e[:, Iy,  ptcl] - (- ZI_elec[Ipx, ptcl] / (me * Omega_e))) / nano,
             ZF_e[:, Ipy, ptcl] / me / kilo)
```


```python
plt.figure()
for ptcl in range(3, num_elec):
    plt.plot((ZF_e[:, Iy,  ptcl] - (- ZI_elec[Ipx, ptcl] / (me * Omega_e))) / nano,
             ZF_e[:, Ipy, ptcl] / me / kilo)
```


```python
plt.figure()
for ptcl in range(2):
    plt.plot( (ZF_e - ZF0_e)[:, Ix,  ptcl] / nano,
             ((ZF_e - ZF0_e)[:, Ipx, ptcl]
              - (-qe) * (-B_mag * (ZF_e - ZF0_e)[:, Iy, ptcl])) / me / one)
```


```python
plt.figure()
for ptcl in range(2,4):
    plt.plot( (ZF_e - ZF0_e)[:, Ix,  ptcl] / nano,
             ((ZF_e - ZF0_e)[:, Ipx, ptcl]
              - (-qe) * (-B_mag * (ZF_e - ZF0_e)[:, Iy, ptcl])) / me / one)
```


```python
plt.figure()
for ptcl in range(4,6):
    plt.plot( (ZF_e - ZF0_e)[:, Ix,  ptcl] / nano,
             ((ZF_e - ZF0_e)[:, Ipx, ptcl]
              - (-qe) * (-B_mag * (ZF_e - ZF0_e)[:, Iy, ptcl])) / me / one)
```


```python
plt.figure()
for ptcl in range(6,7):
    plt.plot( (ZF_e - ZF0_e)[:, Ix,  ptcl] / nano,
             ((ZF_e - ZF0_e)[:, Ipx, ptcl]
              - (-qe) * (-B_mag * (ZF_e - ZF0_e)[:, Iy, ptcl])) / me / milli)
```


```python
plt.figure()
nst = -1
for ptcl in range(2):
    plt.plot((ZF_e - ZF0_e)[:nst, Iy,  ptcl] / nano,
             (ZF_e - ZF0_e)[:nst, Ipy, ptcl] / me / one)
```


```python
plt.figure()
nst = -1
for ptcl in range(2,4):
    plt.plot((ZF_e - ZF0_e)[:nst, Iy,  ptcl] / nano,
             (ZF_e - ZF0_e)[:nst, Ipy, ptcl] / me / one)
```


```python
plt.figure()
nst = -1
for ptcl in range(4,6):
    plt.plot((ZF_e - ZF0_e)[:nst, Iy,  ptcl] / pico,
             (ZF_e - ZF0_e)[:nst, Ipy, ptcl] / me / one)
```


```python
plt.figure()
nst = -1
for ptcl in range(6,7):
    plt.plot((ZF_e - ZF0_e)[:nst, Iy,  ptcl] / pico,
             (ZF_e - ZF0_e)[:nst, Ipy, ptcl] / me / milli)
```


```python
pRatio = np.asarray(map(lambda zi, ze: pertubationRatio(zi, toGuidingCenter(ze)), ZF_i, ZF_e))
bDivR = np.asarray(map(lambda zi, ze: impactParamOverRL(zi, toGuidingCenter(ze)), ZF_i, ZF_e))
```


```python
print pRatio.shape
print bDivR.shape
print nsteps
```

    (4001, 7)
    (4001, 7)
    4001



```python
plt.figure()
n0 = 0
n1 = nsteps
times = delta_Te * np.asarray(range(nsteps))
for ptcl in range(2):
    plt.plot(times[n0:n1],
             pRatio[n0:n1, ptcl])
```


```python
plt.figure()
n0 = 0
n1 = nsteps
times = delta_Te * np.asarray(range(nsteps))
for ptcl in range(2,4):
    plt.plot(times[n0:n1],
             pRatio[n0:n1, ptcl])
```


```python
plt.figure()
n0 = 0
n1 = nsteps
times = delta_Te * np.asarray(range(nsteps))
for ptcl in range(4,6):
    plt.plot(times[n0:n1],
             pRatio[n0:n1, ptcl])
```


```python
plt.figure()
n0 = 0
n1 = nsteps
times = delta_Te * np.asarray(range(nsteps))
for ptcl in range(6,7):
    plt.plot(times[n0:n1],
             pRatio[n0:n1, ptcl])
```


```python
plt.figure()
for ptcl in range(2):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             -bDivR[:, ptcl])
```


```python
plt.figure()
for ptcl in range(2,4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             20 - np.abs(bDivR[:, ptcl]))
```


```python
plt.figure()
for ptcl in range(4,6):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             20 - np.abs(bDivR[:, ptcl]))
```


```python
plt.figure()
for ptcl in range(6,7):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             20 - np.abs(bDivR[:, ptcl]))
```


```python
# Delta-P_ion, calculated from the explicit map
(ZF_i[-1, [Ipx, Ipy, Ipz], :] - ZF_i[0, [Ipx, Ipy, Ipz], :]).T
```




    array([[  1.08030578e-26,   5.72808106e-26,   1.45007842e-27],
           [  4.96361379e-28,   1.78198085e-26,   1.73582683e-27],
           [  2.00967145e-29,   5.67317879e-27,   6.60213989e-29],
           [  7.85620796e-31,   1.58638251e-27,   2.60010730e-30],
           [  1.98872502e-32,   3.65764404e-28,   6.59269280e-32],
           [  2.60902308e-34,   6.46440746e-29,   8.65070568e-34],
           [  8.15538018e-36,   1.61631831e-29,   2.70410697e-35]])




```python
ZI_elec_gc = toGuidingCenter(ZI_elec)
ZFgc_i, ZFgc_e = apply_MapZgc(delta_Te, num_steps, ZI_ion, ZI_elec_gc)
```


```python
print (ZFgc_i[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i[0, [Ipx, Ipy, Ipz], :]).T
```

    [[  1.06736119e-26   5.69619435e-26   1.21445745e-27]
     [  4.93897920e-28   1.77771766e-26   1.73546012e-27]
     [  1.99971566e-29   5.65986907e-27   6.59754096e-29]
     [  7.81449851e-31   1.58239983e-27   2.59765932e-30]
     [  1.97803771e-32   3.64835112e-28   6.58618687e-32]
     [  2.59498236e-34   6.44796043e-29   8.64212236e-34]
     [  8.11148695e-36   1.61220560e-29   2.70142291e-35]]



```python
print "T_intxn = ", T_intxn
print "T_intxn / delta_Te = ", T_intxn / delta_Te
print "num_steps = ", num_steps
print "num_steps_per_gyro = ", num_steps_per_gyro
```

    T_intxn =  [  3.57238673e-09   3.57238673e-09   3.57238673e-09   3.57238673e-09   3.57238673e-09
       3.57238673e-09   3.57238673e-09]
    T_intxn / delta_Te =  [ 4000.  4000.  4000.  4000.  4000.  4000.  4000.]
    num_steps =  4000
    num_steps_per_gyro =  40



```python
ZFgc_i_x1, ZFgc_e_x1 = apply_MapZgc(40 * delta_Te, num_steps // 40, ZI_ion, ZI_elec_gc)
print (ZFgc_i_x1[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x1[0, [Ipx, Ipy, Ipz], :]).T
```

    [[  1.05836708e-26   5.67504416e-26   8.51768805e-28]
     [  4.92999849e-28   1.77610316e-26   1.73512631e-27]
     [  1.99946520e-29   5.65951484e-27   6.59695843e-29]
     [  7.81439841e-31   1.58238970e-27   2.59745819e-30]
     [  1.97804149e-32   3.64835461e-28   6.58568793e-32]
     [  2.59499163e-34   6.44797194e-29   8.64146934e-34]
     [  8.11151685e-36   1.61220857e-29   2.70121883e-35]]



```python
ZFgc_i_x2, ZFgc_e_x2 = apply_MapZgc(2 * 40*delta_Te, num_steps//40 // 2, ZI_ion, ZI_elec_gc)
print (ZFgc_i_x2[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x2[0, [Ipx, Ipy, Ipz], :]).T
```

    [[  2.14817657e-27   2.61350551e-26  -2.18295742e-27]
     [  4.90448210e-28   1.77150755e-26   1.74567430e-27]
     [  1.99870730e-29   5.65844276e-27   6.59521178e-29]
     [  7.81409557e-31   1.58235906e-27   2.59685491e-30]
     [  1.97805278e-32   3.64836502e-28   6.58419133e-32]
     [  2.59501942e-34   6.44800646e-29   8.63951055e-34]
     [  8.11160656e-36   1.61221749e-29   2.70060664e-35]]



```python
ZFgc_i_x4, ZFgc_e_x4 = apply_MapZgc(4 * 40*delta_Te, num_steps//40 // 4, ZI_ion, ZI_elec_gc)
print (ZFgc_i_x4[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x4[0, [Ipx, Ipy, Ipz], :]).T
```

    [[  5.20096310e-26   1.16988782e-25   1.29139359e-26]
     [  4.73388458e-28   1.74041833e-26   2.62007515e-27]
     [  1.99553678e-29   5.65395568e-27   6.56948668e-29]
     [  7.81288077e-31   1.58223614e-27   2.59691645e-30]
     [  1.97810315e-32   3.64841150e-28   6.58241371e-32]
     [  2.59513811e-34   6.44815393e-29   8.63690068e-34]
     [  8.11198910e-36   1.61225550e-29   2.69978481e-35]]



```python
ZFgc_i_x5, ZFgc_e_x5 = apply_MapZgc(5 * 40*delta_Te, num_steps//40 // 5, ZI_ion, ZI_elec_gc)
print (ZFgc_i_x5[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x5[0, [Ipx, Ipy, Ipz], :]).T
```

    [[  4.79887796e-27   3.89088607e-26  -1.65627962e-26]
     [  5.34180299e-28   1.84878855e-26   1.19444709e-27]
     [  1.99450155e-29   5.65249092e-27   6.68397054e-29]
     [  7.81119753e-31   1.58206575e-27   2.55524304e-30]
     [  1.97799548e-32   3.64831221e-28   6.50514063e-32]
     [  2.59504580e-34   6.44803924e-29   8.53964576e-34]
     [  8.11171128e-36   1.61222789e-29   2.66947425e-35]]



```python
ZFgc_i_x10, ZFgc_e_x10 = apply_MapZgc(10 * 40*delta_Te, num_steps//40 // 10, ZI_ion, ZI_elec_gc)
print (ZFgc_i_x10[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x10[0, [Ipx, Ipy, Ipz], :]).T
```

    [[  2.23422346e-28   8.45718951e-27  -1.53901390e-26]
     [  1.00821652e-26   7.93120021e-26   1.42130883e-26]
     [  1.98058416e-29   5.63275859e-27  -1.41768453e-28]
     [  7.57240857e-31   1.55769655e-27  -3.26299395e-30]
     [  1.92800945e-32   3.60191899e-28  -6.12956798e-32]
     [  2.53116752e-34   6.36818402e-29  -7.71914090e-34]
     [  7.91241310e-36   1.59229916e-29  -2.40580996e-35]]



```python
(ZFgc_i[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i[0, [Ipx, Ipy, Ipz], :]).T
```




    array([[  1.06736119e-26,   5.69619435e-26,   1.21445745e-27],
           [  4.93897920e-28,   1.77771766e-26,   1.73546012e-27],
           [  1.99971566e-29,   5.65986907e-27,   6.59754096e-29],
           [  7.81449851e-31,   1.58239983e-27,   2.59765932e-30],
           [  1.97803771e-32,   3.64835112e-28,   6.58618687e-32],
           [  2.59498236e-34,   6.44796043e-29,   8.64212236e-34],
           [  8.11148695e-36,   1.61220560e-29,   2.70142291e-35]])




```python
(ZFgc_i_x5[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x5[0, [Ipx, Ipy, Ipz], :]).T
```




    array([[  4.79887796e-27,   3.89088607e-26,  -1.65627962e-26],
           [  5.34180299e-28,   1.84878855e-26,   1.19444709e-27],
           [  1.99450155e-29,   5.65249092e-27,   6.68397054e-29],
           [  7.81119753e-31,   1.58206575e-27,   2.55524304e-30],
           [  1.97799548e-32,   3.64831221e-28,   6.50514063e-32],
           [  2.59504580e-34,   6.44803924e-29,   8.53964576e-34],
           [  8.11171128e-36,   1.61222789e-29,   2.66947425e-35]])




```python
magnus_deltaPIon(T_intxn, ZI_ion, ZI_elec_gc)
```




    array([[ -0.00000000e+00,   1.61222627e-25,  -1.01894869e-40],
           [ -0.00000000e+00,   4.03058781e-26,  -0.00000000e+00],
           [ -0.00000000e+00,   6.44894524e-27,  -0.00000000e+00],
           [ -0.00000000e+00,   1.61223638e-27,  -0.00000000e+00],
           [ -0.00000000e+00,   3.65586483e-28,  -0.00000000e+00],
           [ -0.00000000e+00,   6.44894556e-29,   6.67782327e-44],
           [ -0.00000000e+00,   1.61223639e-29,  -0.00000000e+00]])




```python
deltaP_exp = np.array(ZF_i[-1, [Ipx, Ipy, Ipz], :] - ZF_i[0, [Ipx, Ipy, Ipz], :]).T
deltaP_avg = np.array(ZFgc_i[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i[0, [Ipx, Ipy, Ipz], :]).T
deltaP_avg1 = np.array(ZFgc_i_x1[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x1[0, [Ipx, Ipy, Ipz], :]).T
deltaP_avg2 = np.array(ZFgc_i_x2[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x2[0, [Ipx, Ipy, Ipz], :]).T
deltaP_avg5 = np.array(ZFgc_i_x5[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x5[0, [Ipx, Ipy, Ipz], :]).T
deltaP_mgn = magnus_deltaPIon(T_intxn, ZI_ion, ZI_elec_gc)

dDeltaP_mx = deltaP_mgn - deltaP_exp
dDeltaP_ax = deltaP_avg - deltaP_exp
dDeltaP_a1x = deltaP_avg1 - deltaP_exp
dDeltaP_a2x = deltaP_avg2 - deltaP_exp
dDeltaP_a5x = deltaP_avg5 - deltaP_exp

relErr_mgn = (np.linalg.norm(dDeltaP_mx, axis = 1)
             / np.linalg.norm(deltaP_exp, axis = 1))
relErr_avg = (np.linalg.norm(dDeltaP_ax, axis = 1)
             / np.linalg.norm(deltaP_exp, axis = 1))
relErr_avg1 = (np.linalg.norm(dDeltaP_a1x, axis = 1)
              / np.linalg.norm(deltaP_exp, axis = 1))
relErr_avg2 = (np.linalg.norm(dDeltaP_a2x, axis = 1)
              / np.linalg.norm(deltaP_exp, axis = 1))
relErr_avg5 = (np.linalg.norm(dDeltaP_a5x, axis = 1)
              / np.linalg.norm(deltaP_exp, axis = 1))

print "|ΔP_mgn - ΔP_exp| =\n", np.linalg.norm(dDeltaP_mx, axis = 1)
print "|ΔP_avg - ΔP_exp| =\n", np.linalg.norm(dDeltaP_ax, axis = 1)
print "|ΔP_avg1 - ΔP_exp| =\n", np.linalg.norm(dDeltaP_a1x, axis = 1)
print "|ΔP_avg2 - ΔP_exp| =\n", np.linalg.norm(dDeltaP_a2x, axis = 1)
print "|ΔP_avg5 - ΔP_exp| =\n", np.linalg.norm(dDeltaP_a5x, axis = 1)

print 
print "|ΔP_exp| =\n", np.linalg.norm(deltaP_exp, axis = 1)
print 
print "relErr_mgn / % =\n", 100 * relErr_mgn
print "relErr_avg / % =\n", 100 * relErr_avg
print "relErr_avg1 / % =\n", 100 * relErr_avg1
print "relErr_avg2 / % =\n", 100 * relErr_avg2
print "relErr_avg5 / % =\n", 100 * relErr_avg5

print 
print "log_10(relErr_mgn) =\n", np.log10(relErr_mgn)
print "log_10(relErr_avg) =\n", np.log10(relErr_avg)
print "log_10(relErr_avg1) =\n", np.log10(relErr_avg1)
print "log_10(relErr_avg2) =\n", np.log10(relErr_avg2)
print "log_10(relErr_avg5) =\n", np.log10(relErr_avg5)

print 
print "b / μm =\n", impactParameterB / micro
```

    |ΔP_mgn - ΔP_exp| =
    [  1.04511769e-25   2.25584307e-26   7.78830074e-28   2.59961647e-29   1.90782436e-31
       1.54621646e-31   4.08192405e-32]
    |ΔP_avg - ΔP_exp| =
    [  4.17072784e-28   4.27045990e-29   1.33101652e-29   3.98268013e-30   9.29292757e-31
       1.64470326e-31   4.11271376e-32]
    |ΔP_avg1 - ΔP_exp| =
    [  8.29093685e-28   5.88770743e-29   1.36644307e-29   3.99280893e-30   9.28943509e-31
       1.64355196e-31   4.10974159e-32]
    |ΔP_avg2 - ΔP_exp| =
    [  3.25294329e-26   1.05360960e-28   1.47365968e-29   4.02345011e-30   9.27902249e-31
       1.64009960e-31   4.10082720e-32]
    |ΔP_avg5 - ΔP_exp| =
    [  2.64204915e-26   8.60725969e-28   2.07045972e-29   4.31698974e-30   9.33183291e-31
       1.63682154e-31   4.09042036e-32]
    
    |ΔP_exp| =
    [  5.83086619e-26   1.79110313e-26   5.67359853e-27   1.58638483e-27   3.65764411e-28
       6.46440746e-29   1.61631831e-29]
    
    relErr_mgn / % =
    [  1.79238840e+02   1.25947134e+02   1.37272680e+01   1.63870482e+00   5.21599234e-02
       2.39189202e-01   2.52544565e-01]
    relErr_avg / % =
    [ 0.71528444  0.23842624  0.23459829  0.25105385  0.25406867  0.25442444  0.25444949]
    relErr_avg1 / % =
    [ 1.42190484  0.32871962  0.2408424   0.25169233  0.25397318  0.25424634  0.25426561]
    relErr_avg2 / % =
    [ 55.78833721   0.58824619   0.25973986   0.25362384   0.2536885    0.25371229   0.25371408]
    relErr_avg5 / % =
    [ 45.31143511   4.80556341   0.36492884   0.27212752   0.25513234   0.25320519   0.25307022]
    
    log_10(relErr_mgn) =
    [ 0.25343213  0.10018829 -0.86241589 -1.78549927 -3.28266306 -2.62125843 -2.59766197]
    log_10(relErr_avg) =
    [-2.14552122 -2.62264595 -2.62967515 -2.60023312 -2.59504889 -2.59444117 -2.59439841]
    log_10(relErr_avg1) =
    [-1.84712947 -2.48317437 -2.61826706 -2.59913002 -2.59521214 -2.59474528 -2.59471238]
    log_10(relErr_avg2) =
    [-0.25345658 -2.23044088 -2.5854614  -2.59580993 -2.59569922 -2.5956585  -2.59565542]
    log_10(relErr_avg5) =
    [-0.34379218 -1.31825569 -2.43779182 -2.56522753 -2.59323449 -2.59652739 -2.59675895]
    
    b / μm =
    [   1.    2.    5.   10.   21.   50.  100.]



```python
eV_closest_approach = (invFourPiEps0 * Z_ion * qe ** 2 / impactParameterB)
E_kinetic_e = (me / 2) * (V_perp ** 2 + V_z ** 2)
print "E_pot_ca =\n", eV_closest_approach
print "E_kin =\n", E_kinetic_e
print 
print "E_pot_ca / E_kin =\n", (eV_closest_approach / E_kinetic_e)
print 
print "log_10(E_pot_ca / E_kin) =\n", np.log10(eV_closest_approach / E_kinetic_e)
```

    E_pot_ca =
    [  2.30707751e-22   1.15353876e-22   4.61415503e-23   2.30707751e-23   1.09860834e-23
       4.61415503e-24   2.30707751e-24]
    E_kin =
    [  3.87932171e-23   1.55172869e-22   9.69830428e-22   3.87932171e-21   1.71078088e-20
       9.69830428e-20   3.87932171e-19]
    
    E_pot_ca / E_kin =
    [  5.94711572e+00   7.43389465e-01   4.75769257e-02   5.94711572e-03   6.42167770e-04
       4.75769257e-05   5.94711572e-06]
    
    log_10(E_pot_ca / E_kin) =
    [ 0.77430639 -0.1287836  -1.32260362 -2.22569361 -3.1923515  -4.32260362 -5.22569361]



```python
mae_compare = plt.figure()

# my_font_12 = {'fontname' : 'Arial', 'size' : '12'}
# my_font_14 = {'fontname' : 'Arial', 'size' : '14'}
# plt.rc('font', **my_font_12)
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_mgn, 'rs')
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_avg, 'bo')
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_avg1, 'gd')
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_avg2, 'rd')
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_avg5, 'yd')

plt.xlabel("$eV_{\mathrm{ca}} / E_{\mathrm{kin}}$")
plt.ylabel("relative error in $\Delta{P}_{\mathrm{ion}}$")

#plt.annotate(s="bar",xy=(0.001,1))
#plt.text(0.1,0.01, "foo", rotation = 45)

plt.show()
```


```python
mae_compare.savefig("/Users/dabell/RadiaSoft/MCool/MgnAvgExpCompare.pdf")
```


```python
# will have to fix sign issues
deltaP_exp = np.array(ZF_i[-1, [Ipx, Ipy, Ipz], :] - ZF_i[0, [Ipx, Ipy, Ipz], :]).T
deltaP_avg = (ZFgc_i_x1[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i_x1[0, [Ipx, Ipy, Ipz], :]).T
deltaP_mgn = magnus_deltaPIon(T_intxn, ZI_ion, ZI_elec_gc)

dDeltaP_ax = deltaP_avg - deltaP_exp
dDeltaP_mx = deltaP_mgn - deltaP_exp

relErr_avg = (np.linalg.norm(dDeltaP_ax, axis = 1)
             / np.linalg.norm(deltaP_exp, axis = 1))
relErr_mgn = (np.linalg.norm(dDeltaP_mx, axis = 1)
             / np.linalg.norm(deltaP_exp, axis = 1))

print "|ΔP_avg - ΔP_exp| =\n", np.linalg.norm(dDeltaP_ax, axis = 1)
print "|ΔP_mgn - ΔP_exp| =\n", np.linalg.norm(dDeltaP_mx, axis = 1)
print 
print "|ΔP_exp| =\n", np.linalg.norm(deltaP_avg, axis = 1)
print 
print "relErr_avg / % =\n", 100 * relErr_avg
print "relErr_mgn / % =\n", 100 * relErr_mgn
print 
print "log_10(relErr_avg) =\n", np.log10(relErr_avg)
print "log_10(relErr_mgn) =\n", np.log10(relErr_mgn)
print 

print "b / μm =\n", impactParameterB / micro
```

    |ΔP_avg - ΔP_exp| =
    [  8.29093685e-28   5.88770743e-29   1.36644307e-29   3.99280893e-30   1.02364248e-30
       1.64355196e-31   4.10974159e-32]
    |ΔP_mgn - ΔP_exp| =
    [  1.04511769e-25   2.25584307e-26   7.78830074e-28   2.59961647e-29   1.09581428e-31
       1.54621646e-31   4.08192405e-32]
    
    |ΔP_exp| =
    [  5.77351905e-26   1.78523936e-26   5.65993463e-27   1.58239202e-27   4.02100997e-28
       6.44797194e-29   1.61220857e-29]
    
    relErr_avg / % =
    [ 1.42190484  0.32871962  0.2408424   0.25169233  0.25392704  0.25424634  0.25426561]
    relErr_mgn / % =
    [  1.79238840e+02   1.25947134e+02   1.37272680e+01   1.63870482e+00   2.71830143e-02
       2.39189202e-01   2.52544565e-01]
    
    log_10(relErr_avg) =
    [-1.84712947 -2.48317437 -2.61826706 -2.59913002 -2.59529104 -2.59474528 -2.59471238]
    log_10(relErr_mgn) =
    [ 0.25343213  0.10018829 -0.86241589 -1.78549927 -3.56570239 -2.62125843 -2.59766197]
    
    b / μm =
    [   1.    2.    5.   10.   20.   50.  100.]



```python
eV_closest_approach = (invFourPiEps0 * Z_ion * qe ** 2 / impactParameterB)
E_kinetic_e = (me / 2) * (V_perp ** 2 + V_z ** 2)
print "E_pot_ca =\n", eV_closest_approach
print "E_kin =\n", E_kinetic_e
print 
print "E_pot_ca / E_kin =\n", (eV_closest_approach / E_kinetic_e)
print 
print "log_10(E_pot_ca / E_kin) =\n", np.log10(eV_closest_approach / E_kinetic_e)
```

    E_pot_ca =
    [  2.30707751e-22   1.15353876e-22   4.61415503e-23   2.30707751e-23   1.15353876e-23
       4.61415503e-24   2.30707751e-24]
    E_kin =
    [  3.87932171e-23   1.55172869e-22   9.69830428e-22   3.87932171e-21   1.55172869e-20
       9.69830428e-20   3.87932171e-19]
    
    E_pot_ca / E_kin =
    [  5.94711572e+00   7.43389465e-01   4.75769257e-02   5.94711572e-03   7.43389465e-04
       4.75769257e-05   5.94711572e-06]
    
    log_10(E_pot_ca / E_kin) =
    [ 0.77430639 -0.1287836  -1.32260362 -2.22569361 -3.1287836  -4.32260362 -5.22569361]



```python
llres = plt.figure()
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_avg, 'bo')
plt.loglog(eV_closest_approach / E_kinetic_e, relErr_mgn, 'rs')
plt.xlabel("$eV_{\mathrm{cl.appr.}}$")
plt.ylabel("$eV_{\mathrm{cl.appr.}}$")
plt.annotate(s="Averging",xy=(1.e-1,1.e-3))
plt.annotate(s="Magnus",xy=(1.e-2,5.e-1))
plt.show()

llres.savefig("/Users/dabell/Desktop/foo.pdf")
```

# Explore other values of the perturbation parameters


```python
Z_ion = 1
M_ion = mp
e_temp = 300. # Kelvin
B_mag = 1. # Tesla

N_gyro = 100    # enforce adequate averaging
tan_alpha = 5.0 # and an adequate opening angle

""" derived quantities """
Omega_e = omega_Larmor(me, B_mag, Z = -1)
T_e = (2 * pi) / abs(Omega_e)
```


```python
def minImpactParam(magB, EkOverEV, bOverRL):
    numer = 2 * EkOverEV * (invFourPiEps0 * Z_ion * me)
    denom = ((1 / bOverRL) ** 2 + (tan_alpha / (N_gyro * pi)) ** 2) * magB ** 2
    return (numer / denom) ** (1/3)
```


```python
bOverLarmorR = 20.0 # 20 ==> max perturbation ratio of about 0.1
EkinOverVmax = 10.0 # 10 ==> eV_closest_approach / E_kinetic = 0.1
print minImpactParam(B_mag, EkinOverVmax, bOverLarmorR) / micro
```

    3.90333173895



```python
decades = 2
pts_per_decade = 6
logBs = np.linspace(0., 1. * float(decades), num = 1 + pts_per_decade * decades, endpoint = True)
bvals = np.power(10, logBs)
impactParameterB = micro * bvals

print "b / μm =    ", impactParameterB / micro
print "b_min / μm =", minImpactParam(B_mag, EkinOverVmax, bOverLarmorR) / micro
```

    b / μm =     [   1.            1.46779927    2.15443469    3.16227766    4.64158883    6.81292069   10.
       14.67799268   21.5443469    31.6227766    46.41588834   68.12920691  100.        ]
    b_min / μm = 3.90333173895



```python
""" (the following depend on the impact parameter) """
LarmorR = impactParameterB / bOverLarmorR  # (this version is defined positive)
V_perp = - LarmorR * Omega_e
L_intxn = tan_alpha * impactParameterB #[-1] * np.ones_like(impactParameterB)
V_z = 2 * L_intxn / (N_gyro * T_e)
T_intxn = 2 * L_intxn / V_z
```


```python
num_steps_per_gyro = 40
delta_Te = T_e / num_steps_per_gyro
print T_intxn / delta_Te
num_steps = int(np.floor(T_intxn[0] / delta_Te))
print "num_steps =", num_steps
```

    [ 4000.  4000.  4000.  4000.  4000.  4000.  4000.  4000.  4000.  4000.  4000.  4000.  4000.]
    num_steps = 4000



```python
""" initial condition for the ion --- just one, for now """
QVion = np.array([
          ( 0.0 * rho_gc,  0.000 * V_th, 0.0 * rho_gc,  0.000 * V_th, 0.0,  0.000 * V_th)
        ]).transpose()
Zion = QVtoZion.dot(QVion)

""" initial conditions for the electrons """
QVelec = np.asarray([
           np.zeros(impactParameterB.shape),
           V_perp,
           impactParameterB - LarmorR,
           np.zeros(impactParameterB.shape),
           - L_intxn,
           L_intxn * abs(Omega_e) / (pi * N_gyro)
        ])
Zelec = QVtoZe.dot(QVelec)

num_elec0 = Zelec.shape[1]
num_ion0 = Zion.shape[1]

"""
=== create arrays of electrons and ions ===
Here we arrange them so that we can pair each ion with each
electron, and compute the \Delta{p} for each interaction.
"""
ZI_elec = np.hstack([Zelec for e in range(num_ion0)])
ZI_ion = Zion[:, np.arange(num_ion0 * num_elec0) // num_elec0]

num_elec = ZI_elec.shape[1]
num_ion = ZI_ion.shape[1]
```


```python
ZF_i, ZF_e = apply_MapZ(delta_Te, num_steps, ZI_ion, ZI_elec)
```


```python
navg = 1 # number of gyrotron oscillations over which to average; set to 1, 2, 4, 5, or 10
ZI_elec_gc = toGuidingCenter(ZI_elec)
ZFgc_i, ZFgc_e = apply_MapZgc(navg * 40*delta_Te, num_steps//40 // navg, ZI_ion, ZI_elec_gc)
```


```python
deltaP_exp = np.array(ZF_i[-1, [Ipx, Ipy, Ipz], :] - ZF_i[0, [Ipx, Ipy, Ipz], :]).T
deltaP_avg = np.array(ZFgc_i[-1, [Ipx, Ipy, Ipz], :] - ZFgc_i[0, [Ipx, Ipy, Ipz], :]).T
deltaP_mgn = magnus_deltaPIon(T_intxn, ZI_ion, ZI_elec_gc)

print deltaP_exp
print
print deltaP_avg
print
print deltaP_mgn
```

    [[  1.08030578e-26   5.72808106e-26   1.45007842e-27]
     [  1.11981759e-27   2.29054810e-26   4.74338761e-27]
     [  4.03458322e-28   1.66767979e-26   1.38486416e-27]
     [  1.21157181e-28   1.10760987e-26   3.98681026e-28]
     [  2.76032939e-29   6.40599505e-27   9.05956974e-29]
     [  4.96244634e-30   3.29085593e-27   1.63717477e-29]
     [  7.85620796e-31   1.58638251e-27   2.60010730e-30]
     [  1.18279915e-31   7.45748787e-28   3.91944864e-31]
     [  1.75039675e-32   3.47567664e-28   5.80271684e-32]
     [  2.57594045e-33   1.61537214e-28   8.54061442e-33]
     [  3.78408946e-34   7.50099341e-29   1.25468090e-33]
     [  5.55573638e-35   3.48210818e-29   1.84212636e-34]
     [  8.15538018e-36   1.61631831e-29   2.70410697e-35]]
    
    [[  1.05836708e-26   5.67504416e-26   8.51768805e-28]
     [  1.10752456e-27   2.27810318e-26   4.74122323e-27]
     [  4.00940540e-28   1.66262914e-26   1.38431451e-27]
     [  1.20545450e-28   1.10493620e-26   3.98478768e-28]
     [  2.74641067e-29   6.39066436e-27   9.05285901e-29]
     [  4.93653338e-30   3.28272166e-27   1.63565389e-29]
     [  7.81439841e-31   1.58238970e-27   2.59745819e-30]
     [  1.17645959e-31   7.43858940e-28   3.91532568e-31]
     [  1.74099270e-32   3.46684860e-28   5.79654896e-32]
     [  2.56209059e-33   1.61126617e-28   8.53150612e-33]
     [  3.76373889e-34   7.48192282e-29   1.25334141e-33]
     [  5.52585575e-35   3.47325458e-29   1.84015906e-34]
     [  8.11151685e-36   1.61220857e-29   2.70121883e-35]]
    
    [[ -0.00000000e+00   1.61222627e-25  -1.01894869e-40]
     [ -0.00000000e+00   7.48332357e-26   8.78107061e-41]
     [ -0.00000000e+00   3.47345583e-26  -3.78365367e-41]
     [ -0.00000000e+00   1.61223607e-26  -0.00000000e+00]
     [ -0.00000000e+00   7.48333795e-27  -0.00000000e+00]
     [ -0.00000000e+00   3.47345794e-27  -3.02692474e-42]
     [ -0.00000000e+00   1.61223638e-27  -0.00000000e+00]
     [ -0.00000000e+00   7.48333841e-28   5.61989612e-43]
     [ -0.00000000e+00   3.47345800e-28  -0.00000000e+00]
     [ -0.00000000e+00   1.61223639e-28  -0.00000000e+00]
     [ -0.00000000e+00   7.48333842e-29  -4.49591691e-44]
     [ -0.00000000e+00   3.47345801e-29  -1.93723187e-44]
     [ -0.00000000e+00   1.61223639e-29  -0.00000000e+00]]



```python
dDeltaP_ax = deltaP_avg - deltaP_exp
dDeltaP_mx = deltaP_mgn - deltaP_exp

relErr_avg = (np.linalg.norm(dDeltaP_ax, axis = 1)
             / np.linalg.norm(deltaP_exp, axis = 1))
relErr_mgn = (np.linalg.norm(dDeltaP_mx, axis = 1)
             / np.linalg.norm(deltaP_exp, axis = 1))

eV_closest_approach = (invFourPiEps0 * Z_ion * qe ** 2 / impactParameterB)
E_kinetic_e = (me / 2) * (V_perp ** 2 + V_z ** 2)
eVcaOverEkin = eV_closest_approach / E_kinetic_e

llres = plt.figure()
plt.loglog(eVcaOverEkin, relErr_avg, 'bo')
plt.loglog(eVcaOverEkin, relErr_mgn, 'rs')
plt.annotate(s="Averging",xy=(1.e-1,1.e-3))
plt.annotate(s="Magnus",xy=(1.e-2,5.e-1))

plt.xlabel("$eV_{\mathrm{ca}} / E_{\mathrm{kin}}$")
plt.ylabel("relative error in $\Delta{P}_{\mathrm{ion}}$")

plt.show()

llres.savefig("/Users/dabell/Desktop/foo.pdf")
```


```python
num_steps_per_gyro = 80
delta_Te = T_e / num_steps_per_gyro
num_steps = int(np.floor(T_intxn / delta_Te))
print "num_steps =", num_steps
```


```python
ZF0_i, ZF0_e = apply_MapZ_0(delta_Te, num_steps, ZI_ion, ZI_elec)
```


```python
ZF_i, ZF_e = apply_MapZ(delta_Te, num_steps, ZI_ion, ZI_elec)
```


```python
ZF0gc_e = np.asarray(map(toGuidingCenter, ZF0_e))
ZFgc_e = np.asarray(map(toGuidingCenter, ZF_e))
```


```python
nsteps = ZF0_e.shape[0]
```


```python
H2b_0 = np.asarray(map(lambda zi, ze: H_twobody_0(zi, ze), ZF0_i, ZF0_e))
H2b   = np.asarray(map(lambda zi, ze: H_twobody(zi, ze), ZF_i, ZF_e))

m2r2_0 = np.asarray(map(lambda zi, ze: minus2RLsq(zi, ze), ZF0_i, ZF0gc_e))
p2r2_0 = np.asarray(map(lambda zi, ze: plus2RLsq(zi, ze), ZF0_i, ZF0gc_e))

m2r2 = np.asarray(map(lambda zi, ze: minus2RLsq(zi, ze), ZF_i, ZFgc_e))
p2r2 = np.asarray(map(lambda zi, ze: plus2RLsq(zi, ze), ZF_i, ZFgc_e))
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b_0[:, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b_0[:, ptcl] - H2b_0[0, ptcl])
```


```python
plt.figure()
for ptcl in range(0, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b_0[:, ptcl] - H2b_0[0, ptcl])
```


```python
plt.figure()
for ptcl in range(1, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b_0[:, ptcl] - H2b_0[0, ptcl])
```


```python
plt.figure()
for ptcl in range(2, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b_0[:, ptcl] - H2b_0[0, ptcl])
```


```python
plt.figure()
for ptcl in range(3, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b_0[:, ptcl] - H2b_0[0, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(0, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(1, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(2, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(3, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             H2b[:, ptcl] - H2b[0, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             ZF0gc_e[:, Ipx, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             ZF0gc_e[:, Ipx, ptcl] - ZF0gc_e[0, Ipx, ptcl])
```


```python
plt.figure()
for ptcl in range(0, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             ZF0gc_e[:, Ipx, ptcl] - ZF0gc_e[0, Ipx, ptcl])
```


```python
plt.figure()
for ptcl in range(1, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             ZF0gc_e[:, Ipx, ptcl] - ZF0gc_e[0, Ipx, ptcl])
```


```python
plt.figure()
for ptcl in range(2, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             ZF0gc_e[:, Ipx, ptcl] - ZF0gc_e[0, Ipx, ptcl])
```


```python
plt.figure()
for ptcl in range(3, 16, 4):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             ZF0gc_e[:, Ipx, ptcl] - ZF0gc_e[0, Ipx, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             p2r2_0[:, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             p2r2[:, ptcl] - p2r2_0[:, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             m2r2[:, ptcl] - m2r2_0[:, ptcl])
```


```python



plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             m2r2_0[:, ptcl])
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             m2r2[:, ptcl])
```


```python
plt.figure()
ptcl = 0
plt.plot(delta_Te * np.asarray(range(nsteps)),
         m2r2[:, ptcl])
plt.plot(delta_Te * np.asarray(range(nsteps)),
         p2r2[:, ptcl])
```


```python
plt.figure()
ptcl = 4
plt.plot(delta_Te * np.asarray(range(nsteps)),
         m2r2[:, ptcl])
plt.plot(delta_Te * np.asarray(range(nsteps)),
         p2r2[:, ptcl])
```


```python
plt.figure()
ptcl = 7
plt.plot(delta_Te * np.asarray(range(nsteps)),
         m2r2[:, ptcl])
plt.plot(delta_Te * np.asarray(range(nsteps)),
         p2r2[:, ptcl])
```


```python
ZI_elec_gc = toGuidingCenter(ZI_elec)
ZFgc_i, ZFgc_e = apply_MapZgc(delta_Te, num_steps, ZI_ion, ZI_elec_gc)
```


```python
(ZFgc_i[-1, [Ipx, Ipy], :] - ZFgc_i[0, [Ipx, Ipy], :]).T
```


```python
(ZFgc_i[-1, [Ipx, Ipy], :] - ZFgc_i[0, [Ipx, Ipy], :]).T
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(ZFgc_e[:, Iy,  ptcl] / micro,
             ZFgc_e[:, Ipy, ptcl] / me / mega)
```


```python
nsteps = ZFgc_i.shape[0]
nsteps
```


```python
hgc = []
for i in range(nsteps):
    hgc.append(H_gc(ZFgc_i[i], ZFgc_e[i]))
hgc = np.asarray(hgc)
```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             hgc[:, ptcl])
```


```python
plt.figure()
ptcl = 6
plt.plot(delta_Te * np.asarray(range(nsteps)),
         hgc[:, ptcl])
```


```python
plt.figure()
ptcl = 6
plt.plot(delta_Te * np.asarray(range(nsteps)),
         hgc[:, ptcl])
```


```python

```


```python

```


```python
plt.figure()
for ptcl in range(num_elec):
    plt.plot(delta_Te * np.asarray(range(nsteps)),
             hgc[:, ptcl])
```


```python

```


```python

```


```python

```


```python

```

bfaaf9af6dfc1162341a220c7ad6a96802a05d60
ipynb
Jupyter Notebook
examples/notebooks/mcool.ipynb
radiasoft/rsdynfric
0d905e1ef85063ec6fe93df64751f7b8c19da160
Apache-2.0
examples/notebooks/mcool.ipynb
radiasoft/rsdynfric
0d905e1ef85063ec6fe93df64751f7b8c19da160
Apache-2.0
2022-03-23T22:53:10.000Z
2022-03-28T19:02:40.000Z
examples/notebooks/mcool.ipynb
radiasoft/rsdynfric
0d905e1ef85063ec6fe93df64751f7b8c19da160
Apache-2.0
2017-01-06T04:24:29.000Z
2020-02-16T01:17:20.000Z
Qwen/Qwen-72B
1. YES
2. YES