Mathematica Code for Integration of the ODEs for an Electrostatic Spherical Deflector in Laboratory Coordinates
November 2017
Eremey Valetov

Supplementary material for :
E. Valetov, M. Berz, and K. Makino,
Validation of Transfer Map Calculation for Electrostatic Deflectors in the Code COSY INFINITY,
Int. Journal of Modern Physics A (2019)

(* Integration of the ODEs of motion for a charged particle of kinetic energy 1 MeV,
mass 1 amu, and charge 1 e in an electrostatic spherical deflector with the inner
sphere charge such that the reference particle would have a circular orbit of radius 1 m. *)
v00=13891388.79714028;
r00=1;
x=0.3; (* Initial x beamline coordinate *)
a=0; (* Initial a beamline coordinate *)
r0={x+r00,0,0}; (* Initial radius vector *)
v0={a v00,v00 Sqrt[1-(a^2)],0}; (* Initial velocity vector *)
mu=r00  (v00^2); (* Parameter μ=αm in the paper *)
vr0=v0[[1]]; 
rr0=r0[[1]];
omega0=v0[[2]]/rr0; (* Initial angular velocity *)
theta0=0; (* Initial polar angle *)
system={ (* ODE's and the initial condition *)
vr'[t]==-mu/(omega[t](((rr[t]))^2))+rr[t] omega[t],
rr'[t]==vr[t]/omega[t],
omega'[t]==-2rr0 v0[[2]]vr[t]/(omega[t](((rr[t]))^3)),
theta'[t]==1,
vr[0]==vr0,
rr[0]==rr0,
omega[0]==omega0,
theta[0]==theta0
};

sol=NDSolve[system,{vr,rr,omega,theta},{t,0,2π}] ;(* Numerical solution of the ODEs *)

rsol[t_]=First[Evaluate[rr[t]/.sol]]; (* Radius magnitude as a function of polar angle theta *)
vsol[t_]=First[Evaluate[vr[t]/.sol]]; (* Velocity magnitude as a function of polar angle theta *)

Print[PolarPlot[{rsol[t],1},{t,0,2π}, LabelStyle->(FontSize->18),AxesLabel -> {"(Subscript[(OverVector[i]), e])","(Subscript[(OverVector[i]), p])"},ImageSize->Large,PlotStyle->{{ColorData[97,"ColorList"][[1]]},{Red,Thin}},Epilog->{FontSize->18,Text[0,Offset[{-8.4,-13.6},{0,0}]]}]]; 
(* Plot one turn of the orbit *)

x1[t_]=rsol[t]-r00; (* x coordinate at polar angle theta *)
a1[t_]=vsol[t]/v00; (* a coordinate at polar angle theta *)

Print["x_f=",x1[Pi/4]]; 
Print["a_f=",a1[Pi/4]];