This contains the main PyZgoubi classes: Line and Results
The Line object holds a series of elements, to represent an accelerator lattice.
Add some extra input files to the directory where zgoubi is run. This is useful for field map files. file_paths can be a string:
l.add_input_files('map')
an iterable, for example a list:
l.add_input_files(['map1', 'map2', 'map3'])
To add many files use a pattern eg:
l.add_input_files(pattern="maps/*")
Will use symlinks when avaliable (Linux/UNIX), falls back to copying otherwise.
This class lets you analyse the results after running a line.
returns a list of coordinates specified, and multiply them by the multiplier if not none eg get_trac(‘plt’, [‘X’,’Y’,’element_label’],[0.01,0.01,None])
If all the columns requested are numerical, and new headered data formats are being used then this function will return a numpy array
A bunch object to hold the coordinates for many particles
Note that all values are in SI units, m, rad, eV, s
Object to store a bunch of particles efficiently using numpy. All values are in SI units, m, rad, eV, s Can be use to create a bunch of particles with all there coordinates set to zero (appart from the D coordinate which is set to 1):
my_bunch = Bunch(nparticles=100, ke=1e6, mass=PROTON_MASS, charge=1)
It can also be used to create a bunch from some existing coordinates stored in a numpy array:
my_bunch = Bunch(ke=1e6, mass=PROTON_MASS, charge=1, particles=existing_coords)
Generate a Gaussian bunch in x-xp (Y-T) and in y-yp (Z-P). S and D are set to 0 and 1 respectively. example:
my_bunch = Bunch.gen_kv_x_xp_y_yp(1000, 1e-3, 1e-3, 4, 5, 1e-3, 2e-2, ke=50e6, mass=PROTON_MASS, charge=1)
creates a Gaussian bunch called my_bunch with 1000 particles of the given parameters.
Generate a halo bunch, i.e. an elipse outline in x-xp (Y-T) and in y-yp (Z-P). S and D are set to 0 and 1 respectively. example:
my_bunch = Bunch.gen_halo_x_xp_y_yp(1000, 1e-3, 1e-3, 4, 5, 1e-3, 2e-2, ke=50e6, mass=PROTON_MASS, charge=1)
creates a halo bunch called my_bunch with 1000 particles of the given parameters.
Generate a uniform (KV) bunch, i.e. a filled elipse in x-xp (Y-T) and in y-yp (Z-P). S and D are set to 0 and 1 respectively. example:
my_bunch = Bunch.gen_kv_x_xp_y_yp(1000, 1e-3, 1e-3, 4, 5, 1e-3, 2e-2, ke=50e6, mass=PROTON_MASS, charge=1)
creates a KV bunch called my_bunch with 1000 particles of the given parameters.
Some physical constants.
Algorithm by J. Scott Berg to find the smallest circle that enclose a set of ellipses that lie on the midplane.
Class to allow one to find a circle enclosing a set of ellipses. append((a, b, c)) adds an ellipse to the set of ellipses, where a is the half width, b is the half height, and c is the horizontal center.
radius(z) returns the radius of a circle with horizontal center at z which encloses all the ellipses.
get_circle() returns a sequence (z, r) giving the horizontal center z and radius r of the smallest circle enclosing all of the ellipses
Exceptions raised by pyzgoubi
Module to handle reading and writing of Zgoubi files
All values in energy eV, momentum eV/c, mass ev/c**2, speed c, charge elementry charge
Various useful functions and utilities
make a list of coordinate aranged in a grid in a generalised space eg:
coords_grid(min_c=[0,0],max_c=[1,1],step=[10,10])
will give you 100 points in 2 dimensions between 0 and 1
Given some emittance in horizonal and vertical space
If ncoords <= 1 return points where phase space ellipse crosses the y,y’ and z,z’ axis.
where (a,b) are major and minor radii and phi is the tilt of the major radius w.r.t horizontal axis. A (z,p) distribution is similarly calculated.
Can use these points, returned in coords_YTZP to start tracking at the desired emittance (assuming that the optical functions don’t change with amplitude).
Emittances in both the horizontal and vertical planes may be supplied. Twiss parameters beta and gamma in both places may be determined calling get_twiss_parameters beforehand i.e.:
twissparam = r.get_twiss_parameters()
betayz = [twissparam[0],twissparam[3]]
gammayz = [twissparam[2],twissparam[5]]
Find a closed orbit for the line. can optionally give a list of initial coordinates, init_YTZP, eg: find_closed_orbit(line, init_YTZP=[1.2,2.1,0,0]) otherwise [0,0,0,0] are used.
The line is expected to have a OBJET2 and either a FAISCNL at the end or, using FAISTORE, a MARKER at the end of the cell identified by fai_label. The latter option allows the zgoubi.fai file to contain coordinates throughout the lattice while using just those points with fai_label in this function.
It is recommend to have a REBELOTE, with several laps. The area of the phase space ellipse is approximated from the coordinates from the FAISCNL (or MARKER with fai_label), and the center is used for the new coordinates. Once the relative variation between iterations is less that the tolerance, the function returns the closed orbit coordinates. If a coordinate is close to zero (less than the tolerance) then it is compared absolutely instead of relatively.
Same as find_closed_orbit, but if init_YTZP is unstable, will generate a bunch or particles, with a spread of range_YTZP, and if any of those are stable will do a closed orbit search with that particle:
init_YTZP=[5,0,0,0], range_YTZP=[10,0,0,0], count_YTZP[50,0,0,0]
will vary the Y coordinate from -5 to 15 cm in 50 steps:
range_YTZP=[10,5,0,0], count_YTZP=[10,10,0,0]
would create 100 particles in a grid.
flatten(sequence) -> list
Returns a single, flat list which contains all elements retrieved from the sequence and all recursively contained sub-sequences (iterables).
http://kogs-www.informatik.uni-hamburg.de/~meine/python_tricks
Examples: >>> [1, 2, [3,4], (5,6)] [1, 2, [3, 4], (5, 6)] >>> flatten([[[1,2,3], (42,None)], [4,5], [6], 7, MyVector(8,9,10)]) [1, 2, 3, 42, None, 4, 5, 6, 7, 8, 9, 10]
Calculate tune using FFT. nfourierturns determines the number of passes through the lattice. Can supply set of horizontal and vertical coordinates in coords = [ycoords,zcoords], otherwise
routine will calculate coordinates
Set plot_fourier= True to show Fourier spectrum. Default is False
Generate Gaussian distribution about mean with standard deviation sigma npoints - number of points in distribution required sigma_cutoff - the cutoff factor i.e. points outside sigma*sigma_cutoff eliminated
Seed - Optional parameter to set random seed
Find smallest circle that encloses a set of ellipses centred on the midplane. The ellipses are defined by their horizontal and vertical radii and by their centre along the horizontal axis (a,b,c). The algorithm optimised both the centre of the enclosing circle and its radius.
Ellipse algorithm by J. Scott Berg (BNL)
Reference - J. Scott Berg. ‘Finding the Circular Magnet Aperture which Encloses and Arbitrary number of midplace-centered Beam Ellipses’, Proc. EPAC04, 5-9 July, Lucerne Congress Centre, Lucerne, Switzerland
Calculates the twiss parameters at all points written to zgoubi.plt. 11 particle trajectories are used to calculate the transfer matrix, just as is done in zgoubi. The code mirrors that found in mat1.f, mksa.f. The twiss parameters are first calculated at the end of the cell using either input_twiss_parameters (format [beta_y, alpha_y, gamma_y, beta_z, alpha_z, gamma_z]) or, if this is not supplied, it assumes the cell is periodic and uses get_twiss_parameters to find the boundary condition.
The twiss parameters are then mapped to all points in the magnets using the transfer matrix calculated at each point. The results are stored in list twiss_profiles where each row represents a point in the zgoubi.plt file and has format
[s_coord, label, mu_y, beta_y, alpha_y, gamma_y, mu_z, beta_z, alpha_z, gamma_z]
Requires an OBJET type 5, and a MATRIX element.
Note - This calculation uses trajectories as measured in the local coordinate system of the magnet.
Misalign the set of elements identified by element_indices. The misalignment is Gaussian with sigma, sigma_cutoff and mean among the input parameters. If sigma = 0.0, the misalignment is constant and given by mean
Alternatively, misalignments specified by list misalign_dist
Returns the misailgnment distribution (misalign_dist), unit meters
element_indices can be determined using line.find_elements(element_name)
mean and sigma should be input in meters
Plots multiple sets of data where the X and Y coordinates are each specified in a list of lists. Should also work if a single set of X, Y data is specified or if one X is supplied with multiple Y data points (as long the dimensions of Y equals that of X in all cases).
Required input - data_x_list - Set of X data. May be a single list or a list of lists each with the same dimension. data_y_list - Set of Y data. May be a single list or a list of lists each with the same dimension. Can have a single list of X data and multiple lists of data_y_list. filename - Plot saved to filename.png
Optional parameters - labels = [“Plot Title”,”X axis label”,”Y axis label”] style - The plot style can be supplied as a list or as a single value. If no style is specified, will cycle over a predefined style list. legend - Legends can be supplied as a list of strings or as a single string. legend_location - Default of location of legend box is ‘best’, otherwise can select ‘upper right’ or numerical code for position. See matplotlib documentation. xlim = [lower_value, upper_value] - Limit horizontal axis ylim = [lower_value, upper_value] - Limit vertical axis
# PURPOSE: read an array from a file. Skip empty lines or lines # starting with the comment delimiter (defaulted to ‘#’). # # OUTPUT: a float numpy array # # EXAMPLE: >>> data = readArray(“/home/albert/einstein.dat”) # >>> x = data[:,0] # first column # >>> y = data[:,1] # second column from http://www.scipy.org/Cookbook/InputOutput
Note: please switch to using numpy functions: numpy.savetxt(), numpy.loadtxt() or numpy.genfromtxt()
A scaling FFAG has field given by B/B0 = (r/r0)^k where k is the scaling factor. The Taylor expansion about r0 yields
B(r)/B0 = (k!/(n!(k-n)!)) * (r-r0)^n/r0^n for each multipole n
In DIPOLES the magnetic field (ignoring fringe field factor) is given by
B(r)/B0 = bcoef[n] (r-r0)^n/r0^n for each multipole n
It follows that setting bcoef[n] = (k!/(n!(k-n)!)) will create a DIPOLE magnet whose field approximates the scaling law. Routine returns bcoef up to number of terms specified.
If r0 is scaled by scale_factor, the coefficients are scaled accordingly. If r0 is shifted by d_r0, the routine will update B0 (b0_in) according to the scaling law
In a scaling FFAG, the magnetic field follows the scaling law given by B=B0*(r/r0)^k where r is the radial coordinate, B0 is the field at r=r0 and k is the scaling factor. This def fits a polynomial of a given order to the scaling field at points in the region rmin < r < rmax, where the number of points is determined by the increment step. The coefficients of polynomial are returned.
Required input - b0 - the field at r0 r0 - the reference radius rmin,rmax - region of fit step - step size in (rmin,rmax) at which scaling field is evaluated
Optional input- order - Order of polynomial fit (default 5)
Author - S. Sheehy
Check a list of emittances (emit_list, units Pi m rad) to see if tracking succeeds. Can be used to establish the dynamic aperture. If the elements in emit_list are increasing then will stop tracking once it finds the lowest emittance where a particle is lost. On the other hand, if the elements in emit_list are decreasing then will stop tracking once it reaches the first point where tracking succeeds without loss.
emit_list - List of emittances to check
closedorb_YTZP - Closed orbit coordinates as returned by find_closed_orbit
npass - Number of passes through lattice
D_mom - Momentum factor p/pref
beta_gamma_input - If this is supplied then the emittances supplied in emit_list are assumed to be normalised and a conversion t o geometrical emittance is made. Otherwise the emittances are assumed to be geometrical.
ellipse_coords - If greater than 1 will test uniformly distributed set of coordinates around around phase space ellipse. Otherwise will take single point where phase space cuts y (or z) axis. Can also specify ellipse coords = [n,m] – distribute n points around ellipse but only test point m of these.
twiss_parameters - Can supply twiss parameters in format [[alpha_y,beta_y,gamma_y],[alpha_z,beta_z,gamma_z]], i.e the output of r.get_twiss_parameters(). If not supplied an attempt will be made to calculate the twiss parameters in this function.
plot_data - If True, creates phase space plots at all emittances scanned in both transverse planes.
If a particle is lost, returns index_lost in emit_list where the loss occurs. Otherwise index_lost remains 0.
Plot a list of tunes on a tune diagram with resonance line up to a given order.
At the moment, since the tune diagram covers the range [0,1] just fractional tunes can be shown. The default value of order is 3. Written by Y. Giboudet