from .transfermatrix import TransferMatrix
from .elements import Collimator
from .elements import TMElement
from .settings import PLOTBACKEND
import sys
import os
import copy
import subprocess
import contextlib
import logging
import numpy as np
import scipy.constants as const
from scipy.spatial.transform import Rotation
import matplotlib.patches as patches
import matplotlib.lines as lines
logger = logging.getLogger(__name__)
@contextlib.contextmanager
def cd(newdir):
prevdir = os.getcwd()
os.chdir(os.path.expanduser(newdir))
try:
yield
finally:
os.chdir(prevdir)
__all__ = ['GICOSYElement', 'GICOSYDipole', 'GICOSYDrift', 'GICOSYQuadrupole', 'setGICOSYPath', 'writeFringeFields', 'setFringeField']
_gicosy_path = None
[docs]def setGICOSYPath(path):
"""Define the path to the GICOSY installation directory.
The GICOSYElement and the derived classes require access to GICOSY
program for producing the elements matrices. This function is used to set
the path to the GICOSY program manually. The access to GICOSY can
be defined with this function or by default the PIOL searches for
the presence of either a GICOSYPATH environment variable or a
folder with 'gicosy' in the PATH. If neither can be found, an
exception will be raised.
Parameters
----------
path : str
Path to GICOSY installation directory.
"""
global _gicosy_path
_gicosy_path = path
def readGICOSYPath():
"""Reads the environment variables GICOSYPATH and PATH and defines GICOSY installation directory.
The GICOSYElement and the derived classes require access to GICOSY
program for producing the elements matrices. This function
searches for the presence of either a GICOSYPATH environment
variable or a folder with 'gicosy' in the PATH. If neither can be
found, an exception will be raised.
This function is automatically called by the GICOSYElement as needed.
"""
global _gicosy_path
if _gicosy_path is None:
if 'GICOSYPATH' in os.environ:
_gicosy_path = os.environ['GICOSYPATH']
logger.info('Found environment variable GICOSYPATH: {:s}'.format(_gicosy_path))
else:
path = os.environ['PATH'].split(';')
for p in path:
if 'gicosy' in p.lower():
_gicosy_path = p
logger.info('Found gicosy location in PATH: {}'.format(_gicosy_path))
break
else:
logger.info('\n'.join(path))
raise Exception('No GICOSYPATH environment variable set, nor location found in PATH')
[docs]class GICOSYElement(TMElement):
"""Base class for all transfer matrix elements using gicosy.
"""
def __init__(self, refple, name='', description=''):
super().__init__(refple, name, description)
self.collimator = None
self.parts = 1
self.tm_ffentrance = None
self.tm_part = None
self.tm_ffexit = None
self.tm_full = None
self.part_ds = 0.0
self.part_dr = np.zeros(3)
self.part_dx = np.zeros(3)
self.part_dy = np.zeros(3)
[docs] def transfer(self, driver):
"""
Transfers the particles saved in the driver through the element.
Parameters
----------
driver : Driver
Driver of the system
"""
driver.onElementEnter(self)
# Getting particles as a copy.
pleset = driver.getParticlesAsCopy()
# Creating a new dictionary but just copying the values as references.
coords = copy.copy(pleset)
# Building gicosy coordinates G and D and adding to the newly created dict.
refple = self.parameters["refple"]
coords['g'] = pleset['m']/pleset['q']*(refple.charge/refple.mass) - 1.0
coords['d'] = pleset['K']/pleset['q']*(refple.charge/refple.energy) - 1.0
# Entrance fringe field if needed
if self.parameters['ff_entrance']:
pleset = driver.getParticlesAsCopy()
pleset['g'] = pleset['m']/pleset['q']*(refple.charge/refple.mass) - 1.0
pleset['d'] = pleset['K']/pleset['q']*(refple.charge/refple.energy) - 1.0
self.tm_ffentrance.apply(pleset)
del pleset['g']
del pleset['d']
driver.registerTransfer(pleset, 0.0)
# Running collimator before actual element transfer matrices
# and then after each transfer matrix.
if self.collimator:
self.collimator.transfer(driver)
# Looping over parts
for i in range(self.parts):
pleset = driver.getParticlesAsCopy()
pleset['g'] = pleset['m']/pleset['q']*(refple.charge/refple.mass) - 1.0
pleset['d'] = pleset['K']/pleset['q']*(refple.charge/refple.energy) - 1.0
self.tm_part.apply(pleset)
del pleset['g']
del pleset['d']
driver.registerTransfer(pleset, self.part_ds, self.part_dr, self.part_dx, self.part_dy)
# Collimating after each part
if self.collimator:
self.collimator.transfer(driver)
# Exit fringe field if needed
if self.parameters['ff_exit']:
pleset = driver.getParticlesAsCopy()
pleset['g'] = pleset['m']/pleset['q']*(refple.charge/refple.mass) - 1.0
pleset['d'] = pleset['K']/pleset['q']*(refple.charge/refple.energy) - 1.0
self.tm_ffexit.apply(pleset)
del pleset['g']
del pleset['d']
driver.registerTransfer(pleset, 0.0)
# Running collimator after the fringe field.
if self.collimator:
self.collimator.transfer(driver)
driver.onElementExit(self)
[docs] def getTransferMatrix(self):
"""
Returns full transfer matrix of the system.
If full transfer matrix is found already then returning it.
Otherwise calculating it.
"""
#if self.tm_full != None:
# return self.tm_full
matrices = []
matrices.append(self.tm_ffentrance)
matrices.extend([self.tm_part] * self.parts)
matrices.append(self.tm_ffexit)
matrices = list(reversed([matrix for matrix in matrices if matrix is not None]))
total = matrices[0]
for matrix in matrices[1:]:
total = total.multiply( matrix, order=max(total.order,matrix.order))
return total
[docs]class GICOSYDipole(GICOSYElement):
"""Dipole element
Interact with a magnetic or electric dipole to bend a beam with the given radius
and deflection angle.
Parameters
----------
refple : ReferenceParticle
Reference particle for the element creation.
ftype : {'electric', 'magnetic'}
String representing the type of dipole.
rho0 : float
Radius of curvature in m.
phi0 : float
Angle of deflection in degrees.
G0 : float
Half gap of dipole in m.
directionAngle : float
Defines the direction of B field in degrees in the dipole in
the optical axis xyz-coordinate system at the entrance. 0 deg
corresponds the case that B field is towards positive
y-direction and positive particles are to be bend towards
negative x-axis (by Lorentz force), 90 deg <-> B field towards
positive x-axis and positive particles are to be bend upwords,
180 deg <-> B-field to neg. y-axis and positive particles are
to be bend towards positive x-axis.
nvalues : list of floats, optional
Field indices of the multipole expansion.
eps_entrance : float, optional
Entrance angle of the field boundary in degree (default zero).
curv_entrance : float, optional
Entrance field curvature relative to rho0 (default zero).
eps_exit : float, optional
Exit angle of the field boundary in degree (default zero).
curv_exit : float, optional
Exit field curvature relative to rho0 (default zero).
ff_entrance : int, optional
Fringe fields table number to use for entrance. Use zero for no fringe field.
ff_exit : int, optional
Fringe fields table number to use for exit. Use zero for no fringe field.
order : int, optional
Order of the transfermatrix. Defaults to 1.
maxSAdvance : float, optional
Maximal advance along the s axis in m. Defaults to 0.2
collimator : float or tuple, optional
Collimator parameters if one is used.
name : string, optional
Name given to the element.
description : string, optional
Description of the element, currently not used.
Attributes
----------
collimator : Collimator or None
Collimator used in the element.
parts : int
Number of times the transfer matrix is to be used.
Notes
-----
The most parameters are saved in the parameters dictionary. A warning will be issued in the logger
if the collimator is larger than it should be.
"""
def __init__(self, refple, ftype, rho0, phi0, G0, directionAngle=0, nvalues=None,
eps_entrance=0.0, curv_entrance=0.0, eps_exit=0.0, curv_exit=0.0,
ff_entrance=1, ff_exit=1, order=1, maxSAdvance=0.2, collimator=None,
name='Dipole', description=''):
super().__init__(refple, name=name, description=description)
if ftype not in ['electric','magnetic']:
raise ValueError('ftype must be \'electric\' or \'magnetic\'')
# Storing all parameters needed for calculation of matrices
self.parameters['ftype'] = ftype
self.parameters['rho0'] = rho0
self.parameters['phi0'] = phi0
self.parameters['G0'] = G0
self.parameters['directionAngle'] = directionAngle
self.parameters['nvalues'] = nvalues
self.parameters['eps_entrance'] = eps_entrance
self.parameters['curv_entrance'] = curv_entrance
self.parameters['eps_exit'] = eps_exit
self.parameters['curv_exit'] = curv_exit
self.parameters['order'] = order
self.parameters['ff_entrance'] = ff_entrance
self.parameters['ff_exit'] = ff_exit
self.parameters['maxSAdvance'] = maxSAdvance
self.parameters['displayfield'] = True
self.collimator = collimator if collimator is None else Collimator(collimator)
self.parameterModified = True
self.rebuild()
def getCalibration(self, massrange=[1, 300], energyrange=[1000, 60000], energysteps=10):
mass = np.arange(min(massrange), max(massrange)+1, 1.0)
energy = np.linspace(min(energyrange), max(energyrange), energysteps)
mass, energy = np.meshgrid(mass, energy)
refple = self.parameters['refple']
kind = self.parameters['ftype']
if kind == 'magnetic':
field = np.sqrt(2*energy*mass*const.u/const.e)/ \
(refple.charge*self.parameters['rho0'])
else:
V0 = energy
r0 = self.parameters['rho0']
r1 = r0 - self.parameters['G0']
r2 = r0 + self.parameters['G0']
field = V0*np.log(r2/r1)
return mass, energy, field
[docs] def field(self, kind):
if kind == self.parameters['ftype']:
refple = self.parameters['refple']
if kind == 'magnetic':
field = np.sqrt(2*refple.energy*refple.mass*const.u/const.e)/ \
(refple.charge*self.parameters['rho0'])
else:
V0 = refple.energy
r0 = self.parameters['rho0']
r1 = r0 - self.parameters['G0']
r2 = r0 + self.parameters['G0']
field = V0*np.log(r2/r1)
return field
else:
return 0.0
def setVoltage(self, field):
if self.parameters['ftype'] != 'electric':
raise TypeError
r0 = self.parameters['rho0']
r1 = r0 - self.parameters['G0']
r2 = r0 + self.parameters['G0']
V0 = field / np.log(r2/r1)
self.parameters['refple'].energy = V0
def length(self):
return self.parameters['rho0'] * np.deg2rad(self.parameters['phi0'])
[docs] def height(self):
return self.parameters['G0'] * 2
[docs] def opticalAxisExit(self):
"""Returns the exiting optical axis.
Returns the coordinate transformation from the entrance plane
to the exit plane. The returned value is a tuple (r,x,y),
where r (vector) is a translation from entrance optical axis
origin to exit optical axis origin; x and y are unit vectors
defining the exit plane coordinate system. The z-axis is the
cross product of x and y, and points towards positive optical
axis.
"""
phi0rad = np.deg2rad(self.parameters['phi0'])
dirRad = np.deg2rad(self.parameters['directionAngle'])
# Optical axis direction at entrance (this is definition, not
# in absolute coordinates).
# Rotation around entrance y-axis
phirot = Rotation.from_rotvec([0.0,-phi0rad,0.0])
# Rotation around entrance z-axis
dirrot = Rotation.from_rotvec([0.0,0.0,dirRad])
# Vector from entrance to exit
dR = dirrot.apply( self.parameters['rho0']*np.array( [-1+np.cos(phi0rad), 0, np.sin(phi0rad) ] ) )
x0 = np.array([1.0,0.0,0.0])
print("### x0 = ", x0 )
y0 = np.array([0.0,1.0,0.0])
z0 = np.array([0.0,0.0,1.0])
z1 = dirrot.apply( phirot.apply( z0 ) )
print("### z1 = ", z1 )
dirrot2 = Rotation.from_rotvec( -dirRad*z1 )
u0 = phirot.apply( x0 )
uu0 = dirrot.apply( u0 )
x1 = dirrot2.apply( uu0 )
v0 = phirot.apply( y0 )
vv0 = dirrot.apply( v0 )
y1 = dirrot2.apply( vv0 )
return np.round(dR,decimals=9), np.round(x1-x0,decimals=9), np.round(y1-y0,decimals=9)
[docs] def plot(self, ax, view, s):
length = self.length()
try:
if 'y' in view:
height = self.collimator.height('y')
else:
height = self.collimator.height('x')
except:
height = max(self.collimator.height())
if not np.isfinite(height):
height = self.height()
try:
thickness = self.parameters['thickness']
except:
thickness = height * 0.1
fill = 'palegoldenrod'
hatch_fill = 'royalblue'
color = 'slategray'
if PLOTBACKEND == 'matplotlib':
from matplotlib import patches
patch = patches.Rectangle((s, -height/2), length, height, fill=True, fc=fill, ec=color, zorder=0)
draw_patches = [patch]
hatched_1 = patches.Rectangle((s, height/2), length, thickness, fill=True, hatch='////', fc=hatch_fill, ec=color, zorder=0)
hatched_2 = patches.Rectangle((s, -height/2-thickness), length, thickness, fill=True, hatch='////', fc=hatch_fill, ec=color, zorder=0)
if ('x' in view and self.parameters['ftype'] == 'electric') or ('y' in view and self.parameters['ftype'] == 'magnetic'):
draw_patches.extend([hatched_1, hatched_2])
try:
displaytext = self.parameters['display']
except:
displaytext = False
for patch in draw_patches:
ax.add_patch(patch)
if displaytext:
for t in [((s+length/2, height/2+thickness+0.001, self.name), {'ha': 'center'})]:
ax.text(*t[0], **t[1])
[docs] def rebuild(self):
if not self.collimator:
if self.parameters['ftype'] == 'magnetic':
self.collimator = Collimator((self.parameters['G0'], self.parameters['G0']))
else:
self.collimator = Collimator((self.parameters['G0'], self.parameters['G0']))
else:
if self.parameters['ftype'] == 'magnetic':
try:
aperture = self.collimator.parameters['size'][1]
except:
aperture = self.collimator.parameters['size']
if aperture > self.parameters['G0']:
logger.warning('Magnetic dipole aperture is larger in the y direction than the pole gap! Results can be unphysical!')
else:
try:
aperture = self.collimator.parameters['size'][0]
except:
aperture = self.collimator.parameters['size']
if aperture > self.parameters['G0']:
logger.warning('Electric dipole aperture is larger in the x direction than the pole gap! Results can be unphysical!')
rho = np.abs(self.parameters['rho0'])
self.totalDeltaS = rho*np.deg2rad(self.parameters['phi0'])
self.parts = int(self.totalDeltaS / self.parameters['maxSAdvance'] + 1)
commands = []
# Need to have a drift length in order to generate the
# fringe field. This matrix is not used nor passed for a user.
commands.append('DRIFT LENGTH 0.1')
commands.append('FRINGING FIELD {} {:g} {:g}'.format(self.parameters['ff_entrance'], self.parameters['eps_entrance'], self.parameters['curv_entrance']))
nvaluesstr = ''
if self.parameters['nvalues'] is not None:
nvaluesstr = ' '.join([str(n) for n in self.parameters['nvalues']])
if self.parameters['ftype']=='magnetic':
commands.append('M S {:g} {:g} {:g} {:s}'
.format(self.parameters['rho0'], self.parameters['phi0']/self.parts,
self.parameters['G0'], nvaluesstr) )
else:
commands.append('E S {:g} {:g} {:g} {:s}'
.format(self.parameters['rho0'],self.parameters['phi0']/self.parts,
self.parameters['G0'], nvaluesstr) )
commands.append('FRINGING FIELD {} {:g} {:g}'.format(self.parameters['ff_exit'], self.parameters['eps_exit'], self.parameters['curv_exit']))
# Need to add drift length in order to generate the fringe field.
commands.append('DRIFT LENGTH 0.1')
refple = self.parameters['refple']
matrices = self.extractMatrices(refple.energy, refple.mass, refple.charge,
commands, self.parameters['order'] )
# Storing the transfer matrices
self.tm_ffentrance, self.tm_part, self.tm_ffexit = (matrices[1], matrices[2], matrices[3])
self.tm_ffentrance.name = 'TransferMatrix for fringe field entrance of dipole'
self.tm_ffexit.name = 'TransferMatrix for fringe field exit of dipole'
self.tm_part.name = 'TransferMatrix for dipole'
self.tm_full = None;
self.tm_full = self.getTransferMatrix()
self.parameterModified = False
self.part_ds = rho*np.deg2rad(self.parameters['phi0'])/self.parts
#### Calculating the change in profile planes for one part
phirad = np.deg2rad(self.parameters['phi0']/self.parts)
phirot = Rotation.from_rotvec([0.0,-phirad,0.0])
# Vector from entrance to exit
self.part_dr = rho*np.array( [-1+np.cos(phirad), 0,
np.sin(phirad) ] )
x0 = np.array([1.0,0.0,0.0])
y0 = np.array([0.0,1.0,0.0])
x1 = phirot.apply( x0 )
y1 = phirot.apply( y0 )
self.part_dx = x1-x0
self.part_dy = y1-y1
[docs] def length(self):
return self.totalDeltaS
[docs]class GICOSYDrift(GICOSYElement):
"""Drift element
Transports the beam through a length without electromagnetic interactions.
Parameters
----------
refple : ReferenceParticle
Reference particle for the element creation.
dl : float
Drift length of the element in m.
order : int
Order of the transfermatrix.
maxSAdvance : float
Maximal advance along the s axis in m.
collimator : float or tuple, optional
Collimator parameters if one is used.
name : string, optional
Name given to the element.
description : string, optional
Description of the element, currently not used.
Attributes
----------
collimator : Collimator or None
Collimator used in the element.
parts : int
Number of times the transfer matrix is to be used.
Notes
-----
The reference particle, the order and the driftlength are saved in the parameters.
"""
def __init__(self, refple, dl, order=1, maxSAdvance=0.2, collimator=None, name='Drift length', description='' ):
super().__init__(refple, name=name, description=description)
if collimator is not None and not isinstance(collimator, Collimator):
self.collimator = Collimator(collimator)
else:
self.collimator = collimator
self.parameters['order'] = order
self.parameters['dl'] = dl
self.parameters['ff_entrance'] = False
self.parameters['ff_exit'] = False
self.parameters['maxSAdvance'] = maxSAdvance
self.parameters['displayfield'] = False
self.parameterModified = True
self.rebuild()
[docs] def length(self):
return self.parameters['dl']
[docs] def height(self):
try:
return self.collimator.height()
except:
return 0.0
[docs] def plot(self, ax, view, s):
if self.collimator is None:
return ([], [], [])
length = self.length()
try:
if 'x' in view:
height = self.collimator.height('x')
elif 'y' in view:
height = self.collimator.height('y')
else:
raise ValueError
except:
height = self.collimator.height()
try:
thickness = self.parameters['thickness']
except:
thickness = height * 0.1
fill = 'palegoldenrod'
hatch_fill = 'royalblue'
color = 'slategray'
if PLOTBACKEND == 'matplotlib':
from matplotlib import patches
patch = patches.Rectangle((s, -height/2-thickness/2), length, height+thickness, fill=True, fc=fill, ec=color, zorder=0)
hatched_1 = patches.Rectangle((s, height/2), length, thickness, fill=True, hatch='////', fc=hatch_fill, ec=color, zorder=0)
hatched_2 = patches.Rectangle((s, -height/2-thickness), length, thickness, fill=True, hatch='////', fc=hatch_fill, ec=color, zorder=0)
draw_patches = [patch, hatched_1, hatched_2]
[docs] def rebuild(self):
if self.parameterModified == False:
return
self.parts = int(self.parameters['dl'] / self.parameters['maxSAdvance'] + 1)
commands = ['DRIFT LENGTH {:g}'.format(self.parameters['dl']/self.parts)]
refple = self.parameters['refple']
matrices = self.extractMatrices(refple.energy, refple.mass, refple.charge, commands, self.parameters['order'])
self.tm_part = matrices[0]
self.tm_full = None;
self.tm_full = self.getTransferMatrix()
self.part_ds = self.parameters['dl']/self.parts
self.part_dr = np.array([0.0,0.0,self.part_ds])
[docs]class GICOSYQuadrupole(GICOSYElement):
"""Quadrupole element
Interact with a magnetic or electric quadrupole with the given field over the given length.
Parameters
----------
refple : ReferenceParticle
Reference particle for the element creation.
ftype : {'electric', 'magnetic'}
String representing the type of quadrupole.
dl : float
Drift length of the element in m.
field : float
Voltage (V) or magnetic (T) field to be used in the quadrupole.
G0 : float
Aperture radius of the quad in m.
ff_entrance : int, optional
Fringe fields table number to use for entrance. Use zero for no fringe field.
ff_exit : int, optional
Fringe fields table number to use for exit. Use zero for no fringe field.
order : int
Order of the transfermatrix.
maxSAdvance : float
Maximal advance along the s axis in m.
collimator : float or tuple, optional
Collimator parameters if one is used.
name : string, optional
Name given to the element.
description : string, optional
Description of the element, currently not used.
Attributes
----------
collimator : Collimator or None
Collimator used in the element.
parts : int
Number of times the transfer matrix is to be used.
Notes
-----
The reference particle, the order and the driftlength are saved in the parameters.
A warning will be sent into the logger if the aperture is larger that it should be.
"""
def __init__(self, refple, ftype, dl, field, G0, ff_entrance=1, ff_exit=1, order=1, maxSAdvance=0.2, collimator=None, name='Quadrupole', description=''):
super().__init__(refple, name=name, description=description)
self.collimator = collimator if collimator is None else Collimator(collimator)
self.parameters['ftype'] = ftype
self.parameters['dl'] = dl
self.parameters['G0'] = G0
self.parameters['field'] = field
self.parameters['refple'] = refple
self.parameters['order'] = order
self.parameters['ff_entrance'] = ff_entrance
self.parameters['ff_exit'] = ff_exit
self.parameters['maxSAdvance'] = maxSAdvance
self.parameters['displayfield'] = True
self.parameterModified = True
self.rebuild()
[docs] def field(self, kind):
if kind == self.parameters['ftype']:
return self.parameters['field']
else:
return 0.0
[docs] def length(self):
return self.parameters['dl']
[docs] def height(self):
return self.parameters['G0'] * 2
[docs] def plot(self, ax, view, s):
length = self.length()
height = self.height()
try:
thickness = self.parameters['thickness']
except:
thickness = self.parameters['G0']*1.148
fill = 'palegoldenrod'
hatch_fill = 'royalblue'
color = 'slategray'
if PLOTBACKEND == 'matplotlib':
from matplotlib import patches
patch = patches.Rectangle((s, -height/2-thickness/2), length, height+thickness, fill=True, fc=fill, ec=color, zorder=0)
if self.parameters['ftype'] == 'electric':
hatched_1 = patches.Rectangle((s, height/2), length, thickness, fill=True, hatch='////', fc=hatch_fill, ec=color, zorder=0)
hatched_2 = patches.Rectangle((s, -height/2-thickness), length, thickness, fill=True, hatch='////', fc=hatch_fill, ec=color, zorder=0)
# draw_patches.extend([hatched_1, hatched_2])
ax.add_patch(patch)
if self.parameters['ftype'] == 'electric':
ax.add_patch(hatched_1)
ax.add_patch(hatched_2)
[docs] def rebuild(self):
if self.parameterModified == False:
return
#print( "rebuilding drift with {:f}".format( self.parameters['dl']) )
self.parts = int(self.parameters['dl'] / self.parameters['maxSAdvance'] + 1)
if not self.collimator:
self.collimator = Collimator(self.parameters['G0'])
else:
try:
if self.collimator.parameters['size'][0] > self.parameters['G0']:
logger.warning('Quadrupole aperture is larger in the x direction than the pole gap! Results can be unphysical!')
elif self.collimator.parameters['size'][1] > self.parameters['G0']:
logger.warning('Quadrupole aperture is larger in the y direction than the pole gap! Results can be unphysical!')
else:
pass
except:
if self.collimator.parameters['size'] > self.parameters['G0']:
logger.warning('Quadrupole aperture is larger than the pole gap! Results can be unphysical!')
else:
pass
commands = []
# Need to have a drift length in order to generate the
# fringe field. This matrix is not used nor passed for a user.
commands.append('DRIFT LENGTH 0.1')
commands.append('FRINGING FIELD {}'.format(self.parameters['ff_entrance']))
s = self.parameters['dl']/self.parts
F = self.parameters['field']
G0 = self.parameters['G0']
refple = self.parameters['refple']
if self.parameters['ftype'] == 'magnetic':
commands.append('M Q {:g} {:g} {:g}'.format(s, F, G0))
else:
commands.append('E Q {:g} {:g} {:g}'.format(s, 1e-3*F, G0))
commands.append('FRINGING FIELD {}'.format(self.parameters['ff_exit']))
# Need to add drift length in order to generate the fringe field.
commands.append('DRIFT LENGTH 0.1')
matrices = self.extractMatrices(refple.energy, refple.mass, refple.charge, commands, self.parameters['order'])
self.tm_ffentrance = matrices[1]
self.tm_part = matrices[2]
self.tm_ffexit = matrices[3]
self.tm_full = None
self.tm_full = self.getTransferMatrix()
self.part_ds = self.parameters['dl']/self.parts
self.part_dr = np.array([0.0,0.0,self.part_ds])
ff_format_str = """*-----------------------------------
{} | {} | | {} {} {} {}
| | | {} {} {} {}
| | | {} {} {}
"""
MSparams = {
1: (0.205133, 0.840972, -0.141308, 0.050050, 0.000076, 0.005197, 0, 0, 0, 0, 0),
2: (2, 0.20513185638, 0.8409723305, -0.141308269, 5.00500213E-02, 7.5965396644E-05, 5.1969989319E-03, 0, 0, 0, 0, 0),
3: (0.38855444, 2.40798992, -0.834568931, 0.19822262, 2.444162E-03, -3.187845E-04, 0, 0, 0, 0, 0),
4: (0.061479475, 1.771726174, -0.158179632, 0.056507703, 0.033236894, 0.005135692, 0, 0, 0, 0, 0),
5: (0.056658, 1.008271, -0.016882, -0.003549, 0.000139, 0.000013, 0, 0, 0, 0, 0),
6: (0.478959, 0.955644, -0.296488, -0.067666, 0.009941, 0.176659, 0, 0, 0, 0, 0),
7: (0.183486031, 2.573644333, -0.207994544, -0.091821834, -0.002378879, 0.003844973, 0, 0, 0, 0, 0),
8: (0.285850000, 1.235990000, -0.291610000, 0.443700000, -0.205160000, 0.055770000, 0, 0, 0, 0, 0),
}
ESparams = {
1: (0.2049, 0.84105, -0.14135, 0.05005, 6.875E-5, 5.1968E-3, 0, 0, 0, 0, 0),
2: (0.75, 0.09, 90, 0, 0, 0, 0, 0, 0, 0, 0),
3: (0.38855444, 2.40798992, -0.834568931, 0.19822262, 2.444162E-03, -3.187845E-04, 0, 0, 0, 0, 0),
4: (0.061479475, 1.771726174, -0.158179632, 0.056507703, 0.033236894, 0.005135692, 0, 0, 0, 0, 0),
7: (0.183486031, 2.573644333, -0.207994544, -0.091821834, -0.002378879, 0.003844973, 0, 0, 0, 0, 0),
}
MQparams = {
1: (0, 3.59463, 0, -0, 0, 0, 0, 0, 0, 0, 0),
2: (0, 3.59463, 0, -0, 0, 0, 0, 0, 0, 0, 0),
3: (0.10391, 6.27098, -1.51225, 3.59925, -2.13239, 1.72300, 0, 0, 0, 0, 0),
4: (0.388554, 2.407990, -0.834569, 0.198223, 0.002444, -0.000319, 0, 0, 0, 0, 0),
5: (0.3091603650, 2.3792544734, -0.4781358927, 3.56933296E-02, 3.980007740E-04, 6.649291231E-04, 0, 0, 0, 0, 0),
6: (0.2964295664, 2.266630258, -0.56775282622, 0.1335785413, -2.274501006530167E-03, 6.95656239929E-04, 0, 0, 0, 0, 0),
9: (0.271647162, 2.83514462, -0.602080979, 2.604809986E-02, 7.275585244E-03, 1.550781439E-03, 0, 0, 0, 0, 0),
}
EQparams = {
1: (0, 3.59276, 0, -0.081520, 0, -0.00177, 0, 0, 0, 0, 0),
2: (-0.77363101, 8.26159995, 0, 0, 0, 0, 0, 0, 0, 0, 0),
3: (0.10391, 6.27098, -1.51225, 3.59925, -2.13239, 1.72300, 0, 0, 0, 0, 0),
4: (0.388554, 2.407990, -0.834569, 0.198223, 0.002444, -0.000319, 0, 0, 0, 0, 0),
5: (0, 2.407990, 0, 0, 0, 0, 0, 0, 0, 0, 0),
6: (0.2964295664, 2.266630258, -0.56775282622, 0.1335785413, -2.274501006530167E-03, 6.95656239929E-04, 0, 0, 0, 0, 0),
7: (0.176659, 3.576539, -0.778279, 0.430539, -0.123546, 0.016877, 0, 0, 0, 0, 0),
9: (0.271647162, 2.83514462, -0.602080979, 2.604809986E-02, 7.275585244E-03, 1.550781439E-03, 0, 0, 0, 0, 0),
}
[docs]def setFringeField(element, tableno, params, force=False):
"""Set fringe field for an element type and table number tableno.
By default PIOL uses a built-in standard table for fringe fields. Currently the fringe field
data contains mostly nonsense.
Parameters
----------
element : str
The element type. One of {'MS', 'ES', 'MQ', 'EQ'}.
tableno : int
The table number to set. Valid table numbers are from 1 to 9.
params : array-like
An array of polynomial coefficients describing the fringe field with the Enge model.
Lowest order first. Maximum of 11th order polynomial is supported.
force : boolean
Forces the overwriting of the GICOSYFF.DAT file. Defaults to 'False'.
"""
params = np.array(params)
if len(params) > 11:
raise ValueError('Enge model supports a maximum of 11th order polynomial, {:d} terms was given!'.format(len(params)))
if params[params.nonzero()][-1] < 0:
raise ValueError('Highest order term is zero in Enge model - not physical!')
mapping = {'MS': MSparams, 'ES': ESparams, 'MQ': MQparams, 'EQ': EQparams}
if element not in mapping.keys():
raise ValueError('Unknown element type {} for fringe field!'.format(element))
tableno = int(tableno)
if (tableno > 9 ) or (tableno < 1):
raise ValueError('Invalid table number {}, use numbers 1-9!'.format(tableno))
params.resize(11)
mapping[element][tableno] = params
writeFringeFields(force=force)
[docs]def writeFringeFields(force=False):
"""Writes the GICOSYFF.DAT fringe field data file.
Overwrites the fringe field data file GICOSYFF.DAT if marked with tag '* PIOL', otherwise
rises an exception. Overwriting the fringe field data file can be forced by setting
force=True.
"""
readGICOSYPath()
global _gicosy_path
mapping = {'MS': MSparams, 'ES': ESparams, 'MQ': MQparams, 'EQ': EQparams}
# Read the GICOSYFF.DAT and check the tag
if not force:
with cd(_gicosy_path):
with open('GICOSYFF.DAT', 'r') as gff:
s = gff.readline()
if s[:6] != '* PIOL':
raise Exception('GICOSYFF.DAT not written by PIOL. Avoiding overwriting.\nBackup if necessary and use writeFringeFields(force=True) to proceed.')
# Writing the GICOSYFF.DAT fringe field data file
s = '* PIOL has written this file\n'
for key in mapping.keys():
if key.lower() == 'custom':
continue
for k in mapping[key].keys():
s += ff_format_str.format(key, k, *mapping[key][k])
with cd(_gicosy_path):
with open('GICOSYFF.DAT', 'w') as gff:
gff.write(s)
def Enge(z, a, D):
""""The Enge fringe field shape evaluator.
Returns the relative field strength as a function of position along the beam line according to the
Enge model, i.e.
E = 1 / (1+exp( a[0] + a[1]*z/D + a[2]*(z/D)^2 + a[3]*(z/D)^3 + ... ))
Parameters:
-----------
z : float
Coordinate along the beam line relative to the physical edge of the element. The
Enge fringe field model describes the field going into the element, i.e. negative
coordinates are outside the element (where this function returns 0) and positive
coordinates are inside (where this function returns 1).
a : ndarray
An array of polynomial coefficients describing the fringe field with the Enge model.
Lowest order first.
D : The physical scaling factor for the element (usually half the element gap size).
"""
z = -z/D
v = np.polyval(a[::-1], z)
E = 1/(1+np.exp(v))
return E
def calculateEnge(element, tableno, face='entrance'):
"""Calculate the Enge fringe field for element
Returns the field shape as a function of z-coordinate near the element edge. The returned
(z,E) contains value pairs, which describe the field shape from field value of 0.0001 to
0.9999.
Parameters:
-----------
element : Element
The element object
tableno : int
The table number to use. Valid table numbers are from 1 to 9.
face : str
Either 'entrance' or 'exit'
Returns:
--------
tuple
(z,E), where z is the z-coordinates and E is the scaled fringe field.
"""
D = element.height()/2
if type(element) == GICOSYDipole:
map = {'electric': ESparams, 'magnetic': MSparams}
elif type(element) == GICOSYQuadrupole:
map = {'electric': EQparams, 'magnetic': MQparams}
else:
raise Exception('Unknown element type')
params = map[element.parameters['ftype']][tableno]
try:
from scipy.optimize import root_scalar
function = lambda z: Enge(z, params, D) - 0.9999
result_max = root_scalar(function, x0=0, method='brentq', bracket=[0, 10*D], maxiter=100)
function = lambda z: Enge(z, params, D) - 0.0001
result_min = root_scalar(function, x0=0, method='brentq', bracket=[-10*D, 0], maxiter=100)
z = np.linspace(result_min.root, result_max.root, 100)
except:
logger.warning('Could not find fringe field limits')
z = np.linspace(-D, D, 100)
y = Enge(z, params, D)
if type(element) == GICOSYDipole:
# Adjust scale along optical axis due to dipole face angle
if face == 'entrance':
z = z/np.cos(np.radians(element.parameters['eps_entrance']))
else:
z = z/np.cos(np.radians(element.parameters['eps_exit']))
return z, y
