# Define the beam parameters
gamma = 8;
radius = 1e-3;
length = 1e-6;
number_of_particles = 1000;
#total_charge of macroparicles
tc = 1e-13
initial_emittance = 1e-10;

#Radius and Current of Current Loop
R=.01;
I=50000;
#Divisor used to scale timestep and lattice
DVR = 10;

# Start the beam
setparticles("beam", number_of_particles, me, qe, tc);
#setrxydist("beam", "u", radius/2, radius);
#setphidist("beam", "u", 0, 2*pi);
setellipse("beam",radius,radius,length);
setGdist("beam", "u", gamma, 0);
#setzdist( "beam", "g", 0, length, 3, 3) ;

# Set its initial emittance to 1nm
setGBxdist("beam", "g", 0.0, 1e-3, 3.0, 3.0);
setGBxemittance("beam", initial_emittance);
setGBydist("beam", "g", 0.0, 1e-3, 3.0, 3.0);
setGByemittance("beam", initial_emittance);

#Current loop halfway through the beamline with radius R and Current I
solenoid("wcs","z",1.5/DVR, R, I);

#Add space charge
Alpha  = 1.2;
Fn     = 0.25;
verror = 0.005;
Nstd   = 8;
spacecharge3Dmesh_caching("Cathode", "MeshNfac", Alpha, "MeshAdapt", Fn, "SolverAcc", verror, "MeshBoxSize", Nstd);

# Compute the emittance components
emittance_component("4d space charge","false");
emittance_component("4d beamline element","false");
emittance_component("4d non-force","false");

# Set the solver tolerance
dtmax = 2.5*1e-13/DVR;
accuracy(6);

# Print touts
tout(0.0, 1.00e-8/DVR, 2e-10/DVR);