diff --git a/ess/lib_tw.py b/ess/lib_tw.py
index dd0a8800da79bbbf12df59e410940a225cd60800..7e760beeda3e3f37afb47a36810397c08bad54fd 100644
--- a/ess/lib_tw.py
+++ b/ess/lib_tw.py
@@ -19,6 +19,7 @@ from itertools import chain
from subprocess import check_output
from os import system
from os.path import isdir,isfile
+from sys import exit
from lib_tw_elem import *
@@ -29,14 +30,14 @@ from lib_tw_elem import *
#---- Lattice and project
class LATTICE:
- '''
+ '''
2015.10.14
'''
def __init__(self,file_name_lat,file_name_fmap,freq=352.21,gamma=1.0):
# In case file_name_fmap is str
if isinstance(file_name_fmap,basestring)==True: file_name_fmap=[file_name_fmap]
-
+
# Elem/comm class dict
dic_cls={'DRIFT':DRIFT,
'QUAD':QUAD,'BEND':BEND,'EDGE':EDGE,'THIN_STEERING':THIN_STEERING,
@@ -44,7 +45,7 @@ class LATTICE:
'APERTURE':APERTURE,'DIAG_POSITION':DIAG_POSITION,
'STEERER':STEERER,'CHOPPER':CHOPPER,
'FREQ':FREQ,'MARKER':MARKER}
-
+
# Field map dict
dic_fmap={}
for file_name_fmap_i in file_name_fmap:
@@ -109,8 +110,8 @@ class LATTICE:
for k in range(i)[::-1]:
try : self.lst[i].apt=self.lst[k].apt; break
except: pass
- except: pass
-
+ except: pass
+
def update_gamma(self):
for i in range(len(self.lst)):
@@ -132,7 +133,7 @@ class LATTICE:
if lin[i]==':' : idx_elem =int(lin[i-1])-1 # "-1" since idx=0,1,2,...
if lin[i]=='BLx=': BLx[idx_elem]=float(lin[i+1])
if lin[i]=='BLy=': BLy[idx_elem]=float(lin[i+1])
-
+
# Assign/update steerers
for i in range(len(self.lst)):
if self.lst[i].idx_elem in BLx.keys():
@@ -146,7 +147,7 @@ class LATTICE:
self.lst[k].Bx=BLx[idx_elem]/self.lst[i].L
self.lst[k].By=BLy[idx_elem]/self.lst[i].L
break
-
+
def get_tw(self,file_name):
with open(file_name,'w') as file:
@@ -156,7 +157,7 @@ class LATTICE:
def get_madx(self,file_name_elem='elem.madx',file_name_seq='seq.madx'):
if self.lst[-1].gamma==1.0: self.update_gamma() # Assign gamma, if not done yet
-
+
with open(file_name_elem,'w') as file:
for lat_i in self.lst:
try : print >>file,lat_i.get_madx()
@@ -173,9 +174,9 @@ class LATTICE:
for lat_i in self.lst:
try : print >>file,lat_i.get_fluka()
except: pass
-
+
def get_mars(self,file_name='elem.dat'):
-
+
if self.lst[-1].gamma==1.0: self.update_gamma() # Assign gamma, if not done yet
with open(file_name,'w') as file:
@@ -187,7 +188,7 @@ class PROJECT:
'''
- This is for running multiple simulations 1-by-1 under 1 project.
- Maybe not very useful...
-
+
2015.10.15
'''
@@ -210,7 +211,7 @@ class PROJECT:
## if self.path_cal!=None:
## if isdir(self.path_cal)==False: system('mkdir '+self.path_cal)
-
+
system(opt_exe)
#---- Data related
@@ -218,14 +219,14 @@ class PROJECT:
class PARTRAN:
'''
- The list not complete. Add parameters as needed.
-
- 2015.11.24
+
+ 2016.01.18
'''
def __init__(self,file_name):
# Consts to convert phs to z.
mass=938.272029
- c =2.99792458
+ c =2.99792458
freq=352.21
# Extract indices.
@@ -233,15 +234,15 @@ class PARTRAN:
for lin in file.readlines():
lin=lin.split()
if lin[0]=='####':
- idx_s =lin.index("z(m)" ); idx_gamma=lin.index("gama-1")
- idx_x =lin.index("x0" ); idx_y =lin.index("y0" ); idx_phs =lin.index("p0" )
- idx_sigx=lin.index("SizeX" ); idx_sigy =lin.index("SizeY" ); idx_sigz=lin.index("SizeZ"); idx_sigp=lin.index("SizeP")
- idx_alfx=lin.index("sxx'" ); idx_alfy =lin.index("syy'" ); idx_alfz=lin.index("szdp" )
- idx_epsx=lin.index("ex" ); idx_epsy =lin.index("ey" ); idx_epsz=lin.index("ezdp" ); idx_epsp=lin.index("ep" )
- idx_halx=lin.index("hx" ); idx_haly =lin.index("hy" ); idx_halp=lin.index("hp" )
- idx_loss=lin.index("Powlost")
+ idx_s =lin.index("z(m)" ); idx_gamma=lin.index("gama-1" )
+ idx_x =lin.index("x0" ); idx_y =lin.index("y0" ); idx_phs =lin.index("p0" )
+ idx_sigx =lin.index("SizeX" ); idx_sigy =lin.index("SizeY" ); idx_sigz=lin.index("SizeZ"); idx_sigp=lin.index("SizeP")
+ idx_alfx =lin.index("sxx'" ); idx_alfy =lin.index("syy'" ); idx_alfz=lin.index("szdp" )
+ idx_epsx =lin.index("ex" ); idx_epsy =lin.index("ey" ); idx_epsz=lin.index("ezdp" ); idx_epsp=lin.index("ep" )
+ idx_halx =lin.index("hx" ); idx_haly =lin.index("hy" ); idx_halp=lin.index("hp" )
+ idx_Nptcl=lin.index("npart" ); idx_loss =lin.index("Powlost")
break
-
+
# Extract data.
with open(file_name) as file:
data=[]; flag=0
@@ -252,12 +253,12 @@ class PARTRAN:
data=array(data).transpose()
# Instances
- self.s =data[idx_s ]
- self.x =data[idx_x ]; self.y =data[idx_y ]; self.phs =data[idx_phs ]
- self.sigx=data[idx_sigx]; self.sigy=data[idx_sigy]; self.sigz=data[idx_sigz]; self.sigp=data[idx_sigp]
- self.epsx=data[idx_epsx]; self.epsy=data[idx_epsy]; self.epsz=data[idx_epsz]; self.epsp=data[idx_epsp]
- self.halx=data[idx_halx]; self.haly=data[idx_haly]; self.halz=data[idx_halp]; self.halp=data[idx_halp]
- self.loss=data[idx_loss]
+ self.s =data[idx_s ]
+ self.x =data[idx_x ]; self.y =data[idx_y ]; self.phs =data[idx_phs ]
+ self.sigx =data[idx_sigx ]; self.sigy=data[idx_sigy]; self.sigz=data[idx_sigz]; self.sigp=data[idx_sigp]
+ self.epsx =data[idx_epsx ]; self.epsy=data[idx_epsy]; self.epsz=data[idx_epsz]; self.epsp=data[idx_epsp]
+ self.halx =data[idx_halx ]; self.haly=data[idx_haly]; self.halz=data[idx_halp]; self.halp=data[idx_halp]
+ self.Nptcl=data[idx_Nptcl]; self.loss=data[idx_loss]
# Additional instances
self.gamma= data[idx_gamma]+1.0
@@ -271,28 +272,28 @@ class PARTRAN:
self.alfy =-data[idx_alfy]/self.epsy*self.beta*self.gamma
self.alfz =-data[idx_alfz]/self.epsz*self.beta*self.gamma**3
self.alfp =-self.alfz
-
+
# Convert to list (not necessary?)
- self.s =self.s.tolist() ; self.gamma=self.gamma.tolist(); self.beta=self.beta.tolist()
- self.x =self.x.tolist() ; self.y =self.y.tolist() ; self.z =self.z.tolist() ; self.phs =self.phs.tolist()
- self.sigx=self.sigx.tolist(); self.sigy =self.sigy.tolist() ; self.sigz=self.sigz.tolist(); self.sigp=self.sigp.tolist()
- self.betx=self.betx.tolist(); self.bety =self.bety.tolist() ; self.betz=self.betz.tolist(); self.betp=self.betp.tolist()
- self.alfx=self.alfx.tolist(); self.alfy =self.alfy.tolist() ; self.alfz=self.alfz.tolist(); self.alfp=self.alfp.tolist()
- self.epsx=self.epsx.tolist(); self.epsy =self.epsy.tolist() ; self.epsz=self.epsz.tolist(); self.epsp=self.epsp.tolist()
- self.halx=self.halx.tolist(); self.haly =self.haly.tolist() ; self.halz=self.halz.tolist(); self.halp=self.halp.tolist()
- self.loss=self.loss.tolist()
-
- def loss_den(self,file_name_dt,dlt_dt=5e-6):
+ self.s =self.s.tolist() ; self.gamma=self.gamma.tolist(); self.beta=self.beta.tolist()
+ self.x =self.x.tolist() ; self.y =self.y.tolist() ; self.z =self.z.tolist() ; self.phs =self.phs.tolist()
+ self.sigx =self.sigx.tolist() ; self.sigy =self.sigy.tolist() ; self.sigz=self.sigz.tolist(); self.sigp=self.sigp.tolist()
+ self.betx =self.betx.tolist() ; self.bety =self.bety.tolist() ; self.betz=self.betz.tolist(); self.betp=self.betp.tolist()
+ self.alfx =self.alfx.tolist() ; self.alfy =self.alfy.tolist() ; self.alfz=self.alfz.tolist(); self.alfp=self.alfp.tolist()
+ self.epsx =self.epsx.tolist() ; self.epsy =self.epsy.tolist() ; self.epsz=self.epsz.tolist(); self.epsp=self.epsp.tolist()
+ self.halx =self.halx.tolist() ; self.haly =self.haly.tolist() ; self.halz=self.halz.tolist(); self.halp=self.halp.tolist()
+ self.Nptcl=self.Nptcl.tolist(); self.loss =self.loss.tolist()
+
+ def loss_den(self,file_name_dt='',dlt_dt=5e-6):
return loss_elem2den(self.s,self.loss,file_name_dt,dlt_dt)
class DST:
'''
Class for a TraceWin's .dst file.
-
+
- TraceWin seems using beta and gamma for each particle
so the conversion to (z,z') is based on this assumption.
-
+
2015.10.06
'''
def __init__(self,file_name,unit_x='cm',unit_px='rad',unit_z='rad',unit_pz='MeV'):
@@ -309,7 +310,7 @@ class DST:
dummy=fromfile(file,dtype= uint8 ,count=1 )
x =fromfile(file,dtype=float64,count=Nptcl*6).reshape(Nptcl,6).transpose()
mass =fromfile(file,dtype=float64,count=1 )[0]
-
+
# Adjust units
gamma=1.0+x[5]/mass; beta=sqrt(1-1/gamma**2)
if unit_x =='mm' : x[0]= x[0]*1e1; x[2]=x[2]*1e1
@@ -319,7 +320,7 @@ class DST:
if unit_pz=='mrad': x[5]=(x[5]-mean(x[5]))/(mass*beta**2*gamma**3)*1e3
# Instances
- self.x =x.transpose()
+ self.x =x.transpose()
self.mass =mass
self.freq =freq
self.Ibeam=Ibeam
@@ -341,15 +342,15 @@ class DENSITY:
Nrun x Nelem: loss
Nelem : loss_num_(..), loss_pow_(..)
-
- Nelem x Nstep: den
-
+
+ Nelem x Nstep: den
+
2015.11.12
'''
def __init__(self,file_name):
#-- Empty arrays
-
+
idx_elem =[]; s =[]
apt =[]; accpt =[]
Nptcl =[]; Ibeam =[]
@@ -364,9 +365,9 @@ class DENSITY:
loss_num_max=[]; loss_pow_max=[]
loss_num_min=[]; loss_pow_min=[]
den =[]; den_pow =[]
-
+
#-- Extract data
-
+
with open(file_name) as file:
while True:
try:
@@ -404,7 +405,7 @@ class DENSITY:
loss_num_rms.append(fromfile(file,dtype= uint64,count=1)[0])
loss_num_min.append(fromfile(file,dtype= uint64,count=1)[0])
loss_num_max.append(fromfile(file,dtype= uint64,count=1)[0])
- loss_pow_ave.append(sum(loss_pow[-1]))
+ loss_pow_ave.append(sum(loss_pow[-1]))
loss_pow_rms.append(fromfile(file,dtype=float64,count=1)[0])
loss_pow_min.append(fromfile(file,dtype=float32,count=1)[0])
loss_pow_max.append(fromfile(file,dtype=float32,count=1)[0])
@@ -429,18 +430,18 @@ class DENSITY:
if Nptcl[0]>0:
loss_num=swapaxes(loss_num,1,0); loss_pow=swapaxes(loss_pow,1,0)
den =swapaxes( den,1,0); den_pow =swapaxes( den_pow,1,0)
-
+
#-- Take care ave and rms
cent_ave=cent_ave/Nrun; cent_rms=sqrt(cent_rms/Nrun)
- sig_ave = sig_ave/Nrun; sig_rms =sqrt( sig_rms/Nrun)
+ sig_ave = sig_ave/Nrun; sig_rms =sqrt( sig_rms/Nrun)
eps_ave = eps_ave/Nrun; eps_rms =sqrt( eps_rms/Nrun)
if Nptcl[0]>0:
loss_num_ave=1.0*array(loss_num_ave)/Nrun; loss_num_rms=sqrt(1.0*array(loss_num_rms)/Nrun)
loss_pow_ave= array(loss_pow_ave)/Nrun; loss_pow_rms=sqrt( array(loss_pow_rms)/Nrun)
#-- Change units, m => mm, pi-m-rad => pi-mm-mrad
-
+
apt *=1e3
eps_ave*=1e6; eps_rms*=1e6
for k in (0,1,4,5):
@@ -450,7 +451,7 @@ class DENSITY:
xmax[k] *=1e3; xmin[k] *=1e3
#-- Define std (around to avoid sqrt(-eps))
-
+
if Nrun>1:
cent_std=sqrt(around(cent_rms**2-cent_ave**2,12))
sig_std =sqrt(around( sig_rms**2- sig_ave**2,12))
@@ -463,7 +464,7 @@ class DENSITY:
loss_pow_std=sqrt(around(loss_pow_rms**2-loss_pow_ave**2,16))
loss_num_std=nan_to_num(loss_num_std) # Replace nan with 0
loss_pow_std=nan_to_num(loss_pow_std) # Replace nan with 0
-
+
#-- Convert to list (just in case...)
apt = apt.tolist(); accpt = accpt.tolist(); accpt.append(accpt[0]); del accpt[0]
@@ -484,20 +485,20 @@ class DENSITY:
eps_std = eps_std.tolist()
if Nptcl[0]>0:
loss_num_std=loss_num_std.tolist()
- loss_pow_std=loss_pow_std.tolist()
+ loss_pow_std=loss_pow_std.tolist()
#-- Outputs
-
+
self.idx_elem=idx_elem
self.s =s
self.apt =apt
self.accpt =accpt
self.Nptcl =Nptcl
self.Ibeam =Ibeam
-
+
self.Nrun =Nrun
self.Nstep=Nstep
-
+
if Nrun==1:
self.cent=cent_ave
self.sig = sig_ave
@@ -517,7 +518,7 @@ class DENSITY:
self.loss_num =loss_num ; self.loss_pow =loss_pow
self.loss_num_ave=loss_num_ave; self.loss_pow_ave=loss_pow_ave
self.loss_num_rms=loss_num_rms; self.loss_pow_rms=loss_pow_rms
- self.loss_num_std=loss_num_std; self.loss_pow_std=loss_pow_std
+ self.loss_num_std=loss_num_std; self.loss_pow_std=loss_pow_std
self.loss_num_max=loss_num_max; self.loss_pow_max=loss_pow_max
self.loss_num_min=loss_num_min; self.loss_pow_min=loss_pow_min
self.den = den ; self.den_pow = den_pow
@@ -533,14 +534,14 @@ class DENSITY:
def x2dst(x,mass,freq,Ibeam,path_file='part_dtl1_new.dst'):
'''
Output a TraceWin's .dst file from x and etc.
-
+
Input: x (Nptcl,6)
For binary characters see https://docs.python.org/2/library/struct.html
2014.10.03
'''
-
+
file=open(path_file,'w')
out =pack('b',125)
out+=pack('b',100)
@@ -551,23 +552,23 @@ def x2dst(x,mass,freq,Ibeam,path_file='part_dtl1_new.dst'):
x=list(chain(*x)) # Flatten x
for x_i in x: out+=pack('d',x_i)
out+=pack('d',mass)
- print >>file,out
+ print >>file,out
file.close()
-
+
def plt2x(path_file):
'''
Extract x and etc from a TraceWin's binary .plt file.
The units are (cm,rad,cm,rad,rad,MeV,loss).
-
+
Output: x (Nelem,Nptcl,7)
For binary characters see https://docs.python.org/2/library/struct.html
2014.10.06
'''
-
+
file=open(path_file,'r')
-
+
# Hetero info
file.read(1); file.read(1) # 125,100
Nelem=unpack('i',file.read(4))[0]
@@ -575,9 +576,9 @@ def plt2x(path_file):
Ibeam=unpack('d',file.read(8))[0]
freq =unpack('d',file.read(8))[0]
mass =unpack('d',file.read(8))[0]
-
+
# Loop for elem
- x=[]; i_unk=[]; i_elem=[]; s=[]; phs=[]; ken=[];
+ x=[]; i_unk=[]; i_elem=[]; s=[]; phs=[]; ken=[];
for i in range(Nelem):
try:
i_unk.append(unpack('b',file.read(1))[0])
@@ -589,7 +590,7 @@ def plt2x(path_file):
except:
break
- file.close()
+ file.close()
return x,mass,freq,Ibeam,i_unk,i_elem,s,phs,ken
@@ -600,10 +601,10 @@ def x2plt(x,mass,freq,Ibeam,i_unk,i_elem,s,phs,ken,path_file='dtl1_new.plt'):
Input: x (Nelem,Nptcl,7)
For binary characters see https://docs.python.org/2/library/struct.html
-
+
2014.10.07
'''
-
+
file=open(path_file,'w')
out =pack('b',125)
out+=pack('b',100)
@@ -625,17 +626,17 @@ def x2plt(x,mass,freq,Ibeam,i_unk,i_elem,s,phs,ken,path_file='dtl1_new.plt'):
file.close()
#---- Data related
-
+
def x2twiss(x):
'''
Calculate twiss from x. Note sig = sqrt(bet*eps).
-
+
Input : x (Nptcl,6)
Output: eps,bet,alf,gam
2015.07.30
'''
-
+
x =[x_i-mean(x_i) for x_i in array(x).transpose()]
sig=[[mean(x_i*x_k) for x_k in x] for x_i in x]
eps=[ sqrt(det(array(sig)[i:i+2][:,i:i+2])) for i in (0,2,4)]
@@ -649,13 +650,13 @@ def x2twiss_halo(x,sig_cut=(None,None,None)):
'''
Calculate twiss of halo particles (r > sig_cut).
Note halo particles are different for each plane.
-
+
Input : x (Nptcl,6), sig_cut
Output: eps,bet,alf,gam
2015.07.30
'''
-
+
eps,bet,alf,gam=x2twiss(x)
for k in (0,2,4):
if sig_cut[k/2]!=None:
@@ -670,35 +671,73 @@ def x2twiss_halo(x,sig_cut=(None,None,None)):
gam[k/2]=gam_tmp[k/2]
return eps,bet,alf,gam
-
-def loss_elem2den(s,loss,path_file_dt,dlt_dt=5e-6):
+
+def loss_elem2den(s,loss,file_name_dt='',dlt_dt=5e-6):
'''
- This function calculates loss density [W/m] from s, loss,
- and the file of the DTL's cell and DT lengths.
+ This function calculates loss density [W/m] from s and loss
+ together with the file of the DTL's cell and DT lengths.
- 2015.01.30
+ 2016.01.15
'''
- # Define l
- l=[s[i+1]-s[i] for i in range(len(s)-1)]; l.insert(0,0.0)
-
- # DTL cell and DT lengths.
- file=open(path_file_dt)
- l_cell,l_dt=array([map(float,lin.split()) for lin in file.readlines()]).transpose()[:2]
- file.close()
-
- # Replace l's in DTL with DT lengths.
- Ndt=0
- for i in range(len(l)):
- dl=list(abs(l_cell-l[i])); dl_min=min(dl)
- if dl_min<dlt_dt:
- l[i]=l_dt[dl.index(dl_min)]; Ndt+=1
- # Check the replacement.
- if Ndt!=len(l_cell):
- print 'drift-tube # unmatched, in file: '+str(len(l_cell))+', from matching: '+str(Ndt)
- print 'review the threshold, exiting...'; exit()
+ # Define L
+ L=[s[i]-s[i-1] for i in range(1,len(s))]; L.insert(0,0.0)
+
+ # Take care DTL part if "file_name_dt" is given
+ try:
+ # DTL cell and DT lengths
+ with open(file_name_dt) as file:
+ L_cell,L_dt=array([map(float,lin.split()) for lin in file.readlines()]).transpose()[:2]
+ # Replace cell lengths with DT lengths
+ Ndt=0
+ for i in range(len(L)):
+ dL=list(abs(L_cell-L[i])); dL_min=min(dL)
+ if dL_min<dlt_dt:
+ L[i]=L_dt[dL.index(dL_min)]; Ndt+=1
+ # Check
+ if Ndt!=len(L_dt):
+ print 'drift-tube # unmatched, in file: '+str(len(L_dt))+', matched: '+str(Ndt)
+ except:
+ pass
+
+ # Take care loss in 0 length elem
+ loss_den=loss[::]
+ for i in range(len(loss_den)):
+ if loss_den[i]!=0.0 and L[i]==0.0:
+ print 'Caution: inf loss density at elem # ',i,'!!!!'
+ for k in range(i)[::-1]:
+ if L[k]!=0.0: loss_den[k]+=loss_den[i]; loss_den[i]=0.0; break
+
+ # Calculate loss den
+ for i in range(len(loss_den)):
+ if loss_den[i]!=0.0: loss_den[i]/=L[i]
+
+ return loss_den
+
+ exit()
- # Treat inf loss den.
+ # Define L
+ l=[s[i+1]-s[i] for i in range(len(s)-1)]; l.insert(0,0.0)
+ #L=[s[i+1]-s[i] for i in range(len(s)-1)]; L.insert(0,0.0)
+
+ # Take care DTL part if "file_name_dt" is given
+ try:
+ # Read DTL cell and DT lengths
+ with open(file_name_dt) as file:
+ l_cell,l_dt=array([map(float,lin.split()) for lin in file.readlines()]).transpose()[:2]
+ # Replace cell lengths with DT lengths
+ Ndt=0
+ for i in range(len(l)):
+ dl=list(abs(l_cell-l[i])); dl_min=min(dl)
+ if dl_min<dlt_dt:
+ l[i]=l_dt[dl.index(dl_min)]; Ndt+=1
+ # Check
+ if Ndt!=len(l_cell):
+ print 'drift-tube # unmatched, in file: '+str(len(l_cell))+', matched: '+str(Ndt)
+ except:
+ pass
+
+ # Take care inf loss den
for i in range(len(l)):
if l[i]==0.0 and loss[i]!=0.0:
if i==1:
@@ -708,15 +747,61 @@ def loss_elem2den(s,loss,path_file_dt,dlt_dt=5e-6):
while l[i-k]==0.0: k+=1
loss[i-k]+=loss[i]; loss[i]=0.0
print 'Caution: Inf loss density at elem #',i,'!!!!'
-
- # Finally, convert to loss density.
+
+ # Finally, convert to loss density
loss_den=[]
for i in range(len(l)):
if loss[i]==0.0: loss_den.append(0.0)
else : loss_den.append(loss[i]/l[i])
-
+
return loss_den
+## def loss_elem2den(s,loss,path_file_dt,dlt_dt=5e-6):
+## '''
+## This function calculates loss density [W/m] from s, loss,
+## and the file of the DTL's cell and DT lengths.
+
+## 2015.01.30
+## '''
+## # Define l
+## l=[s[i+1]-s[i] for i in range(len(s)-1)]; l.insert(0,0.0)
+
+## # DTL cell and DT lengths.
+## file=open(path_file_dt)
+## l_cell,l_dt=array([map(float,lin.split()) for lin in file.readlines()]).transpose()[:2]
+## file.close()
+
+## # Replace l's in DTL with DT lengths.
+## Ndt=0
+## for i in range(len(l)):
+## dl=list(abs(l_cell-l[i])); dl_min=min(dl)
+## if dl_min<dlt_dt:
+## l[i]=l_dt[dl.index(dl_min)]; Ndt+=1
+
+## # Check the replacement.
+## if Ndt!=len(l_cell):
+## print 'drift-tube # unmatched, in file: '+str(len(l_cell))+', from matching: '+str(Ndt)
+## print 'review the threshold, exiting...'; exit()
+
+## # Treat inf loss den.
+## for i in range(len(l)):
+## if l[i]==0.0 and loss[i]!=0.0:
+## if i==1:
+## print 'The 1st elem has 0 length and yet has a loss, exiting...'; exit()
+## else:
+## k=1
+## while l[i-k]==0.0: k+=1
+## loss[i-k]+=loss[i]; loss[i]=0.0
+## print 'Caution: Inf loss density at elem #',i,'!!!!'
+
+## # Finally, convert to loss density.
+## loss_den=[]
+## for i in range(len(l)):
+## if loss[i]==0.0: loss_den.append(0.0)
+## else : loss_den.append(loss[i]/l[i])
+
+## return loss_den
+
def partran_end_all(file_name_in_all,file_name_out):
'''
- Reads multiple partran1.out files and extracts the data at the end.
@@ -725,7 +810,7 @@ def partran_end_all(file_name_in_all,file_name_out):
2015.11.25
'''
-
+
# Extract header from 1 file.
with open(file_name_in_all[0],'r') as file:
header=''
@@ -744,7 +829,7 @@ def partran_end_all(file_name_in_all,file_name_out):
for i in range(1,len(lin)):
if '.' in lin[i]: file_out.write(' %13s'%lin[i])
else : file_out.write(' %10s'%lin[i])
- file_out.write('\n')
+ file_out.write('\n')
#-------- Obsolete classes and functions
@@ -757,7 +842,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## 2015.04.26
## '''
## def __init__(self,typ):
-
+
## # Instances
## self.typ =typ # Need typ_knob and etc ??
## self.i_diag=[]
@@ -765,23 +850,23 @@ def partran_end_all(file_name_in_all,file_name_out):
## # Instances for 'TRAJ'
## self.Rx_inv=None
-## self.Ry_inv=None
-
+## self.Ry_inv=None
+
## class LATTICE:
## '''
## Class to handle a TraceWin lattice file.
-
+
## - For now, only for the steerer and BPM. Extend as needed.
## - For now, only for the MEBT. Extend as needed.
## - For now, reading a file for ken. Implement a method for future ??
-
+
## 2015.04.12
## '''
## def __init__(self,path_file_lat,path_file_ken):
## # Consts, need to define globally ??
## mass=938.2723; clight=2.99792458
-
+
## # Some dict and list (MEBT and DTL are supported for now...)
## dic_class={'QUAD' : QUAD ,
## 'THIN_STEERING' : THIN_STEERING ,
@@ -801,7 +886,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## try : ken.append(float(lin.split()[2]))
## except: pass
## ken.insert(0,ken[0])
-
+
## # Go through the lat file
## with open(path_file_lat) as file:
## lst=[]; i=0; Nelem=0
@@ -822,7 +907,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## Brho =(10.0/clight)*gamma*beta*mass*1e-3; lst[-1].Brho =Brho
## # End at 'END'
## if lst[-1].typ=='END': break
-
+
## # Go through lat again and define matchings, extend/modify as needed
## match={}
## for i in range(len(lst)):
@@ -834,7 +919,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## if lst[i].typ in lst_diag:
## try : match[lst[i].fam].i_diag.append(i)
## except: pass
-
+
## # Instances
## self.lst =lst
## self.match=match
@@ -886,7 +971,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## except:
## break
-## # Index of n-th elem end, starting from 0
+## # Index of n-th elem end, starting from 0
## i_elem=[0]
## for n in range(1,Nelem[-1]):
## k=1
@@ -894,7 +979,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## try : i_elem.append(Nelem.index(n+k)-1); break
## except: k+=1
## i_elem.append(len(Nelem)-1)
-
+
## # Instances
## self.Nelem=Nelem
## self.s =s
@@ -912,7 +997,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## Remove comments, names, and etc from a TraceWin file.
## To convert the output to str, use: [' '.join(l)+'\n' for l in lat].
## The current version not capitalizing.
-
+
## Output: Lattice in a list
## 2015.03.18
@@ -928,7 +1013,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## if lin!=[]: # Avoid empty line
## lat.append(lin)
## if lin[0].upper()=='END': break
-
+
## return lat
## def get_steerer(path_file):
@@ -936,7 +1021,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## '''
## Extract steerer values from "Adjusted_Values.txt" of TraceWin.
## Note the key is the "lat elem #" where commands are not counted.
-
+
## Output is dic of {lat elem #:[val]} or {lat elem #:[val_x,val_y]}
## 2014.10.17
@@ -967,7 +1052,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## 2014.10.17
## '''
-## #--
+## #--
## i_elem=0; dic_steerer={}
## for i in range(len(lat)):
## if lat[i][0] in lst_typ_elem:
@@ -1015,7 +1100,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## for i in range(len(lat)):
## if lat[i][0]=='QUAD':
-
+
## # Adjust format
## if len(lat[i])<5:
## for k in range(5): lat[i].append('0')
@@ -1044,16 +1129,16 @@ def partran_end_all(file_name_in_all,file_name_out):
## '''
## MERGE THIS TO PARTRAN CLASS !!!!
-
+
## Extract output Twiss and halo from partran1.out
## Longitudinal phase space is z-z'.
-
+
## 2015.01.15
## '''
-
+
## file=open(path_file)
-
+
## for lin in file.readlines():
## lin=lin.split()
## if lin[0]=='####':
@@ -1063,7 +1148,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## idx_alf =[lin.index(p) for p in ("sxx'" ,"syy'" ,"szdp" )]
## idx_hal =[lin.index(p) for p in ("hx" ,"hy" ,"hp" )]
## file.close()
-
+
## gamma=float(lin[idx_gamma])+1.0
## eps =array([float(lin[idx]) for idx in idx_eps])
## bet =array([float(lin[idx])**2 for idx in idx_bet]) # Actually Sig_{i,i}
@@ -1073,7 +1158,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## bet= bet/eps*sqrt(gamma**2-1.0); bet[2]*=gamma**2 # z-dp => z-z'
## alf=-alf/eps*sqrt(gamma**2-1.0)
## gam= (1.0+alf**2)/bet
-
+
## return eps,bet,alf,gam,hal
## class PROJECT_SING_CORE:
@@ -1083,7 +1168,7 @@ def partran_end_all(file_name_in_all,file_name_out):
## (Each run should have its own project.)
## - Make sure the project/lattice file is in the calc dir.
## - Changing the calc dir option not working for tracelx64.
-
+
## 2015.10.15
## '''
@@ -1092,10 +1177,10 @@ def partran_end_all(file_name_in_all,file_name_out):
## # Make sure "path_cal" ends with "/".
## if path_cal[-1]!='/': path_cal+='/'
-
+
## # Check the calc dir exists.
## if isdir(path_cal)==False: print 'Calc dir does not exist !!!! Exiting ...'; exit()
-
+
## # Check "tracelx64" is in the calc dir.
## if isfile(path_cal+'tracelx64')==False:
## try : system('cp /opt/cea/tracelx64 '+path_cal)
@@ -1108,10 +1193,10 @@ def partran_end_all(file_name_in_all,file_name_out):
## # Check the project/lattice file is in the calc dir.
## if '/' not in file_name_proj: file_name_proj=path_cal+file_name_proj
-## if '/' not in file_name_lat : file_name_lat =path_cal+file_name_lat
+## if '/' not in file_name_lat : file_name_lat =path_cal+file_name_lat
## if isfile(file_name_proj)==False: print 'project file not in calc dir !!!! Exiting ...'; exit()
## if isfile(file_name_lat) ==False: print 'lattice file not in calc dir !!!! Exiting ...'; exit()
-
+
## # Instances (Add more options as needed.)
## self.path_cal =path_cal
## self.file_name_proj=file_name_proj
@@ -1125,5 +1210,5 @@ def partran_end_all(file_name_in_all,file_name_out):
## if self.file_name_lat!=None: opt_exe+=' dat_file=' +self.file_name_lat
## if self.seed !=None: opt_exe+=' random_seed='+self.seed
## if self.flag_hide !=None: opt_exe+=' hide'
-
+
## system(opt_exe)
diff --git a/ess/lib_tw_elem.py b/ess/lib_tw_elem.py
index fe96ea43346c8692536fc122842e238a117a6ed4..6ce3cbe01f1e67af7c9ab593aae2e9e41e5ef006 100644
--- a/ess/lib_tw_elem.py
+++ b/ess/lib_tw_elem.py
@@ -2,24 +2,24 @@
'''
- Classes for TraceWin elements/commands for external manipulations.
- First double check all the used elements/commands.
-
+
- Element/command types
* Active elements
FIELD_MAP
+ DRIFT
QUAD
THIN_STEERING
GAP
DTL_CEL
-
+
* Passive elements
-
- DRIFT
+
BEND
EDGE
- APERTURE
+ APERTURE
DIAG_POSITION
DIAG_...
@@ -40,7 +40,7 @@
LATTICE_END
SET_ADV
- SET_SIZE
+ SET_SIZE
SET_TWISS
SET_ACHROMAT
MATCH_FAM_GRAD
@@ -64,7 +64,7 @@
PLOT_DST
- 2015.10.14
+ 2016.01.19
'''
#---- Libs
@@ -107,7 +107,7 @@ class ELEM:
# Basic instances
self.name=name
self.typ =typ
- self.para=para
+ self.para=para
# Option instances
self.L = 0.0
@@ -134,11 +134,11 @@ class ELEM:
## def get_madx(self):
- ## i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
+ ## i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
## if self.name!='': lin= self.name +': marker;'
## if self.name=='': lin='elem'+i+'_marker'+': marker;'
-
+
## return lin
class COMM:
@@ -163,7 +163,7 @@ class COMM:
# Update idx
self.idx+=1
-
+
def get_tw(self):
if self.name!='': lin=self.name+': '+self.typ+' '+' '.join(self.para)
@@ -172,7 +172,7 @@ class COMM:
return lin
#---- Field map (need an additional input, dic_fmap)
-
+
class FIELD_MAP(ELEM):
'''
- Calculation of gamma only 1D (Ez) supported.
@@ -194,12 +194,12 @@ class FIELD_MAP(ELEM):
# Option instances
self.data=dic_fmap[self.name_fmap]
-
+
def update_gamma(self):
# Temp error message for field map types
if self.typ_fmap!='100': print 'Gamma calc only supported for 1D, exiting ...', exit()
-
+
# Update gamma (simple ver, closer to TW)
dL=self.data.dL; V=self.Enom*array(self.data.field)*dL; phs=self.phs_rf; C=0.0; S=0.0
for i in range(len(V)):
@@ -247,8 +247,8 @@ class FIELD_MAP(ELEM):
lin+='L=' +str(self.L )+', '
lin+='FREQ='+str(self.freq )+', '
lin+='VOLT='+str(self.E0TL )+', '
- lin+='LAG=' +str((0.5*pi-self.phs_syn)/(2.0*pi))+'; ' # phs_MADX = pi/2-phs_TW
-
+ lin+='LAG=' +str((0.5*pi-self.phs_syn)/(2.0*pi))+'; ' # phs_MADX = pi/2-phs_TW
+
return lin
def get_fluka(self):
@@ -260,18 +260,72 @@ class FIELD_MAP(ELEM):
lin+=str(self.L)+' '+str(self.s)+' '
lin+='0.0 0.0 0.0 '
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
-
+
return lin
-
+
def get_mars(self):
lin ='"RFCAVITY" "" "" '+str(self.s)+' '+str(self.L)+' '
lin+='0 0 0 0 0'
return lin
-
+
#---- Active elements
+class DRIFT(ELEM):
+ '''
+ '''
+ def __init__(self,name,typ,para):
+
+ ELEM.__init__(self,name,typ,para)
+
+ # TW instances
+ self.L =float(para[0])*1e-3 # mm => m
+ self.apt=float(para[1])
+
+ # TW option instances
+ try : self.apty=float(para[2])
+ except: self.apty=0.0
+
+ def get_tw(self):
+
+ if self.name!='': lin=self.name+': '+self.typ+' '
+ if self.name=='': lin= self.typ+' '
+ lin+=str(self.L*1e3)+' '
+ lin+=str(self.apt )+' '
+ lin+=str(self.apty )
+
+ return lin
+
+ def get_madx(self):
+
+ i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
+
+ if self.name!='': lin= self.name +': DRIFT, L='+str(self.L)+';'
+ if self.name=='': lin='ELEM'+i+'_DRIFT'+': DRIFT, L='+str(self.L)+';'
+
+ return lin
+
+ def get_fluka(self):
+
+ i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
+
+ if self.name!='': lin= self.name +' DRIFT '
+ if self.name=='': lin='ELEM'+i+'_DRIFT'+' DRIFT '
+ lin+=str(self.L)+' '+str(self.s)+' '
+ lin+='0.0 0.0 0.0 '
+ if self.apty==0.0: lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
+ if self.apty!=0.0: lin+='RECTANGLE '+str(self.apt*1e-3)+' '+str(self.apty*1e-3)
+
+ return lin
+
+ def get_mars(self):
+
+ lin ='"DRIFT" "" "" '+str(self.s)+' '+str(self.L)+' '
+ lin+='0 0 0 0 0'
+
+ return lin
+
class QUAD(ELEM):
'''
- Tilt not supported.
@@ -279,13 +333,13 @@ class QUAD(ELEM):
def __init__(self,name,typ,para):
ELEM.__init__(self,name,typ,para)
-
+
# TW instances
self.L =float(para[0])*1e-3 # mm => m
self.G =float(para[1])
self.apt =float(para[2])
self.tilt=0.0
-
+
# TW option instances
try : self.b3 =float(para[4])
except: self.b3 =0.0
@@ -297,7 +351,7 @@ class QUAD(ELEM):
except: self.b6 =0.0
try : self.r_good=float(para[8])
except: self.r_good=0.0
-
+
def get_tw(self):
if self.name!='': lin=self.name+': '+self.typ+' '
@@ -309,7 +363,7 @@ class QUAD(ELEM):
lin+=str(self.b3 )+' '
lin+=str(self.b4 )+' '
lin+=str(self.b5 )+' '
- lin+=str(self.b6 )+' '
+ lin+=str(self.b6 )+' '
lin+=str(self.r_good)
return lin
@@ -317,12 +371,12 @@ class QUAD(ELEM):
def get_madx(self):
i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
-
+
if self.name!='': lin= self.name +': QUADRUPOLE, '
if self.name=='': lin='ELEM'+i+'_QUAD'+': QUADRUPOLE, '
lin+='L=' +str(self.L )+', '
- lin+='K1='+str(self.G/Brho(self.gamma))+'; '
-
+ lin+='K1='+str(self.G/Brho(self.gamma))+'; '
+
return lin
def get_fluka(self):
@@ -334,11 +388,11 @@ class QUAD(ELEM):
lin+=str(self.L)+' '+str(self.s)+' '
lin+='0.0 0.0 '+str(self.G*self.L/Brho(self.gamma))+' '
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
-
+
return lin
-
+
def get_mars(self):
-
+
lin ='"QUADRUPOLE" "" "" '+str(self.s)+' '+str(self.L)+' '
lin+='0 '+str(self.G/Brho(self.gamma))+' 0 0 0 '
@@ -352,7 +406,7 @@ class THIN_STEERING(ELEM):
def __init__(self,name,typ,para):
ELEM.__init__(self,name,typ,para)
-
+
# TW instances
self.BLx =float(self.para[0])
self.BLy =float(self.para[1])
@@ -364,16 +418,16 @@ class THIN_STEERING(ELEM):
if self.name!='': lin=self.name+': '+self.typ+' '
if self.name=='': lin= self.typ+' '
lin+=str(self.BLx )+' '
- lin+=str(self.BLy )+' '
+ lin+=str(self.BLy )+' '
lin+=str(self.apt )+' '
lin+= self.typ_EM
return lin
-
+
def get_madx(self):
-
+
i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
-
+
if self.name!='': lin= self.name +': KICKER, L=0, '
if self.name=='': lin='ELEM'+i+'_STEERER'+': KICKER, L=0, '
lin+='HKICK='+str(-self.BLy/Brho(self.gamma))+', '
@@ -390,16 +444,16 @@ class THIN_STEERING(ELEM):
lin+=str(self.L)+' '+str(self.s)+' '
lin+='0.0 0.0 0.0 '
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
-
+
return lin
-
+
def get_mars(self):
lin ='"KICKER" "" "" '+str(self.s)+' '+str(self.L)+' '
lin+='0 0 0 0 0'
return lin
-
+
class GAP(ELEM):
'''
- Note there are E0T0L (for TW) and E0TL (for other codes).
@@ -412,7 +466,7 @@ class GAP(ELEM):
self.E0T0L =float(para[0])*1e-6 # V => MV
self.phs_rf =float(para[1])*pi/180.0 # deg => rad
self.apt =float(para[2])
- self.phs_rf_typ='0' # Relative phase only for now
+ self.phs_rf_typ= int(para[3])
# TW option instances
try : self.beta_s=float(para[4])
@@ -426,9 +480,9 @@ class GAP(ELEM):
# Option instances
self.E0TL=self.E0T0L
-
+
def update_gamma(self):
-
+
# Update gamma (and E0TL)
if self.beta_s!=0.0 and self.T0!=0.0:
gamma_c =self.gamma+0.5*self.E0TL/mass*cos(self.phs_rf)
@@ -436,7 +490,7 @@ class GAP(ELEM):
k =self.beta_s/beta_c
self.E0TL*=(self.T0+self.T1*(k-1.0)+0.5*self.T2*(k-1.0)**2)/self.T0
self.gamma+=self.E0TL/mass*cos(self.phs_rf)
-
+
def get_tw(self):
if self.name!='': lin=self.name+': '+self.typ+' '
@@ -444,7 +498,7 @@ class GAP(ELEM):
lin+=str(self.E0T0L*1e6 )+' '
lin+=str(self.phs_rf*180.0/pi)+' '
lin+=str(self.apt )+' '
- lin+= self.phs_rf_typ +' '
+ lin+=str(self.phs_rf_typ )+' '
lin+=str(self.beta_s )+' '
lin+=str(self.T0 )+' '
lin+=str(self.T1 )+' '
@@ -461,7 +515,7 @@ class GAP(ELEM):
lin+='FREQ='+str(self.freq )+', '
lin+='VOLT='+str(self.E0TL )+', '
lin+='LAG=' +str((0.5*pi-self.phs_rf)/(2.0*pi))+'; ' # phs_MADX = pi/2-phs_TW
-
+
return lin
def get_fluka(self):
@@ -473,7 +527,7 @@ class GAP(ELEM):
lin+=str(self.L)+' '+str(self.s)+' '
lin+='0.0 0.0 0.0 '
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
-
+
return lin
def get_mars(self):
@@ -501,7 +555,7 @@ class DTL_CEL(ELEM):
self.E0T0L =float(para[6])*1e-6 # V => MV
self.phs_rf =float(para[7])*pi/180.0 # deg => rad
self.apt =float(para[8])
- self.phs_rf_typ='0' # Relative phase only for now
+ self.phs_rf_typ= int(para[9])
# TW option instances
try : self.beta_s=float(para[10])
@@ -520,7 +574,7 @@ class DTL_CEL(ELEM):
# Save the gamma_ini (for MADX)
self.gamma_ini=self.gamma
-
+
# Update gamma (and E0TL)
if self.beta_s!=0.0 and self.T0!=0.0:
gamma_c =self.gamma+0.5*self.E0TL/mass*cos(self.phs_rf)
@@ -528,7 +582,7 @@ class DTL_CEL(ELEM):
k =self.beta_s/beta_c
self.E0TL*=(self.T0+self.T1*(k-1.0)+0.5*self.T2*(k-1.0)**2)/self.T0
self.gamma+=self.E0TL/mass*cos(self.phs_rf)
-
+
def get_tw(self):
if self.name!='': lin=self.name+': '+self.typ+' '
@@ -538,11 +592,11 @@ class DTL_CEL(ELEM):
lin+=str(self.L_Q2*1e3 )+' '
lin+=str(self.ds_gap*1e3 )+' '
lin+=str(self.G_Q1 )+' '
- lin+=str(self.G_Q2 )+' '
+ lin+=str(self.G_Q2 )+' '
lin+=str(self.E0T0L*1e6 )+' '
lin+=str(self.phs_rf*180.0/pi)+' '
lin+=str(self.apt )+' '
- lin+= self.phs_rf_typ +' '
+ lin+=str(self.phs_rf_typ )+' '
lin+=str(self.beta_s )+' '
lin+=str(self.T0 )+' '
lin+=str(self.T1 )+' '
@@ -566,13 +620,13 @@ class DTL_CEL(ELEM):
if self.name!='': lin+= self.name +'_GAP: RFCAVITY, L=0, '
if self.name=='': lin+='ELEM'+i+'_DTLCELL'+'_GAP: RFCAVITY, L=0, '
lin+='FREQ='+str(self.freq )+', '
- lin+='VOLT='+str(self.E0TL )+', '
+ lin+='VOLT='+str(self.E0TL )+', '
lin+='LAG=' +str((0.5*pi-self.phs_rf)/(2.0*pi))+'; \n' # phs_MADX = pi/2-phs_TW
if self.name!='': lin+= self.name +'_DRIFT2: DRIFT, '
if self.name=='': lin+='ELEM'+i+'_DTLCELL'+'_DRIFT2: DRIFT, '
lin+='L='+str(0.5*self.L-self.L_Q2+self.ds_gap)+';\n'
-
+
if self.name!='': lin+= self.name +'_QUAD2: QUADRUPOLE, '
if self.name=='': lin+='ELEM'+i+'_DTLCELL'+'_QUAD2: QUADRUPOLE, '
lin+='L=' +str(self.L_Q2 )+', '
@@ -589,7 +643,7 @@ class DTL_CEL(ELEM):
## lin+=str(self.L_Q1)+' '+str(self.s-self.L+self.L_Q1)+' '
## lin+='0.0 0.0 '+str(self.G_Q1*self.L_Q1/Brho(self.gamma_ini))+' '
## lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0\n'
-
+
## if self.name!='': lin+= self.name +'_DRIFT1 DRIFT '
## if self.name=='': lin+='ELEM'+i+'_DTLCELL'+'_DRIFT1 DRIFT '
## lin+=str(0.5*self.L-self.L_Q1-self.ds_gap)+' '+str(self.s-0.5*self.L-self.ds_gap)+' '
@@ -601,7 +655,7 @@ class DTL_CEL(ELEM):
## lin+='0.0 '+str(self.s-0.5*self.L-self.ds_gap)+' '
## lin+='0.0 0.0 0.0 '
## lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0\n'
-
+
## if self.name!='': lin+= self.name +'_DRIFT2 DRIFT '
## if self.name=='': lin+='ELEM'+i+'_DTLCELL'+'_DRIFT2 DRIFT '
## lin+=str(0.5*self.L-self.L_Q2+self.ds_gap)+' '+str(self.s-self.L_Q2)+' '
@@ -619,9 +673,9 @@ class DTL_CEL(ELEM):
lin+=str(self.L)+' '+str(self.s)+' '
lin+='0.0 0.0 0.0 '
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
-
+
return lin
-
+
def get_mars(self):
lin ='"QUADRUPOLE" "" "" '
@@ -633,7 +687,7 @@ class DTL_CEL(ELEM):
lin+=str(self.s-0.5*self.L-self.ds_gap) +' '
lin+=str(0.5*self.L-self.L_Q1-self.ds_gap)+' '
lin+='0 0 0 0 0 \n'
-
+
lin+='"RFCAVITY" "" "" '
lin+=str(self.s-0.5*self.L-self.ds_gap)+' '
lin+='0 '
@@ -642,60 +696,16 @@ class DTL_CEL(ELEM):
lin+='"DRIFT" "" "" '
lin+=str(self.s-self.L_Q2) +' '
lin+=str(0.5*self.L-self.L_Q2+self.ds_gap)+' '
- lin+='0 0 0 0 0 \n'
+ lin+='0 0 0 0 0 \n'
lin+='"QUADRUPOLE" "" "" '
lin+=str(self.s) +' '
lin+=str(self.L_Q2)+' '
lin+='0 '+str(self.G_Q2/Brho(self.gamma))+' 0 0 0 '
-
- return lin
-
-#---- Passive elements
-
-class DRIFT(ELEM):
- '''
- '''
- def __init__(self,name,typ,para):
-
- ELEM.__init__(self,name,typ,para)
-
- # TW instances
- self.L =float(para[0])*1e-3 # mm => m
- self.apt=float(para[1])
-
- # TW option instances
- try : self.apty=float(para[2])
- except: self.apty=0.0
-
- def get_madx(self):
-
- i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
-
- if self.name!='': lin= self.name +': DRIFT, L='+str(self.L)+';'
- if self.name=='': lin='ELEM'+i+'_DRIFT'+': DRIFT, L='+str(self.L)+';'
-
- return lin
-
- def get_fluka(self):
-
- i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
-
- if self.name!='': lin= self.name +' DRIFT '
- if self.name=='': lin='ELEM'+i+'_DRIFT'+' DRIFT '
- lin+=str(self.L)+' '+str(self.s)+' '
- lin+='0.0 0.0 0.0 '
- if self.apty==0.0: lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
- if self.apty!=0.0: lin+='RECTANGLE '+str(self.apt*1e-3)+' '+str(self.apty*1e-3)
return lin
-
- def get_mars(self):
- lin ='"DRIFT" "" "" '+str(self.s)+' '+str(self.L)+' '
- lin+='0 0 0 0 0'
-
- return lin
+#---- Passive elements
class BEND(ELEM):
'''
@@ -707,7 +717,7 @@ class BEND(ELEM):
def __init__(self,name,typ,para):
ELEM.__init__(self,name,typ,para)
-
+
# TW instances
self.angle=float(para[0])*pi/180.0 # deg => rad
self.rho =float(para[1])*1e-3 # mm => m
@@ -729,7 +739,7 @@ class BEND(ELEM):
self.idx +=1
self.idx_elem+=1
self.s +=self.ds # Need an overwrite because of this
-
+
def get_madx(self):
i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
@@ -742,7 +752,7 @@ class BEND(ELEM):
lin+='TILT=' +str(self.tilt )+', '
lin+='HGAP=' +str(self.gap*5e-4)+', '
lin+='FINT=' +str(self.ff_coef1)+'; '
-
+
return lin
def get_fluka(self):
@@ -763,7 +773,7 @@ class BEND(ELEM):
lin+=str(self.angle)+' 0 0 0 '+str(self.tilt)
return lin
-
+
class EDGE(ELEM):
'''
- For MARS, assuming rbend and the export is commented out.
@@ -771,33 +781,33 @@ class EDGE(ELEM):
def __init__(self,name,typ,para):
ELEM.__init__(self,name,typ,para)
-
+
# TW instances
self.angle =float(para[0])*pi/180.0 # deg => rad
- self.rho =float(para[1])*1e-3 # mm => m
+ self.rho =float(para[1])*1e-3 # mm => m
self.gap =float(para[2])
self.ff_coef1=float(para[3])
- self.ff_coef2=float(para[4])
+ self.ff_coef2=float(para[4])
self.apt =float(para[5])
self.plane = para[6]
# Option instances
if self.plane=='0': self.tilt=0.0
if self.plane=='1': self.tilt=0.5*pi
-
+
def get_madx(self):
- i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
-
+ i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
+
if self.name!='': lin= self.name +': DIPEDGE, '
if self.name=='': lin='ELEM'+i+'_EDGE'+': DIPEDGE, '
lin+='E1=' +str(self.angle) +', ' # sbend
#lin+='E1=' +'0' +', ' # rbend
lin+='H=' +str(1.0/self.rho) +', '
- lin+='TILT='+str(self.tilt )+', '
- lin+='HGAP='+str(self.gap*5e-4)+', '
+ lin+='TILT='+str(self.tilt )+', '
+ lin+='HGAP='+str(self.gap*5e-4)+', '
lin+='FINT='+str(self.ff_coef1)+'; '
-
+
return lin
class APERTURE(ELEM):
@@ -812,21 +822,21 @@ class APERTURE(ELEM):
self.apt =float(para[0])
self.apty =float(para[1])
self.typ_apt= para[2]
-
+
def get_madx(self):
i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
if self.name!='': lin= self.name +': DRIFT, L=0;'
if self.name=='': lin='ELEM'+i+'_APT'+': DRIFT, L=0;'
-
+
## if self.name!='': lin= self.name +': collimator, l=0, '
## if self.name=='': lin='ELEM'+i+'_apt'+': collimator, l=0, '
## if self.typ_apt=='0':
## lin+='apertype=rectangle, aperture={'+str(self.apt*1e-3)+','+str(self.apty*1e-3)+'};'
## if self.typ_apt=='1':
## lin+='apertype=CIRCLE, aperture={'+str(self.apt*1e-3)+'};'
-
+
return lin
def get_fluka(self):
@@ -837,10 +847,10 @@ class APERTURE(ELEM):
if self.name=='': lin='ELEM'+i+'_APT'+' DRIFT '
lin+=str(self.L)+' '+str(self.s)+' '
lin+='0.0 0.0 0.0 '
- if self.typ_apt=='0': lin+='RECTANGLE '+str(self.apt*1e-3)+' '+str(self.apty*1e-3)
- if self.typ_apt=='1': lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
+ if self.typ_apt=='0': lin+='RECTANGLE '+str(self.apt*1e-3)+' '+str(self.apty*1e-3)
+ if self.typ_apt=='1': lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
- return lin
+ return lin
def get_mars(self):
@@ -848,7 +858,7 @@ class APERTURE(ELEM):
lin+='0 0 0 0 0'
return lin
-
+
class DIAG_POSITION(ELEM):
'''
'''
@@ -858,11 +868,11 @@ class DIAG_POSITION(ELEM):
def get_madx(self):
- i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
+ i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
if self.name!='': lin= self.name +': MONITOR, L=0;'
if self.name=='': lin='ELEM'+i+'_BPM'+': MONITOR, L=0;'
-
+
return lin
def get_fluka(self):
@@ -875,13 +885,13 @@ class DIAG_POSITION(ELEM):
lin+='0.0 0.0 0.0 '
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
- return lin
+ return lin
def get_mars(self):
lin ='"MONITOR" "" "" '+str(self.s)+' '+str(self.L)+' '
lin+='0 0 0 0 0'
-
+
return lin
class DIAG(ELEM):
@@ -893,11 +903,11 @@ class DIAG(ELEM):
def get_madx(self):
- i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
+ i='0'*(3-int(log10(self.idx_elem+1)))+str(self.idx_elem+1)
if self.name!='': lin= self.name +': INSTRUMENT, L=0;'
if self.name=='': lin='ELEM'+i+'_DIAG'+': INSTRUMENT, L=0;'
-
+
return lin
def get_fluka(self):
@@ -911,16 +921,16 @@ class DIAG(ELEM):
lin+='CIRCLE '+str(self.apt*1e-3)+' 0.0'
return lin
-
+
def get_mars(self):
lin ='"INSTRUMENT" "" "" '+str(self.s)+' '+str(self.L)+' '
lin+='0 0 0 0 0'
-
+
return lin
#---- Active commands
-
+
class STEERER(COMM):
'''
- Only magnetic steerers supported.
@@ -929,7 +939,7 @@ class STEERER(COMM):
def __init__(self,name,typ,para):
COMM.__init__(self,name,typ,para)
-
+
# TW instances
self.Bx=float(para[0])
self.By=float(para[1])
@@ -948,7 +958,7 @@ class STEERER(COMM):
lin+=str(self.By )+' '
lin+=str(self.max )+' '
lin+='0' +' '
- lin+='0' +' '
+ lin+='0' +' '
lin+=str(self.nonlin)
return lin
@@ -980,17 +990,17 @@ class CHOPPER(COMM):
return lin
#---- Passive commands
-
+
class FREQ(COMM):
'''
'''
def __init__(self,name,typ,para):
COMM.__init__(self,name,typ,para)
-
+
# TW instances
self.freq_new=float(para[0])
-
+
def update_idx(self):
self.idx +=1
@@ -1007,7 +1017,7 @@ class MARKER(COMM):
if self.name!='': lin=self.name +': MARKER;'
if self.name=='': lin=self.para[0]+': MARKER;'
-
+
return lin
def get_fluka(self):
@@ -1017,14 +1027,14 @@ class MARKER(COMM):
lin+='0.0 '+str(self.s)
return lin
-
+
## def get_mars(self):
## lin ='"marker" "" "" '+str(self.s)+' 0.0 '
## lin+='0 0 0 0 0'
## return lin
-
+
#---- No longer needed ??
## class ADJUST(COMM):
@@ -1033,10 +1043,10 @@ class MARKER(COMM):
## def __init__(self,name,typ,para):
## COMM.__init__(self,name,typ,para)
-
+
## # TW instances
## self.fam= para[0]
-## self.var=int(para[1])
+## self.var=int(para[1])
## # TW option instances
## try : self.link=self.para[2]