// source: rfq-dip-wiki.lns.mit.edu:/data01/IsoDAR/mariay_old/IsoDAR/tmpl/IsoDAR.tmpl
// Updated input file for OPAL 2.0 compatibility
// add OPTION, VERSION = 20000
// OPTION, PSDUMPLOCALFRAME has been removed, so OPTION, PSDUMPFRAME = BUNCH_MEAN should be used
// H2P has been added, in BEAM set PARTICLE = H2P and remove MASS and CHARGE
// SPACECHARGE in BEAM no longer exists, has to be removed (does nothing)

OPTION, VERSION = 20000; 		// OPAL 2.0
OPTION, ECHO = FALSE;			// default FALSE
OPTION, PSDUMPFREQ=1;
OPTION, PSDUMPFRAME = GLOBAL; //BUNCH_MEAN;	// OPTION, PSDUMPLOCALFRAME = true;
OPTION, SCSOLVEFREQ = 2;
OPTION, REMOTEPARTDEL = -1; // default is -1 (no delete)
OPTION, TELL = TRUE;				// show OPTION settings
//OPTION, ASCIIDUMP = TRUE;             // FALSE for probe .h5 istead of .loss
TITLE, STRING = "IsoDAR Accelerated Orbit (harm=4)";	// title of the simulation, "IsoDAR Track Bunch";


// CYCLOTRON
REAL harmon = 4;
REAL frequency = 32.8; // RF frequency in [MHz], particle oscillation frequency is frequency/harmo
REAL r0= 277.013;
REAL phi0= 0.0;
REAL z0= 0.0;
REAL pr0= 41.0417e-5;
REAL pz0= 0.0;
REAL minz= -40.0;
REAL maxz= 40.0;
REAL minr= 0.0;
REAL maxr= 3000.0;

Inj60: CYCLOTRON, TYPE="CARBONCYCL", CYHARMON=harmon, RFFREQ=frequency,
        FMAPFN="Injmap20120405_Symmetrized.dat", SYMMETRY=1.0,
        RINIT=r0, PHIINIT=phi0, ZINIT=z0, PRINIT=pr0, PZINIT=pz0,
        MINZ=minz, MAXZ=maxz, MINR=minr, MAXR=maxr;


// CAVITIES
REAL cavity_angle = 40.8924;
REAL phi1a = 159.136;
REAL phi1b = phi1a+180;
REAL phi2a = phi1a;
REAL phi2b = phi1a+180;
REAL phi3a = phi1a;
REAL phi3b = phi1a+180;
REAL phi4a = phi1a;
REAL phi4b = phi1a+180;
REAL vol = 0.250; //MV
REAL rmin = 40.0; //mm 
REAL rmax = 2100.0; //mm
REAL gapwidth = 10; //mm

rf1a: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=45.0-cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi1a;
rf1b: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=45.0+cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi1b;
rf2a: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=135.0-cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi2a;
rf2b: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=135.0+cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi2b;
rf3a: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=225.0-cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi3a;
rf3b: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=225.0+cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi3b;
rf4a: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=315.0-cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi4a;
rf4b: RFCavity, VOLT=vol, FMAPFN="InjRFVoltage.dat", TYPE="SINGLEGAP", FREQ=frequency,
        RMIN=rmin, RMAX=rmax, ANGLE=315.0+cavity_angle/2, PDIS = 0.0, GAPWIDTH = Gapwidth, PHI0=phi4b;
rfCavities: LINE=(rf1a,rf1b,rf2a,rf2b,rf3a,rf3b,rf4a,rf4b);


// LINE
l1: LINE = (Inj60,rfCavities);


// DISTRIBUTION
Dist1:DISTRIBUTION, TYPE=GAUSSMATCHED, WRITETOFILE=TRUE, LINE=l1,
    MAGSYM=1, FMAPFN="Injmap20120405_Symmetrized.dat", FMTYPE="CARBONCYCL";
//Dist1:DISTRIBUTION, TYPE=FROMFILE, FNAME="1Particle.dat", INPUTMOUNITS=None;


// FIELDSOLVER
REAL Gridsize= 16;
Fs1:FIELDSOLVER, FSTYPE=FFT, MX=Gridsize, MY=Gridsize, MT=Gridsize,
                 PARFFTX=FALSE, PARFFTY=FALSE, PARFFTT=TRUE,
                 BCFFTX=OPEN, BCFFTY=OPEN, BCFFTT=OPEN;


// BEAM
REAL H2MASS=2.01465054*931.49432/1000.0; //GeV
REAL Edes=1.0*2.01465054/1000;
REAL gamma=(Edes+H2MASS)/H2MASS;
REAL beta=sqrt(1-(1/gamma^2));
REAL P0 = gamma*beta*H2MASS;

beam1: BEAM, PC=P0, PARTICLE=H2P, NPART=1e4, BCURRENT=6.65e-3,
    BFREQ=frequency;


// TRACK
REAL turns = 0;
REAL nstep=2000;

SELECT, LINE=l1;
TRACK,LINE=l1, BEAM=beam1, MAXSTEPS=nstep*turns, STEPSPERTURN=nstep, TIMEINTEGRATOR="RK-4";
RUN, METHOD="CYCLOTRON-T", BEAM=beam1, FIELDSOLVER=Fs1, DISTRIBUTION=Dist1;
ENDTRACK;

// PYTHON SCRIPT FOR OPTIMIZER
//SYSTEM, CMD='python3 ../turn_separation.py';

STOP;
QUIT;