OPTION, ECHO=TRUE;
OPTION, INFO=TRUE;
OPTION, PSDUMPFREQ = 0;     // 6d data never written (h5).
OPTION, STATDUMPFREQ = 10;  // Beam Stats written every 10 time steps (stat).
OPTION, BOUNDPDESTROYFQ = 10; // Delete lost particles, if any
OPTION, VERSION = 20300;    // This indicates that it works with an OPAL version of 2.3

Title, string="Testing Opal Bend Magnets";

//----------------------------------------------------------------------------
//Global Parameters

REAL rf_freq             = 352.2;     //RF frequency. (MHz)
REAL n_particles         = 1E4;      //Number of particles in simulation.  Not used when importing a distribution.

// Beam properties
// from NIST, the following values are found:
// alpha mass = 3.7274 GeV/c^2
// deuteron mass = 1.8756 GeV/c^2 (Opals internal value is -0.7% smaller than this value)
REAL THEMASS   = 1.8756;  // particle mass
REAL THECHARGE = 1;       // number of proton charges

//Initial Momentum Calculation
REAL Edes    = 25.0e-3; //initial energy in GeV
REAL gamma   = (Edes+THEMASS)/THEMASS; 
REAL beta    = sqrt(1-(1/gamma^2));
REAL P0      = gamma*beta*THEMASS;    //inital z momentum
REAL BRHO    = P0 * 1e9 / (THECHARGE * CLIGHT); //initial rigidity (T*m)


REAL sigmat  = CLIGHT * 1E-19; // initial bunch length - s
REAL sigmaz  = beta * CLIGHT * sigmat; // initial bunch length - m

// As of 2.2.1 Opal incorrectly computes the bunch charge in TrackRun.cpp as charge*current
// in order to workaround this problem it was recommended that I halve the current for 
// doubly charged particles, BCURRENT is only used to compute the bunch charge
ABEAM: BEAM, PARTICLE = nusano, pc = P0, NPART = n_particles,
             BFREQ = rf_freq, BCURRENT = 0.001,
             CHARGE = THECHARGE, MASS = THEMASS;

//Printing initial energy and momentum to terminal output.
value , {Edes, gamma, beta, P0, BRHO};

//----------------------------------------------------------------------------
// Drifts


//----------------------------------------------------------------------------
// the bend

REAL switch_angle = 15.0 * Pi / 180.0; // 15-degree switch angle
REAL switch_radius = 1; // 1 meter bending radius

ABEND : SBEND, ANGLE = switch_angle,
		FMAPFN = "1DPROFILE1-DEFAULT",
		ELEMEDGE = 0.25,
		DESIGNENERGY = Edes * 1e3,
		L = 2 * switch_radius * sin(switch_angle / 2.0),
		E1 = 0, E2 = 0, // make a "true" sector magnet
		HAPERT = 0.2,
		GAP = 0.02;


//----------------------------------------------------------------------------
// Quadrupoles

//----------------------------------------------------------------------------
// DEFINE BEAM LINE

// Just a bend
LINEA : LINE = (ABEND);


//----------------------------------------------------------------------------
// INITIAL DISTRIBUTION

Dist: DISTRIBUTION, TYPE = GAUSS,
	SIGMAX = 1E-3, SIGMAY = 1E-3, SIGMAZ = 1E-3,
	SIGMAPX = 0, SIGMAPY = 0, SIGMAPZ = 0,
	EMITTED = FALSE,
	EMISSIONMODEL = NONE,
	EMISSIONSTEPS = 16,
	WRITETOFILE = True;    //Saves the distribution to a text file

//----------------------------------------------------------------------------
// Define Field solvers
// The mesh sizes should be a factor of 2 
// for most efficient space charge calculation.

FS_SC: Fieldsolver, FSTYPE = NONE, 
            MX = 16, MY = 16, MT = 16,
            PARFFTX = False, 
            PARFFTY = False, 
            PARFFTT = True,
            BCFFTX = open, 
            BCFFTY = open, 
            BCFFTT = open,
            BBOXINCR = 4, 
            GREENSF = INTEGRATED;

//----------------------------------------------------------------------------
// Simulate the beamline using TRACK and RUN.

TRACK, LINE = LINEA, BEAM = ABEAM, MAXSTEPS = 50000, 
    DT = 2.1e-11, ZSTOP=1.0;

RUN, METHOD = "PARALLEL-T", BEAM = ABEAM, 
    FIELDSOLVER = FS_SC, DISTRIBUTION = Dist;
ENDTRACK;

Quit;