#**************************************************************************
#**************************************************************************
#          Electron Diffraction Setup - 10/14/2014 C. Gulliford 
#--------------------------------------------------------------------------
# Description: this set-up allows varying of the solenoids + buncher, and
# iris diameter at the center of the first solenoid.  The sample screen is 
# located 1 mm after a fixed size iris (representing the sample crystal)
#**************************************************************************
#**************************************************************************

GPTLICENSE=1251405651;

#--------------------------------------------------------------------------
# MACHINE SETTINGS
#--------------------------------------------------------------------------
gun_voltage=225;
sol_1_current=2.0;
sol_2_current=1.0;

# BUNCHER CAVITY PARAMETERS
buncher_voltage=10;
phiCTB01_0=0.0;
phiCTB01_off=-90.0;
phiCTB01=(phiCTB01_0 + phiCTB01_off)*pi/180;
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# SWITCHES
#--------------------------------------------------------------------------
use_new_buncher = 0;
auto_phase = 0;
space_charge = 1;
space_charge_type = 1;
viewscreens_on =0;
kill_backwards_particles = 0;
optimize_on = 1;
use_gpt_distribution = 0;

use_HO_bzsol = 1;

variable_iris_on = 0;
sample_iris_on = 0;
beampipe_boundary_on = 1;
anode_iris_on = 1;
buncher_boundary_on = 1;
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# INITIAL PARTICLE DISTRIBUTION
#--------------------------------------------------------------------------
xoffset=0;
yoffset=0;

if (auto_phase == 1) {

	setparticles("beam", 10, me, qe, 0);
	setxoffset("beam",xoffset);
	setyoffset("beam",yoffset);
	mc2 = m*c*c/abs(qe);
	KE = 1e-3;
	G  = 1+KE/mc2 ;
	setGdist( "beam", "u", G, 0) ;

} else if (use_gpt_distribution == 1) {
	
        n_particles = 6000;
	Qtot = -1.0e-12; 		# -1.0e-13;
	xyrms = 0.5e-3;  		# [m]
	trms = 8.0e-12;  		# [s]
	MTE = 0.0025;      		# [eV]

	setparticles("beam", n_particles, me, qe, Qtot);

	# Tranversely Gaussian:
	setxdist("beam","g",0,xyrms,5,5);
	setydist("beam","g",0,xyrms,5,5);

        # Flat top longitudinally:
	tlength = sqrt(12)*trms;
	settdist("beam","u",0,tlength);
	
	mc2 = me*c*c/abs(qe);

        sigma_GBx = sqrt(MTE/mc2);
	sigma_GBy = sqrt(MTE/mc2);
	sigma_GBz = sqrt(MTE/mc2);

	setGBxdist("beam","g",0,sigma_GBx,5,5);
        setGBydist("beam","g",0,sigma_GBy,5,5);  
	setGBzdist("beam","u",sigma_GBz,0);

} else {
	setfile("beam", "gpt_particles.gdf");
}
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# ACCURACY AND INTEGRATION PARAMETERS
#--------------------------------------------------------------------------
RadiusMax = 0.04;
# Warning: 2000 particles is converged in B1 section only with GBacc > 5.5
GBacc =5.5;
xacc  =6.5;
accuracy(GBacc, xacc);
dtmin = 1e-16;

if(sample_iris_on==1){
	dtmax = 2e-12;
} else {
	dtmax = 20e-12;
}
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# BEAMLINE DEFINITION
#--------------------------------------------------------------------------
# DC GUN
Zdcgun = 0.0;
anode_z_position = 0.02;
Maxdcgun  = 300.0;
Edcgun    = gun_voltage / Maxdcgun;
Map2D_E("wcs", "z", Zdcgun, "fields/dcgun_analytic_5_20mm_2D.gdf", "R", "Z", "Er", "Ez", Edcgun);

# SOL1
SLAPolarity = -1;
LSLA01 = 0.6;
XSLA01 = 0.0;
YSLA01 = 0.0;
ZSLA01 = 0.3;
#BSLA01 = sol_1_current * SLAPolarity;
#Map1D_B("wcs", XSLA01,YSLA01,ZSLA01,  1,0,0, 0,1,0, "fields/Utoronto_solenoid.gdf", "Z", "Bz", BSLA01);
#Map2D_B("wcs", XSLA01,YSLA01,ZSLA01,  1,0,0, 0,1,0, "fields/scaled_L0_sol.gdf", "R", "Z", "Br", "Bz", BSLA01);
#Map1D_B("wcs", XSLA01,YSLA01,ZSLA01,  1,0,0, 0,1,0, "fields/scaled_L0_sol_1D.gdf", "z", "Bz", BSLA01);

# ANALYTIC MAP:
turns = 27000;
R =  0.0294;
L =  0.0332;
if(use_HO_bzsol==1) {
	bzsolenoidO3("wcs", "z", ZSLA01, R, L, turns*sol_1_current);
} else {
	bzsolenoid("wcs", "z", ZSLA01, R, L, turns*sol_1_current);
}

# BUNCHER
gammaCTB01=1;
XCTB01 = 0.00;
YCTB01 = 0.00;
ZCTB01 = 0.6;

if (use_new_buncher == 1) {
	Master_RF = 3.9e9;
	MaxECTB = 431.3825;  # [MV/m]
	UnitsECTB = 0.1305;  # [MV/m / kV]
	ECTB01 = buncher_voltage / MaxECTB * UnitsECTB;
	map3D_Ecomplex("wcs", XCTB01,YCTB01,ZCTB01, 1,0,0, 0,1,0, "fields/3.9GHzBuncher3D_E.gdf", "x","y","z", "ExRe","EyRe","EzRe", "ExIm","EyIm","EzIm", ECTB01, phiCTB01, 2*pi*Master_RF);
	map3D_Hcomplex("wcs", XCTB01,YCTB01,ZCTB01, 1,0,0, 0,1,0, "fields/3.9GHzBuncher3D_H.gdf", "x","y","z", "HxRe","HyRe","HzRe", "HxIm","HyIm","HzIm", ECTB01, phiCTB01, 2*pi*Master_RF);    
        buncher_length = 0.031;
	buncher_radius = 0.0025;    
} else {
	Master_RF = 3.0e9;
	MaxECTB = 6.2250;  # [MV/m]
	UnitsECTB = 1.0/6.0459; # [MV/m / kV]
	ECTB01 = buncher_voltage / MaxECTB * UnitsECTB;
	Map25D_TM("wcs", XCTB01,YCTB01,ZCTB01,  1,0,0, 0,1,0, "fields/eindhoven_rf_4mm_center.gdf", "R", "Z", "Er", "Ez", "H", ECTB01, 0, phiCTB01, 2*pi*Master_RF);
        #Map1D_TM("wcs", XCTB01,YCTB01,ZCTB01,  1,0,0, 0,1,0, "fields/eindhoven_1D.gdf", "Z", "Ez", ECTB01, phiCTB01, 2*pi*Master_RF);
        buncher_length = 0.04;
	buncher_radius = 0.002;
}


# Minimum distance between buncher and solenoid 2
dZmin_buncher_to_sol_2 = 0.1; 

# SOL2
LSLA02 = 0.6;
XSLA02 = 0.0;
YSLA02 = 0.0;
ZSLA02 = 0.8;
BSLA02 = sol_2_current * SLAPolarity;
#Map1D_B("wcs", XSLA02,YSLA02,ZSLA02,  1,0,0, 0,1,0, "fields/Utoronto_solenoid.gdf", "Z", "Bz", BSLA02);
#Map2D_B("wcs", XSLA02,YSLA02,ZSLA02,  1,0,0, 0,1,0, "fields/scaled_L0_sol.gdf", "R", "Z", "Br", "Bz", BSLA02);
#Map1D_B("wcs", XSLA02,YSLA02,ZSLA02,  1,0,0, 0,1,0, "fields/scaled_L0_sol_1D.gdf", "z", "Bz", BSLA02);

# ANALYTIC MAP:
turns = 27000;
R =  0.0294;
L =  0.0332;
if (use_HO_bzsol==1) {
	bzsolenoidO3("wcs", "z", ZSLA02, R, L, turns*sol_2_current);
} else {
	bzsolenoid("wcs", "z", ZSLA02, R, L, turns*sol_2_current);
}
#-----------------------------------------------------------------------


#-----------------------------------------------------------------------
# APERATURES
#-----------------------------------------------------------------------
beam_pipe_radius = 0.0254/2;  # 1" beampipe

# Anode: WARNING THIS SHOULD MATCH UP WITH THE CORRECT GUN MAP
anode_diameter = 0.005;
forwardscatter("wcs","I","remove",0);

if(anode_iris_on==1) {
	scatteriris("wcs", "z", anode_z_position, anode_diameter/2, beam_pipe_radius) scatter="remove";    
}

z_sample = 1;


# Beam Pipe:
beam_pipe_length = z_sample;

if(beampipe_boundary_on ==1 ) {
	scatterpipe("wcs","I",0,beam_pipe_length,beam_pipe_radius) scatter="remove";  
}

# Buncher Field Map:
#buncher_length = 0.04;
#buncher_radius = 0.002;

if(buncher_boundary_on == 1) {
	scatteriris("wcs", "z", ZCTB01 - 0.5*buncher_length, buncher_radius,beam_pipe_radius) scatter="remove";
	scatterpipe("wcs","I",ZCTB01 - 0.5*buncher_length,ZCTB01 + 0.5*buncher_length,buncher_radius) scatter="remove";
}
#-----------------------------------------------------------------------


#--------------------------------------------------------------------------
# SPACECHARGE PARAMETERS
#--------------------------------------------------------------------------
Alpha  = 1.2;   # 1.2    (manual default is 1.0), (bigger is better)
Fn     = 0.25;  # 0.25   (manual default is 0.5), (smaller is better)
verror = 0.005; # 0.005  (manual default is 0.01), (smaller is better)
Nstd   = 12;    # 12     (manual default is 5), (bigger is better)
if (space_charge == 1) {
	if (space_charge_type == 0) {
		setcharge2Dcircle("beam", total_charge);
		spacecharge2Dcircle();
	}
	if (space_charge_type == 1) {
		spacecharge3Dmesh("Cathode", "MeshNfac", Alpha, "MeshAdapt", Fn, "SolverAcc", verror, "MeshBoxSize", Nstd);
	}
}
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# OUTPUT PARAMETERS
#--------------------------------------------------------------------------
tmax = 70e-9;

if (auto_phase == 1) {
	zmax = z_sample+0.1;
	screen("wcs","I",z_sample);

	ZSTART = -0.005;
	zminmax("wcs","I",ZSTART,zmax);
} else if (optimize_on == 1) {
	zmax = z_sample + 0.1;  # put zmax 1 cm after the output screen
	screen("wcs","I",z_sample);
	ZSTART = -0.005;
	zminmax("wcs","I",ZSTART,zmax);
} else {
	ZSTART = -0.005;
	zminmax("wcs","I",ZSTART,ZSTOP);
	dtout =2e-10;
	tout(0,tmax,dtout);
	if (kill_backwards_particles == 1) {
	#GBzmin("wcs", "I", 5.0e-10, 0.0); # kill particles going backwards (gamma*beta < 0.0) after t=5e-10 s
	}
	# stdxyzmax(0,0,5) ; # remove particles if z > 5*stdz from center of bunch
}
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# AUTOPHASING VARIABLES
#--------------------------------------------------------------------------
phasing_amplitude_0 = buncher_voltage;
phasing_on_crest_0 = phiCTB01_0;
phasing_relative_0 = phiCTB01_off;
phasing_gamma_0 = gammaCTB01;
#--------------------------------------------------------------------------