Test Particle Simulation

Equation of motion of particles can be written as follows \begin{gathered} m\frac{{d{\mathbf{v}}}}{{dt}} = q\left( {{\mathbf{E}} + \frac{{\mathbf{v}}}{c} \times {\mathbf{B}}} \right), \hfill \\ \frac{{d{\mathbf{r}}}}{{dt}} = {\mathbf{v}} \hfill \\ \end{gathered} The equations are solved by the Bunnemann-Boriss method.
I show a particle motion in X-type magnetic fields.

%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import scipy.fftpack as sf
import matplotlib.animation as animation
from tqdm import tqdm
import matplotlib as mpl
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
from numba import jit,f8,u8

def testp(nop,nt,dt,isav):
    global qomdt
    
    q=-1; m=1;
    qomdt=dt*q/2/m;
    up=np.zeros((3,nop))
    vp=np.random.randn(3,nop)*0.1
    gm=1/np.sqrt(1+vp[0,:]**2+vp[1,:]**2+vp[2,:]**2)
    bf=np.zeros((3,nop))
    ef=np.zeros((3,nop))
    
    uphst=np.zeros((nt,3))
    
    for it in range(nt):
        ef[0,:]=0       ;ef[1,:]=0       ;ef[2,:]=0;
        bf[0,:]=-up[1,:];bf[1,:]=-up[0,:];bf[2,:]=0
        
        vp=acc(nop,up,vp,gm,bf,ef)
        up=up+vp*dt
   
         
        if it%isav==0:
            uphst[it//isav,:]=up[:,0]

    return locals()

@jit('f8[:,:](u8,f8[:,:],f8[:,:],f8[:],f8[:,:],f8[:,:])')
def acc(nop,up,vp,gm,bf,ef):
    for ip in range(nop):

        exap=qomdt*ef[0,ip]
        eyap=qomdt*ef[1,ip]
        ezap=qomdt*ef[2,ip]
        bxap=qomdt*bf[0,ip]
        byap=qomdt*bf[1,ip]
        bzap=qomdt*bf[2,ip]
        #---half acceleration
        gvxs=vp[0,ip]*gm[ip]+exap
        gvys=vp[1,ip]*gm[ip]+eyap
        gvzs=vp[2,ip]*gm[ip]+ezap
        #---first half rotation
        gmm=mt.sqrt(1+gvxs**2+gvys**2+gvzs**2)
        #enesave[it,5+isp]=enesave[it,5+isp]+gm[isp,ip]-1
        bxap=bxap/gmm
        byap=byap/gmm
        bzap=bzap/gmm
        vvx=gvxs+gvys*bzap-gvzs*byap
        vvy=gvys+gvzs*bxap-gvxs*bzap
        vvz=gvzs+gvxs*byap-gvys*bxap
        #---second half rotation + second half acceleration
        f=2/(1+bxap**2+byap**2+bzap**2)
        bxap=f*bxap
        byap=f*byap
        bzap=f*bzap
        gvxs=gvxs+vvy*bzap-vvz*byap+exap
        gvys=gvys+vvz*bxap-vvx*bzap+eyap
        gvzs=gvzs+vvx*byap-vvy*bxap+ezap
        #---half acceleration
        gm[ip]  =np.sqrt(1+gvxs**2+gvys**2+gvzs**2)
        vp[0,ip]=gvxs/gm[ip]
        vp[1,ip]=gvys/gm[ip]
        vp[2,ip]=gvzs/gm[ip]
    
    return vp

%%time
nop=1; 
nt=20000; dt=0.1
isav=100
data=testp(nop,nt,dt,isav)
locals().update(data)

CPU times: user 188 ms, sys: 31.2 ms, total: 219 ms
Wall time: 192 ms

### make animation ###
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
def update(it):
    ax.cla()

    ax.plot(uphst[0:it,0],uphst[0:it,1])
    ax.plot(uphst[it,0]  ,uphst[it,1],'ro')
    n = 10
    X, Y = np.mgrid[-n:n, -n:n]
    U, V = -Y, -X
    ax.quiver(X, Y, U, V, 1, alpha=1.)
    ax.set_xlim(-10,10)
    ax.set_ylim(-10,10)

plt.rcParams['animation.html'] = 'html5'
ani = animation.FuncAnimation(fig,update,interval = 100,frames = nt//isav)
plt.close()

ani