# solving electron motion in undulator by Boris method

I am trying to use Boris scheme to solve the electron trajectory in undulator.

The undulator field I used is: $$B_x = b_0\sin(2\pi \tfrac{z}{\lambda_u})$$ where $$b_0 = \dfrac{2\pi c_{0}K}{q m_{e} \lambda_{u}}$$, $$K$$ is undulator strength and $$\lambda_{u}$$ is undulator period.

I test my code with the parameters:

$$v_{0}= (0,0,0.9c_{0})$$

$$K=1.0e^{-5}$$

$$\lambda_{u}=0.01$$

It seems there is a drift effect in the electron trajectory, while the trajectory predicted from theory should be sinusoidal. Does anyone has idea what is the problem?

Bellows is my julia code:

using LinearAlgebra

#### physics const
const me = 9.11e-31
const c0 = 299792458
const h_bar = 1.0545718e-31
const q = -1.6e-19

#### user define
const nstep = 25000
const dt = 1.0e-14

#undulator paramter#
const K = 0.00001
const lambda_u = 0.01
const b0 = 2*pi*(me/1.6e-19)*c0*K/(lambda_u)

function efield(x::Vector{Float64}, t::Float64)
return [0.0;0.0;0.0]
end
function bfield(x::Vector{Float64}, t::Float64)
return [b0;0.0;0.0]*sin(2*pi*x/lambda_u)
end

mutable struct Particle
position::Vector{Float64}
beta::Vector{Float64}
end

function boris_push(p::Particle, t::Float64, dt::Float64)
x = p.position
beta = p.beta

x_half = x + 0.5*c0*dt*beta

u = beta/sqrt(1.0-norm(beta)^2)
u_minus = u + 0.5*q*dt/(me*c0)*efield(x_half, t+0.5*dt)
t_vec = 0.5*dt*q*bfield(x_half, t+0.5*dt) / ( me*sqrt( 1.0 + norm(u_minus)^2 ) )
s_vec = 2.0*t_vec / ( 1.0 + dot(t_vec,t_vec) )
u_plus = u_minus + cross(u_minus + (cross(u_minus,t_vec)), s_vec)
u_final = u_plus + 0.5*q*dt/(me*c0)*efield(x_half, t+0.5*dt)

beta_final = u_final/sqrt( 1+norm(u_final)^2 )
x_final = x_half + 0.5*dt*c0*beta_final

p.position = x_final
p.beta = beta_final
end

p = Particle([0.0;0.0;0.0],[0.0;0.0;0.9])
for i=1:nstep
t = i*dt
boris_push(p, t, dt)
x = p.position
beta = p.beta
println(x, " ",x," ",x, " ", beta, " ", beta, " ", beta)
end