# particle parameters			# Carbon 4+
m0 	= 12 * ma;			# mass [kg]
q0 	= -4 * qe;			# charge [C]
#Eu	= 8e3;				# beam energy is 8 keV/u from ion source [eV] (set in batch file)
E0	= 12 * Eu;			# beam energy is u * keV/u
gamma	= 1-qe*E0/(m0*c*c);  		# corresponding lorentz factor gamma [] 
beta 	= sqrt(1-gamma^-2);  		# corresponding normalized velocity beta []
BG   	= sqrt(gamma^2 - 1);		# corresponding beta-gamma factor []

# simulation parameters 
#npart  = 100;    			# number of particles [] (set in batch file)
zlen	= rfqfac * 4.1;			# length of simulation [m]
pipe	= 0.0036;			# beam pipe 3.6mm radius [m]
sigcut	= 1.8;				# cut gaussian distributions after 1.8 sigma []
rmsratio= sigcut;			# ratio of total parameters to rms parameters []

# beam parameters 
I    	= 5e-4;				# beam current is 0.5 mA [A]
tbeam	= 3e-9;				# beam output pulse is 3 ns [s]
zbunch	= tbeam*(beta*c);		# bunch length [m]
Qtot	= I*tbeam;			# total charge is current x time [C]
dE   	= 0.2;				# energy spread [fraction of E0]
normemit= 0.025e-6;			# normalised emittance is 0.25 pi mm mrad [pi m rad]
rmsemit	= normemit/BG;			# rms emittance
totrad	= 2.2e-3;			# total radius of beam [m]
rmsrad	= totrad/rmsratio;		# rms radius of beam [m]
maxdiv	= -60e-3;			# total radial divergence (-ve for convergence) [rad]
rmsdiv	= maxdiv/rmsratio;		# rms radial divergence [rad]
divx 	= rmsdiv/totrad;		# x divergence [rad]
divy 	= rmsdiv/totrad;		# y divergence [rad]

# start bunch 
setparticles("beam",npart,m0,q0,Qtot);	# particles of mass m0 and charge q0, total charge Qtot, npart macroparticles
setzdist("beam", "u", 
		-zbunch/2, zbunch/2);	# uniform z distribution

setrxydist("beam", "g", 0, rmsrad, 
		0, sigcut);		# gaussian radial distribution
setphidist("beam","u", 0, 2*pi); 	# uniform phi distribution
setGBrxydist("beam", "g", 0, 
		rmsdiv*BG, 0, sigcut);	# gaussian radial velocity distribution
setGBphidist("beam", "u", 0, 2*pi);	# uniform phi velocity distribution
setGBxemittance("beam", rmsemit*BG);	# x emittance
setGByemittance("beam", rmsemit*BG);	# y emittance
setGBzdist("beam", "u", BG, 0);		# no distribution in z velocity
setGdist("beam", "u", gamma, 0);	# no distribution in total velocity
addxdiv("beam", 0, divx*BG);		# x divergence
addydiv("beam", 0, divy*BG);		# y divergence

# space-charge
#spacecharge3D();			# full point-to-point space charge calculation
#setcharge2Dcircle("beam", QTot);	# set charge for spacecharge2Dcircle calculation
#spacecharge2Dcircle();			# cylindrically symmetry space charge calculation
#SCtree2D(npart);			# 2D bunch model
#SCtree3D(npart);			# 3D bunch model

# rfq data
#ffac	= 1.0;				# field factor (increase/decrease field without modifying input file) []  (set in batch file)
phi	= -180 * pi / 180;		# phase of field [rad]
#freq	= 200.0e6;			# frequency of field [Hz] (set in batch file)
w	= freq * 2 * pi;		# frequency of field [Hz rad]
rfqfront= -0.02;			# start rfq before particles to avoid fringe fields [m]
#rfqfac	= 1;				# factor to multiply length of rfq (set in batch file)
rfqlen	= rfqfac*4.05968;		# rfq length [m]
rfqmid	= (rfqlen/2) + rfqfront;	# rfq middle [m]
rfqend	= rfqlen + rfqfront;		# rfq end [m]
map3D_EB				# read in electromagnetic field map from file
 ("wcs", 				#  - location specified in world coordinate system
   0, 0, rfqfront, 			#  - location of start of rfq in wcs
   1, 1, 0, 				#  - direction of field x-axis: x+y (i.e. 45 degrees)
  -1, 1, 0, 				#  - direction of field y-awis: y-x (i.e. 45 degrees)
  "fields.gdf", 			#  - file to read
  "x", "y", "z", 			#  - first three columns are field coordinates
  "Ex", "Ey", "Ez", 			#  - next three columns are E-field
  "Bx", "By", "Bz", 			#  - next three columns are B-field
   ffac, phi, w);			#  - field parameters
rmax("wcs", "z", rfqmid, pipe, rfqlen);	# constrain simulation to beam pipe

# output
tout(0, 8e-7, 4e-9);			# output at fixed time steps
#screen("wcs","I",rfqfront+rfqlen);	# output at a screen at the end of the rfq
screen("wcs","I",0,rfqfront+rfqlen,0.1);# output at a screen every 10cm