The OPAL Discussion Forum

[Jochem Snuverink <jochem.snuverink AT psi.ch>]

Dear Mael,

As Tomas mentioned there is a somewhat experimental standalone tune calculation in OPAL.

In the input file you set:

OPTION, CLOTUNEONLY =true;

In addition you need to specify "FMLOWE", "FMHIGHE" ini the CYCLOTRON command. Those specify the lower and upper limit of the energy in MeV, i.e. FMLOWE=72 and FMHIGHE=590 for the PSI RingCyclotron.

In order to compute the tunes, you then specify the distribution command:

Dist:DISTRIBUTION, TYPE=GAUSSMATCHED, LINE=Line1, NSTEPS=720, SECTOR=FALSE;

SECTOR = FALSE ------> FALSE = use 360 degrees of the machine
NSTEPS=720         -------> number of integration steps of closed orbit finder

Attached an input file example for the Ring Cyclotron.

As mentioned by Tomas the method doesn't always seem to converge, and we are looking into that.

Best regards,
Jochem

On 12/03/19 18:05, Tomas Matlocha wrote:

Hi Mael,

the ic.dat file is generated by FIXPO code which is external code for equilibrium orbit calculation and is not part of OPAL. 

The Tune calculation implementation to OPAL-cycl itself is currently in progress. It works well for machines with higher tunes i.e. PSI Ring, unfortunately for small machines (my case) the results are not fully relevant yet.

I'm sure there will be some skilled OPAL user who can give you more precise answer.

Best regards,

Tomas

 

 

Dne 2019-03-11 16:21, Adelmann Andreas (PSI) napsal:

 

Begin forwarded message:

From: mael.fichou <mael.fichou AT wanadoo.fr>

Subject: Question about OPAL input file

Date: 11 March 2019 at 15:06:51 CET

To: <andreas.adelmann AT psi.ch>, <achim.gsell AT psi.ch>

Cc: Vitcher Sarah <vitcher.sarah AT orange.fr>

Reply-To: mael.fichou <mael.fichou AT wanadoo.fr>

Hello,   

I am Mael Fichou, french engineer working on cyclotron project. In order to evaluate the use of your OPAL software for our products studies, I am trying to run a tune simulation and am wondering how you obtained the ic.dat file, Could you advise on what kind of software I could use to obtain these information ? Checking a file example, I don't understand how you generate the third parameter (Pr).  

Thank you very much for your help,   Sincerely,   Mael Fichou

// Ring template file
// started 21-01-16: Nate Pogue with Andreas Adelmann
//
//---------------------------------------------------------------------------------------------//
//                                      Global Options                        
                 //
//---------------------------------------------------------------------------------------------//
OPTION, VERSION         =10900;         // Version number, acknowledging 
recent changes
OPTION, ECHO            =false;         // If true, command prints an echo of 
the input files on
                                        // the standard error file.
OPTION, PSDUMPFREQ      =1000000000;    // Defines after how many steps the 
phase is dumped into
                                        // the H5hut file, Default is 10. A 
higher number dumps
                                        // less. if 1 then every step, every 
turn. If too low it
                                        // will be time prohibitive, and 
large file sizes.
OPTION, STATDUMPFREQ    =1000000000;  // Defines after how many steps the 
bunch statistics is
                                        // dumped into *.stat file. The 
behaviour is like
                                        // PSDUMPFREQ.
OPTION, SPTDUMPFREQ     =1000000000;    // Defines after how many steps we 
dump the phase space
                                        // for Single Particle simulations. 
It is always useful
                                        // record the trajectory of the 
reference particle or
                                        // some specified particle for a 
primary study.
                                        // Default value is 1. same as above 
but single particle
OPTION, REPARTFREQ      =1000000000;    // Defines after how many time steps 
the aprticles are
                                        // repartitioned to balacne the 
computational load on
                                        // the computer nodes, Default value 
is 10.
OPTION, PSDUMPFRAME     =global; // Determines whether the phase space data 
is dumped in
                                        // global frame or in the local 
frame. If true, run can
                                        // not be restarted as it is in local 
frame. Default is
                                        // false
OPTION, PSDUMPEACHTURN  =false; // Phase space dumped each turn
OPTION, ENABLEHDF5      =false;  // Enable HDF5 (.h5) output
OPTION, CLOTUNEONLY     =true;

//---------------------------------------------------------------------------------------------//
//                                      Physical Parameters                   
                 //
//---------------------------------------------------------------------------------------------//
REAL Ebeam=0.072;                       // injection energy of the beam
REAL gamma              =(Ebeam+PMASS)/PMASS;           // relativistic gamma
REAL beta               =sqrt(1-(1/gamma^2));           // relativistic beta
REAL gammabeta          =gamma*beta;                    // self explanitory
REAL P0                 =PMASS*gammabeta;               // Initial momentum 
of the particle
REAL Brho               =(PMASS*1.0e9*gammabeta)/CLIGHT;// Magnetic rigidity

VALUE,{gamma,beta,gammabeta,Ebeam,Brho,CLIGHT};         // Prints values in 
output to check
                                                        // initial conditions 

//---------------------------------------------------------------------------------------------//
//                                      Cyclotron Configuration               
                 //
//---------------------------------------------------------------------------------------------//
// Below are the series of value and commmands required to set up the 
cyclotron's configuration

// Trim Coils
IF (true) {
  // Fitted Trim coils
  // Fits obtained with https://gitlab.psi.ch/frey_m/TrimCoil
  tc1:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 1000.000, RMAX = 2160.534, 
BMAX=-0.00413729528524,
        COEFNUM   = {-77.6582083367, 147.483554457, -99.356050989, 
28.0510488454, -2.75931563027},
        COEFDENOM = {11087.3671263, -19945.0494828, 14265.1095084, 
-5091.97533455, 907.971669105, -63.9544501666};
  tc2:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 1997.361, RMAX = 2240.789, 
BMAX=0.00218161504809,
        COEFNUM   = {41.5531008893, -28.3951208173, 4.12163299162},
        COEFDENOM = {20064.0752201, -31211.3212264, 18120.6175688, 
-4681.4081011, 459.349240061};
  tc3:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 1934.749, RMAX = 2327.682, 
BMAX=-0.0017577528581,
        COEFNUM   = {2162.9968198, -1475.07674541, 213.49992936},
        COEFDENOM = {320349.135485, -461518.079087, 245420.010296, 
-57350.7810108, 5055.03942358};
  tc4:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 2199.249, RMAX = 2500.693, 
BMAX=-0.00292918145657,
        COEFNUM   = {-0.0311184855434, 0.0304018467119, -0.00734268403986},
        COEFDENOM = {15.5004996028, -22.1852773857, 11.900022194, 
-2.8450478635, 0.257261375807};
  tc5:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 2421.128, RMAX = 2694.377, 
BMAX=0.00323523483991,
        COEFNUM   = {-0.0232820414694, 0.0233635767522, -0.00560765583653},
        COEFDENOM = {27.3085153497, -37.2899052844, 19.1033049888, 
-4.36170628039, 0.375822127348};
  tc6:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 2608.590, RMAX = 2894.288, 
BMAX=0.00300890390296,
        COEFNUM   = {0.0147589213103, -0.00415160359162, -0.000467283940805},
        COEFDENOM = {37.3333889061, -47.5817574034, 22.7123326601, 
-4.81980947963, 0.384648033613};
  tc7:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 2797.717, RMAX = 3095.748, 
BMAX=0.000271243551746,
        COEFNUM   = {0.135330988621, -0.0729640701831, 0.00912649095538},
        COEFDENOM = {105.252930097, -125.050150643, 55.5661776884, 
-10.957992747, 0.810971100419};
  tc8:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 2986.646, RMAX = 3301.785, 
BMAX=-0.00193185058516,
        COEFNUM   = {29.3741065914, -15.9219525812, 2.08824908481},
        COEFDENOM = {13623.8682088, -15049.7382343, 6207.78926511, 
-1134.45457254, 77.6713919544};
  tc9:  TRIMCOIL, TYPE="PSI-PHASE", RMIN = 3181.309, RMAX = 3486.911, 
BMAX=-0.00168781264499,
        COEFNUM   = {7.12521770929, -3.71905051895, 0.473857352798},
        COEFDENOM = {4043.07154766, -4303.44292303, 1715.86127048, 
-304.094875025, 20.2496511554};
  tc10: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 3358.768, RMAX = 3690.307, 
BMAX=0.00114708986366,
        COEFNUM   = {3.74489923333, -1.81983152935, 0.215008386857},
        COEFDENOM = {2095.42645051, -2236.68521207, 902.642501832, 
-163.409794152, 11.2094776696};
  tc11: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 3591.294, RMAX = 3852.300, 
BMAX=0.00186152244639,
        COEFNUM   = {0.325610070561, -0.028811308006, -0.0155305827246},
        COEFDENOM = {1338.24663041, -1552.17121428, 680.400760398, 
-133.278768884, 9.82424577869};
  tc12: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 3792.200, RMAX = 4013.892, 
BMAX=-0.000988980040693,
        COEFNUM   = {0.111326065649, 0.0704809114777, -0.0252481936502},
        COEFDENOM = {1529.85640715, -1705.40431001, 718.894828043, 
-135.437536308, 9.60031844855};
  tc13: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 3943.861, RMAX = 4165.899, 
BMAX=-0.00493430409585,
        COEFNUM   = {2.03548211073, -0.902630357138, 0.098962673965},
        COEFDENOM = {1877.17709245, -1974.80545458, 785.559287661, 
-139.775458035, 9.37049022082};
  tc14: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 4093.048, RMAX = 4439.405, 
BMAX=0.00164689077041,
        COEFNUM   = {0.0336023039334, -0.0134183623486, 0.00130460201021},
        COEFDENOM = {21.2519341602, -21.8526210442, 8.50784116323, 
-1.48233040122, 0.09729913873};
  tc15: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 4324.201, RMAX = 4560.073, 
BMAX=-0.000541564210784,
        COEFNUM   = {-0.0312020990404, 0.0227946756108, -0.00354827255973},
        COEFDENOM = {14.7460286849, -16.9186605846, 7.61516943548, 
-1.53074181639, 0.11384470123};
  tc16: TRIMCOIL, TYPE="PSI-PHASE", RMIN = 4366.692, RMAX = 4533.039, 
BMAX=0.00453576803859,
        COEFNUM   = {-117.893003406, 154.620805581, -74.0713608326, 
15.4791314127, -1.19200157147},
        COEFDENOM = {17960.3792266, -15738.5143071, 3313.45362486, 
668.775387636, -317.372662563, 29.7919862341};

} ELSE {
  // parametrised trim coils
  tc1:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 1000.000, RMAX = 
2160.534, BMAX=-0.00413729528524, SLPTC=1;
  tc2:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 1997.361, RMAX = 
2240.789, BMAX=0.00218161504809, SLPTC=1;
  tc3:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 1934.749, RMAX = 
2327.682, BMAX=-0.0017577528581, SLPTC=1;
  tc4:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 2199.249, RMAX = 
2500.693, BMAX=-0.00292918145657, SLPTC=1./6.;
  tc5:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 2421.128, RMAX = 
2694.377, BMAX=0.00323523483991, SLPTC=1./6.;
  tc6:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 2608.590, RMAX = 
2894.288, BMAX=0.00300890390296, SLPTC=1./6.;
  tc7:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 2797.717, RMAX = 
3095.748, BMAX=0.000271243551746, SLPTC=1./6.;
  tc8:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 2986.646, RMAX = 
3301.785, BMAX=-0.00193185058516, SLPTC=1./6.;
  tc9:  TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 3181.309, RMAX = 
3486.911, BMAX=-0.00168781264499, SLPTC=1./6.;
  tc10: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 3358.768, RMAX = 
3690.307, BMAX=0.00114708986366, SLPTC=1./6.;
  tc11: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 3591.294, RMAX = 
3852.300, BMAX=0.00186152244639, SLPTC=1./6.;
  tc12: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 3792.200, RMAX = 
4013.892, BMAX=-0.000988980040693, SLPTC=1./6.;
  tc13: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 3943.861, RMAX = 
4165.899, BMAX=-0.00493430409585, SLPTC=1./6.;
  tc14: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 4093.048, RMAX = 
4439.405, BMAX=0.00164689077041, SLPTC=1./6.;
  tc15: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 4324.201, RMAX = 
4560.073, BMAX=-0.000541564210784, SLPTC=1./6.;
  tc16: TRIMCOIL, TYPE="PSI-BFIELD-MIRRORED", RMIN = 4366.692, RMAX = 
4533.039, BMAX=0.00453576803859, SLPTC=1./6.;
}

// Note: Trim coil 1 must be larger than trim coil 2 and trim coil 2 must be 
larger than
//       trim coil 3. 

Ring:      CYCLOTRON, TYPE="RINGCYC", CYHARMON=6, PHIINIT=110.154722071, 
PRINIT=-0.0148229176092,
           RINIT=2043.34435172,SYMMETRY=8.0, RFFREQ=50.65, BSCALE=1, 
FMAPFN="s03av.nar",
           FMLOWE=72, FMHIGHE=590, TRIMCOILTHRESHOLD=1.7,
           TRIMCOIL={tc1,  tc2,  tc3,  tc4, tc5, tc6, tc7, tc8, tc9, tc10, 
tc11, tc12, tc13,
                     tc14, tc16};
// Note: CYHARMON - Cyclotron's harmonics, SYMMETRY - symmetry of cyclotron, 
in this case the
//       field map imported is for the whole cyclotron, hence symetrry 1, 
FMAPFN - insert
//       field maps - in this case the sector dipole field map.
// Note: In general there should be gap widths, otherwise no transit is taken 
into account. In
//       this model none are specified by one needs to add this for accurate 
models.

Main1: RFCAVITY, VOLT=(0.795377279641 + 0.25*(0.167803166248)), 
FMAPFN="MainCav", TYPE="SINGLEGAP", FREQ=50.65,
       RMIN=(1003.38+(54.8114134092)), RMAX=(6203.38+(54.8114134092)), 
ANGLE=35.3908241093, PDIS=422.960196999,
       GAPWIDTH=300.0, PHI0=(204);
Main2: RFCAVITY, VOLT=(0.839065409318 + 0.25*(0.167803166248)), 
FMAPFN="MainCav", TYPE="SINGLEGAP", FREQ=50.65,
       RMIN=(1003.38+(37.2499463707)), RMAX=(6203.38+(37.2499463707)), 
ANGLE=((34.8947228347)+90), PDIS=419.016417256,
       GAPWIDTH=300.0, PHI0=((204)+180);
Main3: RFCAVITY, VOLT=(0.869132281998 + 0.25*(0.167803166248)), 
FMAPFN="MainCav", TYPE="SINGLEGAP", FREQ=50.65,
       RMIN=(1003.38+(67.0575751737)), RMAX=(6203.38+(67.0575751737)), 
ANGLE=((35.0434440779)+180), PDIS=426.758463119,
       GAPWIDTH=300.0, PHI0=(204);
Main4: RFCAVITY, VOLT=(0.86018869306 + 0.25*(0.167803166248)), 
FMAPFN="MainCav", TYPE="SINGLEGAP", FREQ=50.65,
       RMIN=(1003.38+(47.9563058633)), RMAX=(6203.38+(47.9563058633)), 
ANGLE=((35.0055171836)+270), PDIS=416.668103168,
       GAPWIDTH=300.0, PHI0=((204)+180);

PhaseProbe: RFCAVITY, VOLT=0, FMAPFN="MainCav", TYPE="SINGLEGAP", 
FREQ=50.65*6,
       RMIN=(1003.38+(54.8114134092)), RMAX=(6203.38+(54.8114134092)), 
ANGLE=170, PDIS=889,
       GAPWIDTH=1.0, PHI0=114;     
       
// Flat Top Cavities

REAL ftfreq             =(50.65)*3;                     // 3rd harmonic 
frequency for FT Cavity

FT1:   RFCAVITY, VOLT=0.542111074151, FMAPFN="FlatCav", TYPE="SINGLEGAP", 
FREQ=ftfreq,
       RMIN=(1845.45+(-82.5690011764)), RMAX=(4477.27+(-82.5690011764)), 
ANGLE=260.005454297, PDIS=451.707502839,
       GAPWIDTH=250.0, PHI0=142.942358367;

//This is the Higher Order Mode Flat top that Flat tops on the first half of 
the cavity and
//then accelerates the beam on the second half. Reduces the turns from ~185 
to 165. 
//FT1:   RFCAVITY, VOLT=-0.542111074151, FMAPFN="NateFlatTopHOM.dat", 
TYPE="SINGLEGAP", FREQ=ftfreq,
//       RMIN=1845.45+(-82.5690011764)), RMAX=(4477.27+(-82.5690011764)), 
ANGLE=260.005454297, PDIS=451.707502839,
//       GAPWIDTH=250.0, PHI0=142.942358367;

//---------------------------------------------------------------------------------------------//
//                                      Probes  - Use Rudolf's positions - 
note comments       //
//---------------------------------------------------------------------------------------------//
// RRE4 - extraction probe, contains the last ~10 turns, is a wire with 
stopping block

   // Bi's position
   // RRE4: PROBE, XSTART= 3809.5, XEND= 4140.4, YSTART= -1610.2, YEND= 
-1668.5, STEP=1.0;

   // Martin Humbel positions
   // RRE4: PROBE, XSTART= 3800, XEND= 4150, YSTART= -1590.85, YEND= 
-1674.88, STEP=1.0;

   // Rudolf Doelling position
   RRE4: PROBE, XSTART= 3877.22614, XEND= 4172.66847, YSTART= -1622.1169, 
YEND= -1674.2113, STEP=1.0;
   
   RRE4new: PROBE, XSTART= 3882.41822, XEND= 4177.86055, YSTART= -1592.6711, 
YEND= -1644.7656,STEP=1.0;
   // Note that in this position the orbits are shifted by 96.7 mm. The 
current reading from
   // the probe are shifted by this amount. Shift to the old position matches 
the simulation,
   // thus one needs to shift over to the old position as it seems to 
accurately describe the
   // machine.
   // Until then here is an alignment that is close.
   //RRE4: PROBE, XSTART= 4123.78, XEND= 4271.509, YSTART= -1087.61, YEND= 
-1113.65, STEP=1.0;

// RRI2 - injections probe, contains about the first 15 turns in the machine

   // Bi's position
   // RRI2: PROBE, XSTART= -1956.8, XEND= -2477.0, YSTART= 285.61, YEND= 
358.727, STEP=1.0;

   // Martin Humbel position
   // RRI2: PROBE, XSTART= -1956, XEND= -2500.0, YSTART= 278.38, YEND= 355.8, 
STEP=1.0;

   // Rudolf Doelling position

   REAL rri2length = 500;
   REAL rri2ss  = 1000 * (2.02870163815);
   REAL rri2a   = 10 * (0.906056562159);
   REAL rri2phi = 100 * (0.0300740456504);
   REAL rri2x0  =   rri2a * SIN(rri2phi);
   REAL rri2y0  = - rri2a * COS(rri2phi);
   REAL rri2xs  = rri2x0 + rri2ss * COS(rri2phi);
   REAL rri2ys  = rri2y0 + rri2ss * SIN(rri2phi);
   REAL rri2xe  = rri2x0 + (rri2ss + rri2length) * COS(rri2phi);
   REAL rri2ye  = rri2y0 + (rri2ss + rri2length) * SIN(rri2phi);

   RRI2: PROBE, XSTART= rri2xs, XEND= rri2xe, YSTART= rri2ys, YEND= rri2ye, 
STEP=1.0;
   //RRI2: PROBE, XSTART= -1979.0748, XEND= -2474.2089, YSTART= 288.744017, 
YEND= 358.330567, STEP=1.0;

// RRL1 - long probe, covers close to all 186 turns except for the injection.

   // Bi's position
   // RRL1: PROBE, XSTART= -1800, XEND= -4202.4, YSTART= 892.31, YEND= 
1315.9, STEP=1.0;
   
   // Martin Humbel positions
   // RRL1: PROBE, XSTART= -1800, XEND= -4500, YSTART= 920.77, YEND= 1300.23, 
STEP=1.0;

   // Rudolfs position (start at r = 2128 mm (turn 8), impact parameter a = 
566.2)
   // RRL1: PROBE, XSTART= -1997.3513, XEND= -4234.8345, YSTART= 927.1215, 
YEND= 1321.65013, STEP=1.0;
   // THis one as well needs to be shifted by 32.9 mm as described above, if 
not shifted,
   // the machine can not reach 590 MeV, and the last three turns are missing.

   // Start from turn 1


   REAL rrllength = 4400;
   REAL rrlss  = 1000 * (0.34439075226);
   REAL rrla   = 10 * (56.9641442937);
   REAL rrlphi = 100 * (0.0294156599045);
   REAL x0   =   rrla * SIN(rrlphi);
   REAL y0   = - rrla * COS(rrlphi);
   REAL xs   = x0 + rrlss * COS(rrlphi);
   REAL ys   = y0 + rrlss * SIN(rrlphi);
   REAL xe   = x0 + (rrlss + rrllength) * COS(rrlphi);
   REAL ye   = y0 + (rrlss + rrllength) * SIN(rrlphi);

   RRL1: PROBE, XSTART= xs, XEND= xe, YSTART= ys, YEND= ye, STEP=1.0;

   //RRL1: PROBE, XSTART= 98.3, XEND= -4234.8345, YSTART= 557.6, YEND= 
1321.65013, STEP=1.0;


//---------------------------------------------------------------------------------------------//
//                                      Distributions                         
                 //
//---------------------------------------------------------------------------------------------//

// Note: Distribution in eV need to be specied as EV in the parameter 
INPUTMOUNITS (see Opal
//       Manual). Some from other sources are in BETA*GAMMA, thus have no 
units.

//For Single Particle
DistSP: DISTRIBUTION, TYPE = fromfile, FNAME = "spdist.opal",
        INPUTMOUNITS = EV;
//For Multiple Particles
DistSC: DISTRIBUTION, TYPE = fromfile, FNAME = "scdist.opal",
        INPUTMOUNITS = NONE;
//For Tune 
DistT:  DISTRIBUTION, TYPE = fromfile, FNAME = "tdist.opal";

//Binomial Distribution
DistB: DISTRIBUTION, TYPE = binomial,
       SIGMAX = 3.495e-03, SIGMAPX = 44.39, CORRX = 0, MX = 2,
       SIGMAY = 19.8e-03, SIGMAPY = 77.12, CORRY = 0, MY = 2,
       SIGMAZ = 1.56e-03, SIGMAPZ = 112.3, CORRZ = 0, MZ = 2,
       R61 = 0, R62 = 0, R51 = 0, R52 = 0, PT = 0,
       CUTOFFX=3, CUTOFFY=3, CUTOFFR=3, CUTOFFLONG=3,
       CUTOFFPX=3, CUTOFFPY=3, CUTOFFPZ=3;

//Gauss Distribution
DistG: DISTRIBUTION, TYPE = binomial,
       SIGMAX = 3.495e-03, SIGMAPX = 44.39, CORRX = 0,
       SIGMAY = 19.8e-03, SIGMAPY = 77.12, CORRY = 0,
       SIGMAZ = 1.56e-03, SIGMAPZ = 112.3, CORRZ = 0,
       R61 = 0, R62 = 0, R51 = 0, R52 = 0, PT = 0,
       CUTOFFX=3, CUTOFFY=3, CUTOFFR=3, CUTOFFLONG=3,
       CUTOFFPX=3, CUTOFFPY=3, CUTOFFPZ=3;

//---------------------------------------------------------------------------------------------//
//                                      Choice of Simulation                  
                 //
//---------------------------------------------------------------------------------------------//
REAL RF_MODE            =0;     // Selects Single Particle mode
REAL SC_MODE            =0;     // Selects Space Charge mode
REAL TUNE_MODE          =1;     // Selects Tune Calculation
//---------------------------------------------------------------------------------------------//
//                                      Tune Mode - Easiest to Start with.    
                 //
//---------------------------------------------------------------------------------------------//
// Calculation of the betatron oscillation frequency for a cyclotron. Calls 
tuning.bash and
// works with FIXPOEORBIT output data.
// Requires a distribution file with two particles: first particle is the 
reference orbit, the
// second is displaced off center with a dx and dz as follows:  (note dx 
radial, dz vertical)
// 2
// 0.0  0.0     0.0     0.0     0.0     0.0
// 0.01 0.0     0.0     0.0     0.01    0.0
// When the total number of particles is equal to 2, SEO mode is triggered 
automatically,
// regardless of the existance of a cavity or other components.

IF (TUNE_MODE>0){

TITLE, STRING           ="Tune Calculation of Cyclotron";

Line1: LINE = (Ring);
// the number of particles must match the number of particles stated in the 
distribution. Thus
// for tune calculations the number is just hard coded here.

REAL numparticles = 2;

Dist:DISTRIBUTION, TYPE=GAUSSMATCHED, LINE=Line1, NSTEPS=720, SECTOR=FALSE;

FS1:  FIELDSOLVER, FSTYPE = NONE, MX = 16, MY = 16, MT = 16,
      PARFFTX = true, PARFFTY = true, PARFFTT = true,
      BCFFTX  = open,      BCFFTY  = open,      BCFFTT  = open,
      BBOXINCR = 5;
}
//---------------------------------------------------------------------------------------------//
//                                      Single Particle Mode                  
                 //
//---------------------------------------------------------------------------------------------//
// Sends a single particle through the machine. Allows for quick assessment 
of the system and
// whether it works properly or as expected.

IF (RF_MODE>0){

TITLE, STRING           ="Single Particle through the Cyclotron";

Line1: LINE = (Ring, Main1, Main2, Main3, Main4, FT1, RRE4, RRL1, RRI2);

// the number of particles must match the number of particles stated in the 
distribution. Thus
// for single particle calculation the number is just hard coded here.

REAL numparticles = 1;

Dist: DistSP;
FS1:  FIELDSOLVER, FSTYPE = NONE, MX = 16, MY = 16, MT = 16,
      PARFFTX = true, PARFFTY = true, PARFFTT = true,
      BCFFTX  = open,      BCFFTY  = open,      BCFFTT  = open,
      BBOXINCR = 5;

}
//---------------------------------------------------------------------------------------------//
//                                      Space Charge Mode                     
                 //
//---------------------------------------------------------------------------------------------//
// Large Multiparticle Runs where Space Charge must be used to get correct 
dynamics.

IF (SC_MODE>0){

TITLE, STRING           ="Space Charge Simulation in the Cyclotron";

REAL Sigmacutter = 0;

IF (Sigmacutter > 0){
OPTION, REMOTEPARTDEL   = Sigmacutter;// Scrapes off the beam larger than the 
simga specified
}

Line1: LINE = (Ring, Main1, Main2, Main3, Main4, FT1, PhaseProbe, RRE4, 
RRE4new, RRL1, RRI2);

// remember the number of particles has to match what is inserted into the 
distribution.

REAL numparticles = 1;

Dist: DistB;  // Or switch to DistG for general distribution.

FS1:  FIELDSOLVER, FSTYPE = NONE, MX = 16, MY = 16, MT = 16,
      PARFFTX = true, PARFFTY = true, PARFFTT = true,
      BCFFTX  = open,      BCFFTY  = open,      BCFFTT  = open,
      BBOXINCR = 2;

}
//---------------------------------------------------------------------------------------------//
//                                      Tracking and Computation              
                 //
//---------------------------------------------------------------------------------------------//

Beam1: BEAM, PARTICLE = PROTON, PC = P0, NPART = numparticles, BCURRENT = 
.00219, CHARGE = 1.0,
      BFREQ = 50.65;

SELECT, LINE = Line1;

REAL nstep1 =2880;              //dummy for duplicate variable

TRACK, LINE = Line1, BEAM = Beam1, MAXSTEPS = ((nstep1)*(182)), STEPSPERTURN 
= nstep1,
       TIMEINTEGRATOR = "RK-4";

RUN, METHOD = "CYCLOTRON-T", BEAM = Beam1, FIELDSOLVER = FS1, DISTRIBUTION = 
Dist;

ENDTRACK;

STOP;