This contains the main PyZgoubi classes: Line and Results
The Line object holds a series of elements, to represent an accelerator lattice.
Add an elements to the line. Can also be used to add one line into another.
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.
Check that line has OBJET or MCOBJET at the start, and an END at the end. Gives warnings otherwise. Called by run() if in debug mode.
clean up temp directories
Iterator for elements in line, including sub lines.
Returns all the positions of element in line, treats sub lines as linear list
Enable full tracking on magnetic elements. This works by setting IL=2 for any element with an IL parameter. use line.full_tracking(False) to disable tracking drift_to_multi=True replaces drifts with weak (B_1 = 1e-99) multipoles
Find the OBJET element
Insert elements into the line before position given by index, treats sub lines as linear list
Generate the zgoubi.dat file, and return it as a string
Add a elements to the start of the line
Remove element at index, treats sub lines as linear list
removes any REBELOTE elements from the line
Replace an element in the line. setting select_index to n will replace the nth occurence of that item. If select index is not set, the first occurence is replaced.
Note: elementold must be the same instance as an element in the Line
Run zgoubi on line. If xterm is true, stop after running zgoubi, and open an xterm for the user in the tmp dir. From here zpop can be run. Returns a Results object
Track a bunch through a Line, and return the bunch. This function will uses the OBJET_bunch object, and so need needs a Line that does not already have a OBJET. If binary is true then particles are sent to zgoubi in binary (needs a version of zgoubi that supports this)
This function should be used identically to the track_bunch function, apart from the addition of the n_threads argument. This will split the bunch into several slices and run them simultaneously. Set n_threads to the number of CPU cores that you have. max_particle can be set to limit how many particles are sent at a time.
This class lets you analyse the results after running a line.
It is created automatically and returned by Line.run()
return binary fai file as string
return file handle for binary fai file
return binary plt file as string
return file handle for binary plt file
read whole file. for each chunk, check if it is an element with bounds, if yes check if it is inside record crashes, by particle id min_bounds and max_bounds should be of the form {‘label’: bound, ‘label’: bound} returns numpy arrays note: lost particles dont crash
clean up temp directory
return dat file as string
return file handle for dat file
return fai file as string
return file handle for fai file
Read all the data out of the file. Set file to plt or fai if using the new headered data formats will return a numpy array with named columns, otherwise returns a list of dicts. each dict has the particle coordinates at a point
“Get back a bunch object from the fai file. It is recommended that you put a MARKER before the last FAISCNL, and pass its label as end_label, so that only the bunch at the final position will be returned. All but the final lap is ignored automatically. Optionally the an old_bunch can be passed to the function, its mass and charge will be copyed to the new bunch.
Get the max and min positions of a certain coordinate, in a certain element
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
Returns a transfer matrix of the line in (MKSA units). Needs a beam line is an OBJET type 5, and a MATRIX element.
Returns a tuple (NU_Y, NU_Z) of the tunes. Needs a beam line is an OBJET type 5, and a MATRIX element.
The results are stored a in structured numpy array containing the following elements
#Outputs beta_y, alpha_y, gamma_y: horizontal Twiss parameters beta_y, alpha_y, gamma_y: vertical Twiss parameters disp_y, disp_py : horizontal dispersion and dispersion-prime disp_z, disp_pz : vertical dispersion and dispersion-prime
Needs a beam line is an OBJET type 5, and a MATRIX element.
return impdev file as string
return file handle for impdev file
check if particle exceeded bounds in this element. pass bounds in zgoubi’s default unit for that coordinate
Returns False if no losses, otherwise returns a summery of losses
loss = res.loss_summary(file='plt')
#or
all = res.get_all('plt')
loss = res.loss_summary(all) # if you already have got the coordinates
Parse data from output of MATRIX element
Used by get_tune(), get_transfer_matrix() and get_twiss_parameters()
return plt file as string
return file handle for plt file
return res file as string
return file handle for res file
Checks that zgoubi completed
save binary fai file to path
save binary plt file to path
save dat file to path
save fai file to path
save impdev file to path
save plt file to path
save res file to path
save spn file to path
show the particle info, a good check of energies, mass etc
return spn file as string
return file handle for spn file
Return true if end of REBELOTE procedure reported Needs a REBELOTE element.
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)
Check that the bunch is not empty, and contains finite values
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 Gaussian bunch in transverse and longitudinal phase space emit_y, emit_z : horizontal and vertical plane geometric emittance (1 sigma) beta_y, beta_z : horizontal and vertical betatron function alpha_y, alpha_z : horizontal and vertical alpha mom_spead : sigma of momentum spread (percentage) bunch_length : sigma of bunch length (metres) disp, disp_prime : dispersion and dispersion prime (D’) in horizontal plane
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 surface of a 4D hypershere in x-xp-y-yp (Y-T-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.
Generate a waterbag bunch, i.e. a filled hypersphere in x-xp-y-yp (Y-T-Z-P). S and D are set to 0 and 1 respectively. example:
my_bunch = Bunch.gen_waterbag_x_xp_y_yp(1000, 1e-3, 1e-3, 4, 5, 1e-3, 2e-2, ke=50e6, mass=PROTON_MASS, charge=1)
creates a waterbag bunch called my_bunch with 1000 particles of the given parameters.
Get bunch kinetic energy
Get bunch rigidity
Returns the center of the bunch in each dimension Y,T,Z,P,S,D
return emittance h and v in m rad. Uses the bunch full width, so should only be used for a hard edge distribution
return emittance h and v in m rad. Uses the RMS quantities
Returns the minimum rigidity of the bunch
Returns the twiss valuse Beta_h, Alpha_h, Beta_v, Alpha_v, calculated from bunch width extent
Returns the rms twiss valuse Beta_h, Alpha_h, Beta_v, Alpha_v, calculated from bunch rms width
Returns the width of the bunch in each dimension Y,T,Z,P,S,D
Returns the rms width of the bunch in each dimension Y,T,Z,P,S,D
Returns the numpy array that holds the coordinates, as a structured numpy array
Plot a bunch, if no file name give plot is displayed on screen. lims can be used to force axis limits eg [lY,lT,lZ,lP,lX,lD] would plot limit plot from -lY to +lY in Y, etc. Additional bunches can be passed, as add_bunch, to overlay onto the same plot. fmt can be a list of formats in matplotlib style, eg [‘rx’, ‘bo’]
Returns the numpy array that holds the coordinates, as a 2d numpy array and a list of comumn names. It is possible that the column names may change or reorder, so if you use this you may want to protect against it with something like
assert(pbunch.raw_particles()[1] == ["D","Y","T","Z","P","S","TOF","X"])
Read in a bunch from a file. assumes columns are Y, T, Z, P, X, D, separated by white space:: my_bunch = Bunch.read_YTZPSD(“mybunch.dat”, ke=1e9, mass=ELECTRON_MASS, charge=-1)
Set bunch kinetic energy
Set bunch rigidity
Split a bunch into n_slices smaller bunches, or more if they would have too many particles in.
Output a bunch, compatible with read_YTZPSD
Common functions that can be used anywhere in pyzgoubi
Make a directory. Make parents if needed. Don’t raise error if it already exists
From: http://stackoverflow.com/questions/600268/mkdir-p-functionality-in-python code by tzotzioy http://stackoverflow.com/users/6899/
Pass either a filename or file handle like object. Returns a file like object
If mode is ‘w’ or ‘a’ then can set mkdir=True to make sure the directory is created if needed.
Helper function to make the twiss param array used as input to get_twiss_profiles() or get_cell_properties_nonperiodic()
tp = twiss_param_array(beta_y=2, alpha_y=-1 ...) all unspecified will be left at 0
Some physical constants.
Atomic mass unit (eV)
Electron anomalous magnetic moment
Electron charge (C)
Electron mass (eV)
Electron lifetime
Muon anomalous magnetic moment
Muon charge (C)
Muon mass (eV)
Muon lifetime
Pion charge (C)
Pion mass (eV)
Pion lifetime
Proton anomalous magnetic moment
Proton charge (C)
Proton mass (eV)
Proton lifetime
Speed of light (m/s)
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
Add an ellipse e to the set of ellipses
Return (z,r), where z is the center and r is the radius of the circle with the smallest radius that encloses all the ellipses
Get the radius of a circle, centered at z, enclosing all of the ellipses
Add the interval ivl=(min, max) to the interval list l
Exceptions raised by PyZgoubi
Can’t parse file format
Line not sufficient for this function
File contains no records
No particle track
Old file format passed to function that expects new format
Zgoubi hit a fatal error
Subclass of numpy.ndarray to add an info field
Quick function for printing parameter tables, eg:
print gcp.cell_properties_table(data_arc, [“KE”, “stable”, “Y”, “T”, “NU_Y”, “NU_Z”])
Get the closed orbits and basic properties of a periodic cell.
cell: A PyZgoubi Line object containing the beamline elements min_ke, max_ke, ke_steps: kinetic energy in eV. For a single step just set min_ke particle: “p”, “e”, “mu-”, “mu+”, or “Bismuth” tol: tolerance when finding closed orbit stop_at_first_unstable: If an energy is found to be unstable, give up closed_orbit_range, closed_orbit_range_count, closed_orbit_init_YTZ: see zgoubi.utils.find_closed_orbit_range() reuse_co_coords: use closed orbit from previous energy to start search for next energy closed_orbit_debug: output debuging information full_tracking=True is required in order get minimum and maximum magnetic fields along orbit
returns orbit_data, an array with ke_steps elements, with the following data KE : particle KE in eV stable, found_co, stable_tm_YT, stable_tm_ZP : boolean stability flags Y, T, Z, P : closed orbit at start of cell, in cm and mrad BETA_Y, BETA_Z, ALPHA_Y, ALPHA_Z, GAMMA_Y, GAMMA_Z, DISP_Y, DISP_Z, DISP_PY, DISP_PZ, NU_Y, NU_Z: periodic twiss functions tof, S: time of flight and path length along closed orbit matrix, matrix_trace_YT, matrix_trace_ZP: transfer matrix and traces MAX_BY, MIN_BY, MAX_BZ, MIN_BZ: minimum and maximum fields seen along closed orbit (requires full_tracking=True)
Get the basic properties of a non-periodic cell.
Works similarly to get_cell_properties(), but rather than finding a periodic solution for closed orbit and twiss parameters, takes them as input: init_YTZP: list of the Y, T, Z and P at the start of the cell init_twiss: initial twiss params, use twiss_param_array() to create
cell: A PyZgoubi Line object containing the beamline elements min_ke, max_ke, ke_steps: kinetic energy in eV. For a single step just set min_ke particle: “p”, “e”, “mu-”, “mu+”, or “Bismuth” full_tracking=True is required in order get minimum and maximum magnetic fields along orbit
returns orbit_data, an array with ke_steps elements, with the following data KE : particle KE in eV stable : boolean stability flags, always true, but useful for controlling which orbits are processed by other functions Y0, T0, Z0, P0 : position and angles at start of cell, in cm and mrad Y, T, Z, P : position and angles at end of cell, in cm and mrad BETA_Y0, BETA_Z0, ALPHA_Y0, ALPHA_Z0, GAMMA_Y0, GAMMA_Z0, DISP_Y0, DISP_Z0, DISP_PY0, DISP_PZ0, NU_Y0, NU_Z0: twiss functions at start of cell BETA_Y, BETA_Z, ALPHA_Y, ALPHA_Z, GAMMA_Y, GAMMA_Z, DISP_Y, DISP_Z, DISP_PY, DISP_PZ, NU_Y, NU_Z: twiss functions at end of cell tof, S: time of flight and path length along closed orbit matrix, matrix_trace_YT, matrix_trace_ZP: transfer matrix and traces MAX_BY, MIN_BY, MAX_BZ, MIN_BZ: minimum and maximum fields seen along closed orbit (requires full_tracking=True)
Get tracks along the closed orbit for values in data cell: periodic cell data: the data structure return from get_cell_properties() particle: see get_cell_properties() full_tracking: record all steps in the magnet
tracks are added in the data structures ftrack and ptrack fields
Get Dynamic Aperture.
cell: the cell to run data: data structure containing the closed orbit to be used initial conditions, see get_cell_properties()
tol: tolerance 0.01 give result to 1% precision npass: number of passes through the cell quick_mode: “+y+z” 1 particle at [dY,0,dZ,0]
“+t+p” 1 particle at [0,dT,0,dP] False:16 combinations of plus and minus offsets are used
nangles: number of realspace angles to use. 1 will give the 45degree DA, i.e. with equal horizontal and vertical excitation, 2 will give pure horizontal (0 deg) and pure vertical (90deg), for other numbers will use numpy.linspace(0,pi/2,nangles). island_avoid: when an unstable amplitude is found take a small step up, to see if it just a small island. set to zero to disable debug_log: file name to write debug information to start: starting emittance
From the starting emittance a search for the stability boundary is made.
Track particle for npass turns
Run a particle through a cell from a give starting point and return track from fai and plt files.
This is mostly used by other functions in this module, but can be useful for debugging a lattice
Get the periodic Twiss (Courant and Snyder) parameters for the cell. Tracks a bunch of particles containing pairs offset in each plane.
Must pass in data returned from get_cell_properties() to give initial conditions. The profile is added to the ‘twiss_profile’ key in the data structure. If full_tracking is enabled, the Twiss parameters are every integration step are recorded, otherwise parameters are recorded at the end of each magnetic element.
Look up a particle by name and return the zgoubi PARTICUL, mass and charge sign
Produce a set of standard plots from the data structure
data: the data structure returned by get_cell_properties()
Plot particle track through cell, starting from in the closed orbits stored in data.
output_file: file to save plot show: display plot to screen plot_unstable: show tracks for initial coordinates marked as unstable in data structure draw_field_midplane: show the magnet field on the magnet midplane sector_width: width to draw sector magnets aspect: “equal”, use same scale for x and y, “auto” stretch to fit style: style for plotting, see lab_plot.LabPlot draw_tracks: draw particle tracks min_y, max_y, y_steps, angle: starting coordinates for test particles for sampling the midplane field plot_extents: list of extents [left, right, bottom, top]
Make phase space plots showing dynamic aperture
Plots radial, longitudinal and 2d mid-plane magnetic field for each element, using high rigidity rays to sample field.
min_y, max_y, y_steps, angle set the starting positions and angle for the test rays If angle=None, then 0 will be used for rectangular magnets, and AT/2 for sector magnets
output_prefix_rad, output_prefix_long, output_prefix_long should contain a ‘%s’ that is replaced with the element label, eg “results/mag_field_rad_%s.pdf” output_prefix_long and output_prefix_2d will be skipped if they are left blank The longitudinal profile is along the straight lines of the test rays, not the closed orbit.
For the radial plot method can be: interp: points along a radial line are interpolated from a 2d profile max: use the max field seen by the ray
extra_test_range factor to extend the range (min_y, max_y) for test particles, to improve interpolation edges
In 2d plots the trick of making rectangular magnets by using large radii is recognised for radii over 1e6cm
Plot the Twiss profiles found with get_twiss_params()
fields: list of fields to plot, beta_y, beta_z, alpha_y, alpha_z, gamma_y, gamma_z, disp_y, disp_z, disp_py, disp_pz
Returns a 1d transverse horizontal (radial) profile of the magnet. At mid point, or ACN for DIPOLES. Used by plot_element_fields()
returns field, ys
Get a 2d interpolated field profile. Used by plot_element_fields
returns a y_steps by y_steps grid of interpolated field values, and the extents of that grid:
returns int_field, (xmin,xmax,ymin,ymax)
Internal function used by profile1d() and profile2d()
Module to handle reading and writing of Zgoubi files
Read header from a file and determine formating. Returns a dict that describes the file
Replace all occurrences of ‘old’ with ‘new’ in ‘l’
Read a zgoubi output file. Return a numpy array with named column headers. The format is automatically worked out from the header information.
Read 1 record from a fortran file
Run define file on set of files, and output a code block that can be put into the to of this file to save rerunning define_file at program runtime
Write a record, adds record length to start and end
Module for drawing plots of a lattice in lab coordinates
A plotter for beam lines and tracks.
Add tracks from plt or fai files.
Draw the plot in memory, then use show() to display to screen or save() to save to file.
plot_extents: list of extents [left, right, bottom, top]
Style is a nested dictionary, that updates the defaults, eg:
style = {"track":{"color":"g"}}
Elements to be styled include “track”, “magnet_outline”, “element_outline”, “reference”. Each can take a “color”, “linestyle” and “linewidth”, using matplotlib notations.
All values in energy eV, momentum eV/c, mass ev/c**2, speed c, charge elementry charge
Convert speed to lorentz factor
Convert speed to kinetic energy
Convert lorentz factor to speed
Convert lorentz factor to kinetic energy
Convert kinetic energy to speed
Convert kinetic energy to lorentz factor
Convert kinetic energy to momentum
Convert kinetic energy to magnetic rigitity
Convert kinetic energy to total energy
Convert momentum to kinetic energy
Convert momentum to magnetic rigitity
Convert momentum to total energy
Convert magnetic rigitity to kinetic energy
Convert magnetic rigitity to momentum
Convert total energy to kinetic energy
Convert total energy to momentum
can’t handle noise. finds closest and furthest point from centre, assumes they are a and b
Given a phase slip and Lorentz gamma, find momentum compaction factor and transition gamma
Calculate the phase advance (mu_y,mu_z) from a transfer matrix. Either use Results.get_transfer_matrix() or calc_transfer_matrix() to get matrix.
Calculate phase slip at relative momentum D (D=1 by default). Reference TOF must be supplied (tof_ref). Optionally set del_p: momentum shift tol_co: adjust tolerance of find_closed_orbit calculation. init_YTZP: set initial guess in find_closed_orbit search.
Track a bunch generated with OBJET5 through a line, and pass the start and end bunch to this function to calculate the twiss matrix. No unit conversion is done, so units match the passed bunch. NOT COMPLETE: only use cells in top 4 rows
Calculate the twiss parameters (beta_y, alpha_y, gamma_y, beta_z, alpha_z, gamma_z) from a transfer matrix. Either use Results.get_transfer_matrix() or calc_transfer_matrix() to get matrix.
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[5]]
gammayz = [twissparam[2],twissparam[7]]
Note coords are returned in ‘Zgoubi units’, i.e. cm and mrad
find centre by using mean of coords, this assumes that point are evenly distributed around ellipse
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.
record_fname is used to record search details to a file that can be used with plot_find_closed_orbit(). extra_iterations allows some extra iterations to improve accuracy, even after tolerance reached
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.
Find all occurences of list_element in list. Return the indices.
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
get things formated like key=value from the commandline. returns the requested value
get things formated like key=value from the commandline. returns true for: yes, on, true, 1 returns false for: no, off, false, 0 case insensitve
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 given below) or, if this is not supplied, it assumes the cell is periodic and uses get_twiss_parameters to find the boundary condition.
The optional input_twiss_parameters should be supplied in the form of a numpy structured array containing the following information at the start of the cell [beta_y, alpha_y, gamma_y, disp_y, disp_py, beta_z, alpha_z, gamma_z, disp_z, disp_pz] This is the same format as the output of get_twiss_parameters allowing one to conveniently obtain twiss parameters at the end of one section of beam line and apply them as starting conditions for the next section, i.e.
r = line1.run() twiss_line1 = r.get_twiss_parameters() twiss_profiles = get_twiss_profiles(line2, input_twiss_parameters = twiss_line1)
Dispersion and dispersion-prime are also calculated if calc_dispersion is True (default).
interpolate: If a azimuthal coordinate X exists, interpolate onto the reference X. X is azimuthal angle in the case of radial elements.
The results are stored a in structured numpy array called twiss_profiles containing the following elements at every point tracked in the magnets (i.e. a point in zgoubi.plt) [s_coord, label, mu_y, beta_y, alpha_y, gamma_y, disp_y, disp_py, mu_z, beta_z, alpha_z, gamma_z, disp_z, disp_p]
The results are also stored in file specified by file_result (optional)
Requires an OBJET type 5, and a MATRIX element.
Note - This calculation uses trajectories as measured in the local coordinate system of the magnet.
ke in eV, mass in eV/c^2
ke in eV, mass in eV/c^2
ke in eV, mass in eV/c^2 gives result in kGauss cm, which is handy for OBJET
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
mom in eV/c gives results in eV
mom in eV/c gives result in kGauss cm, which is handy for OBJET
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 majorLocator = horizontal axis tick marks defined as multiple of majorLocator
When the closed orbit search fails it can be useful to see what happened. find_closed_orbit() and find_closed_orbit_range() can take an optional argument record_fname. This causes them to write a log file, which can be read by this function and plotted.
closed_orbit = find_closed_orbit(tline, init_YTZP=[1,20,0,0], record_fname="res/closedorbit.log")
plot_find_closed_orbit(data_fname="res/closedorbit.log", outfile="res/closedorbit.png")
This will show you the succession of guesses. If using plot_find_closed_orbit_range() it will also show a plot of all the coordinates that survive being tracked through the lattice.
From the plots it may be clear whether the lattice is unstable, or if the convergence is just too slow. In the latter case it may help to repeat the cell using REBELOTE, or to increase the number of iterations allowed.
# 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_h - List of emittances in horizontal plane to check emit_list_v - List of emittances in vertical plane 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 to geometrical emittance is made. Otherwise the emittances supplied are assumed to be geometrical.
ellipse_coords - If greater than 1 will test uniformly distributed set of coordinates around around phase space ellipse. The points are generated by emittance_to_coords. ellipse_coords = 1 (or -1) will take single point where phase space cuts either the y, t, z or p axes (decided by coord_pick). ellipse_coords = -1 selects the negative point in the phase space ellipse. The default settings ellipse_coords = 1, coord_pick = 0 selects the point where the phase space ellipse cuts the y, z axis to the outside of the closed orbit point.
coord_pick - Choose one of ellipse coords, i.e. ellipse_coords[coord_pick] is selected. If ellipse_coords equal 1 or -1, then coords_pick = 0 selects the points where the phase space ellipse crosses the Y, Z axes. If coord_pick = 1 then selects the point where the phase space ellipse crosses the T, P axes. If coord_pick = None, will select all points generated by emittance_to_coords.
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 [beta_y, alpha_y, gamma_y, disp_y, disp_py, beta_z, alpha_z, gamma_z, disp_z, disp_pz], 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.
Returns [index_lost, coord_index], fourier_tune_emit, coords_YTZP_ini] where index_lost is the index in emit_list where particle is lost (0 if no loss). coord_index indicates at which coord in ellipse_coords the particle is lost. fourier_tune_emit is a list of tunes found at each emittance using Fourier analysis. coords_YTZP_ini_list is a list of all initial coordinates tested.
If a particle is lost, returns index_lost in emit_list where the loss occurs. Otherwise index_lost remains 0.
shows a file. mode can be: cat - dumps to screen less - scrollable view win - scollable new win (less in an xterm)
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
Returns a new line where every element has a unique label.
Checks for repeating labels and appends a count to each that repeats. Ignores elements with blank labels.
Note: New line contains deep copies of original elements.