diff --git a/.idea/Beamline Editor.iml b/.idea/Beamline Editor.iml index 554af79..7f2553c 100644 --- a/.idea/Beamline Editor.iml +++ b/.idea/Beamline Editor.iml @@ -1,7 +1,9 @@ - + + + diff --git a/beamline_editor/table_utils.py b/beamline_editor/table_utils.py new file mode 100644 index 0000000..fa8515b --- /dev/null +++ b/beamline_editor/table_utils.py @@ -0,0 +1,139 @@ +""" +table_utils.py — Utility helpers for populating Qt table widgets. + +populate_table(table, records) + Fill a QTableWidget from a list of dicts. + + Columns (fixed, in this order): + Domain, Prefix, Suffix, Index, s (m), z (m), x (m), y (m), L (m), + Reference, Group, Baugruppe, Variante, Roll, PS-Ch. + + "Index" is always rendered as a zero-padded 3-digit integer (e.g. 7 → "007"). + + Row background colour is determined by the "Group" value: + "Section" → light grey (#D0D0D0) + "Magnet" / "Magnets" → orange (#FFAA55) + "RF" → teal (#66BBAA) + "Diagnostics" → yellow (#FFEE88) + "Vacuum" → light green (#90EE90) + "ID" → light purple (#CC99FF) + anything else → white (#FFFFFF) + + Font colour is always black. + + Any key absent from a dict produces a blank cell; existing table + content is replaced. +""" + +from __future__ import annotations +from typing import Any + +from PyQt5.QtWidgets import QTableWidget, QTableWidgetItem, QHeaderView +from PyQt5.QtGui import QColor, QBrush +from PyQt5.QtCore import Qt + + +# ── Column definition ───────────────────────────────────────────────────────── + +COLUMNS: list[str] = [ + "Domain", + "Prefix", + "Suffix", + "Index", + "s (m)", + "z (m)", + "x (m)", + "y (m)", + "L (m)", + "Reference", + "Group", + "Baugruppe", + "Variante", + "Roll", + "PS-Ch.", +] + +# ── Row colour map (Group value → background hex) ───────────────────────────── + +_GROUP_COLORS: dict[str, str] = { + "Section": "#D0D0D0", # light grey + "Magnet": "#FFAA55", # orange + "Magnets": "#FFAA55", # orange (alternate spelling) + "RF": "#66BBAA", # teal + "Diagnostics": "#FFEE88", # yellow + "Vacuum": "#90EE90", # light green + "ID": "#CC99FF", # light purple +} +_DEFAULT_COLOR = "#FFFFFF" # white +_FONT_COLOR = QColor("#000000") # black text for all rows + + +# ── Formatting helpers ──────────────────────────────────────────────────────── + +_POSITION_COLS = {"s (m)", "z (m)", "x (m)", "y (m)", "L (m)"} + + +def _format_value(key: str, value: Any) -> str: + """Return the display string for a single cell value.""" + if key == "Index": + if value is None or value == "": + return "" + try: + return f"{int(value):03d}" + except (ValueError, TypeError): + return str(value) + if key in _POSITION_COLS: + if value is None or value == "": + return "" + try: + return f"{float(value):.4f}" + except (ValueError, TypeError): + return str(value) + return "" if value is None else str(value) + + +# ── Public function ──────────────────────────────────────────────────────────── + +def populate_table( + table: QTableWidget, + records: list[dict[str, Any]], +) -> None: + """ + Populate *table* from *records*. + + Parameters + ---------- + table : QTableWidget + The widget to fill. Column count and headers are set + automatically; any previous content is cleared. + records : list[dict] + Each dict may contain any subset of the known column keys. + Missing keys produce an empty cell. Extra keys are ignored. + """ + # ── setup ───────────────────────────────────────────────────────────────── + table.setUpdatesEnabled(False) + table.clearContents() + table.setRowCount(len(records)) + table.setColumnCount(len(COLUMNS)) + table.setHorizontalHeaderLabels(COLUMNS) + + header = table.horizontalHeader() + header.setSectionResizeMode(QHeaderView.ResizeToContents) + header.setStretchLastSection(True) + + # ── fill rows ───────────────────────────────────────────────────────────── + for row, record in enumerate(records): + group_value = record.get("Group", "") + bg_hex = _GROUP_COLORS.get(group_value, _DEFAULT_COLOR) + bg_brush = QBrush(QColor(bg_hex)) + + for col, key in enumerate(COLUMNS): + raw = record.get(key) # None when key absent + text = _format_value(key, raw) + cell = QTableWidgetItem(text) + cell.setFlags(Qt.ItemIsEnabled | Qt.ItemIsSelectable) + cell.setBackground(bg_brush) + cell.setForeground(QBrush(_FONT_COLOR)) + table.setItem(row, col, cell) + + table.setUpdatesEnabled(True) diff --git a/containers.py b/containers.py index a681a80..9edf450 100644 --- a/containers.py +++ b/containers.py @@ -164,7 +164,7 @@ class LineContainer: for ele in self.Element: ds = self.position(ele) - Last app.writeDrift(ds) - name = ele.writeElement(app, EM) + name = ele.writeElement(app) Last = self.position(ele) + ele.getResLength() # if name<>0: # not valid in python3.x if name != 0: @@ -175,7 +175,7 @@ class LineContainer: return baseseq else: for ele in self.Element: - name = ele.writeLattice(app, EM) + name = ele.writeLattice(app) # if name<>0: # not valid in python3.x if name != 0: seq.append(name) diff --git a/elements.py b/elements.py index 30e51e1..7fffb78 100644 --- a/elements.py +++ b/elements.py @@ -50,7 +50,9 @@ class SimpleContainer: prop.update({'DisplacementY': 0.}) self.__dict__.update(prop) + self.adjustLength() + def adjustLength(self): if self.LengthRes <= 0: self.LengthRes = self.Length self.sRef = 0. @@ -70,6 +72,7 @@ class SimpleContainer: def getResLength(self): # return the reserved length for that element + self.adjustLength() return self.LengthRes def getLength(self): @@ -302,7 +305,7 @@ class Diagnostic(SimpleContainer): if not ('enable' in prop): # prop.has_key('enable')): self.__dict__.update({'enable': 0}) - def writeElement(self, app, EM): + def writeElement(self, app): app.writeDiagnostic(self) return self.Name diff --git a/main.py b/main.py index d9443a4..978a744 100644 --- a/main.py +++ b/main.py @@ -7,9 +7,12 @@ from PyQt5 import QtWidgets from ui.BeamlineGUI import Ui_BeamlineGUI from beamline_editor import BeamlineEditorWidget +from beamline_editor.table_utils import populate_table + from elements import * from containers import * from propertytypes import getPropertyList +from protoliste import ProtoListe class BeamlineEditor(QtWidgets.QMainWindow, Ui_BeamlineGUI): def __init__(self): @@ -44,6 +47,9 @@ class BeamlineEditor(QtWidgets.QMainWindow, Ui_BeamlineGUI): self.savefile='' self.lines={} + self.PL=ProtoListe(0) + + def saveas(self): options = QtWidgets.QFileDialog.Options() @@ -114,7 +120,7 @@ class BeamlineEditor(QtWidgets.QMainWindow, Ui_BeamlineGUI): case "Reference": ele = Marker({"Tag":"QREF"}) case "Start": - ele = Marker({"Tag":"QBEG"}) + ele = Marker({"Tag":"QBEG","Group":"Start"}) case "Branch": ele = Marker({"Tag":"MKBR"}) case "Stop": @@ -186,41 +192,46 @@ class BeamlineEditor(QtWidgets.QMainWindow, Ui_BeamlineGUI): self.lines.clear() type = self.editor.getElementType(node_id) if type == 'Line': - name = self.elementDB[node_id].Name - resL = self.elementDB[node_id].LengthEnding - elif type == 'Start': - name = self.elementDB[node_id].Name - resL=0 + line = self.elementDB[node_id] else: - name='' - resL=0 - if len(name)<1: - name='XXX' - print('Flattening Beamline %s ...' % name) - self.lines[node_id]=LineContainer(namein=name,Lin=resL) - self.unwrap(name,node_id,node_id) + if type == 'Start': + name = self.elementDB[node_id].Name + if len(name) < 1: + name='XXX' + line=LineContainer(namein=name,Lin=0) + print('\nUnwrapping beamline...\n') + self.unwrap(name,node_id,line) + self.PL.generateLayout(line) + plrecs = [self.PL.info[ele] for ele in self.PL.order] + populate_table(self.UIProtoList,plrecs) - def unwrap(self, name, node_id,line_id): + + def unwrap(self, name, node_id,line): print('Unwrapping Line: %s' % name) tree=self.editor.flatten(node_id) for ele in tree.path: + self.report(ele, name) + sRef=self.elementDB[ele].OffsetS + Ref='absolute' + if self.elementDB[ele].relOff: + Ref='relative' + line.append(self.elementDB[ele],sRef=sRef,Ref=Ref) + self.elementDB[ele].Prefix=name type = self.editor.getElementType(ele) if type == 'Line': - name1 = self.elementDB[ele].Name - resL = self.elementDB[ele].LengthEnding - ref = self.elementDB[ele].firstElementID - + subline = self.elementDB[ele] + name1 = subline.Name + ref=subline.firstElementID if ref in self.elementDB.keys(): - self.lines[ele] = LineContainer(namein=name+name1, Lin=resL) - self.unwrap(name+name1, ref, ele) + self.unwrap(name+name1, ref, subline) else: print('Undefined reference in Line Container:', name+name1) else: - self.report(ele, name) - self.lines[line_id].append(self.elementDB[ele]) + self.elementDB[ele].prefix=name + if tree.primary: - self.unwrap(name,tree.primary.node_id,line_id) + self.unwrap(name,tree.primary.node_id,line) @@ -228,11 +239,17 @@ class BeamlineEditor(QtWidgets.QMainWindow, Ui_BeamlineGUI): def report(self,node_id,prefix=''): type = self.editor.getElementType(node_id) name = self.elementDB[node_id].Name - if len(prefix)>0: - name = prefix+'-'+name + if type=='Start': name=prefix - print('FLATTEN: Adding %s Element with Name %s' % ( type,name)) + elif type == 'Line': + name=prefix+name + else: + name = prefix+'-'+name + if type == 'Line': + print('FLATTEN: Adding Line Container: %s' % name) + else: + print('FLATTEN: Adding %s Element with Name %s' % ( type,name)) if __name__ == '__main__': diff --git a/onlinemodel/core/apptemplate.py b/onlinemodel/core/apptemplate.py new file mode 100644 index 0000000..5ee9859 --- /dev/null +++ b/onlinemodel/core/apptemplate.py @@ -0,0 +1,59 @@ +class ApplicationTemplate: + def __init__(self): + self.switch=0 + self.path= '..' + self.avoidPreset=0 + self.MapIndx=[] + self.MapIndxSave=[] + + def isType(self,key): + return 0 + + def demandMapID(self): + return 0 + + + def writeLine(self,line,seq): + return + + def writeDrift(self,ele): + return + + def writeVacuum(self,ele): + self.writeMarker(ele) + return + + def writeAlignment(self,ele): + self.writeMarker(ele) + return + + def writeBend(self,ele): + return + + def writeQuadrupole(self,ele): + return + + def writeCorrector(self,ele): + return + + + def writeSextupole(self,ele): + return + + def writeRF(self,ele): + return + + + def writeUndulator(self,ele): + return + + def writeDiagnostic(self,ele): + return + + + def writeMarker(self,ele): + return + + def writeSolenoid(self,ele): + return + diff --git a/onlinemodel/core/facility.py b/onlinemodel/core/facility.py new file mode 100644 index 0000000..4f462de --- /dev/null +++ b/onlinemodel/core/facility.py @@ -0,0 +1,376 @@ +from copy import deepcopy +from string import * +import math +import re +import numpy as np +import sys + + +from onlinemodel.core.facilityaccess import FacilityAccess, LineOrder +#from onlinemodel.layout.layout import SwissFEL +#from onlinemodel.core.energymanager import EnergyManager +#from onlinemodel.core.omtype import * + +class Facility(FacilityAccess): + def __init__(self,init=1,alt=0,layout=None): + + # init=0 : no initialization of dynamic parameters + # init=1 : initialization of dynamic parameters as in OMLayout.initialize() + + # Note that FacilityContainer.BeamPath() set the branch bend to the design angle + + # alt = 0 : select current layout in OMLayout + # alt <> 0: select alternative Layout in OMLayout + + # layout = None : loads the curren tlayout + # layout = '10.3.2': loads SwissFEL from the module file OMLayout_10_3_2.py + + Layout=SwissFEL(alt) # instanciate the layout + if layout is not None and Layout.Version is not layout: + module_name='OMLayout_%s' % layout.replace('.','_') + module=__import__(module_name) + Layout=module.SwissFEL(alt) + print('Importing Layout Module:',module_name) + self.Version=Layout.Version + print('Initializing Lattice',self.Version,flush=True); + + self.PartsList=Layout.build() # generate the lattice and return the beamlines + self.EM=EnergyManager() + self.branch = None + + # Note that the subsections have also a direct reference to the elements but if elements are altered + # in ElementDB then the SectionDB is automatically updated with the new values. + self.ElementDB={} # dictionary to hold all elements of the facility + self.SectionDB={} # dictionary of subsections with 7 Letter ID + self.IntegratedCorrDB={} # Database for integrated correctors + + # flatten beamlines so that only one outer container exists with the section container, e.g. SINSB01 + + for p in self.PartsList: + line=LineContainer('') # + p[0].flatten('S',line,self.ElementDB,self.SectionDB,self.IntegratedCorrDB) + p[0]=line # replace the current nested line container with the list of atomic subsection + + # initiate the element values with default values + if init: + Layout.initialize(self.ElementDB) + + # define the mapping of the element, defining their position and asignment to a branch + self.mapping() + + # setup a dictionary of bending magnets. + BCdict={} + self.BranchBend=[] + for k in self.ElementDB.keys(): + if isinstance(self.ElementDB[k],Dipole): + try: + self.ElementDB[k].branch +# print(k) + self.BranchBend.append(k) + except: + None + try: + self.ElementDB[k].BC + if self.ElementDB[k].BC in BCdict: + # It must be the element itself (not Element.Name) becuase the new feature altSection allows the same element name in the altered section... + BCdict[self.ElementDB[k].BC].append(self.ElementDB[k]) + else: + BCdict[self.ElementDB[k].BC]=[self.ElementDB[k]] + except: + None + + + # get braching points and bunch compressor magnets + + self.BC=[] + self.BCBPM={} + for k in BCdict.keys(): + BC=BCdict[k] + BCsorted=[0]*len(BC) + for i in range(0,len(BC)): + judge=[] + for j in range(0,len(BC)): + if i!=j: + judge.append(self.isUpstream(BC[i].Name,BC[j].Name)) + jsum=sum(judge) + BCsorted[len(BC)-jsum-1]=BC[i] + + + + if len(BCsorted)==4: + sect=BCsorted[3].Name.split('.')[0] + for s in self.SectionDB.keys(): + if s==sect: + #if BCsorted[3].altSection==s.altSection: # altSection is turned off for now... 29.05.2015 + #arm=self.SectionDB[s].position(self.ElementDB[BCsorted[3].Name])-self.SectionDB[s].position(self.ElementDB[BCsorted[2].Name]) + arm=self.SectionDB[s].position(BCsorted[3])-self.SectionDB[s].position(BCsorted[2]) + BCsorted.append(arm) + self.BCBPM[BCsorted[0].Name]=['','',0] + for ele in self.SectionDB[s].Element: + if 'DBPM' in ele.Name: # This does not detect Laser heater... + if self.isUpstream(BCsorted[0].Name,ele.Name) and self.isUpstream(ele.Name,BCsorted[1].Name): + self.BCBPM[BCsorted[0].Name][0]=ele # BPM between the first and second dipoles + if self.isUpstream(BCsorted[2].Name,ele.Name) and self.isUpstream(ele.Name,BCsorted[3].Name): + self.BCBPM[BCsorted[0].Name][1]=ele # BPM between the third and fourth dipoles + self.BCBPM[BCsorted[0].Name][2]=self.SectionDB[s].position(BCsorted[3])+BCsorted[3].Length/2-self.SectionDB[s].position(ele)-ele.Length/2 + + self.BC.append(BCsorted) + + + # get list of first nd last magnetin bunch compressor. + self.BCstart={} + self.BCend={} + for ele in self.BC: + self.BCstart[ele[0]]=1 + self.BCend[ele[-1]]=1 + + + #initialize all element with the energy. All elements need to be registered first to avoid an error message of not registered elements + for ele in self.ElementDB.values(): + self.EM.register(ele) + self.writeFacility(self.EM) + self.forceEnergyAt('SINLH01.DBAM010', 140e6) # set reference energy to 140 MeV at Laser Heater + + + def turnOffBends(self): + for ele in self.ElementDB.values(): + if isinstance(ele,Dipole): + ele.angle=0.0 + + def turnOffBranchBends(self): + for ele in self.ElementDB.values(): + if isinstance(ele,Dipole) and ('branch' in ele.__dict__): + ele.angle=0.0 + +# end of initialization +#---------------------------------------------------------------------------------- + + def exportBranch(self,exporter): + if self.branch: + self.branch.writeLattice(exporter,self) + + def setBranchingMagnets(self,dest): + # update fixed angle of bending magnets for the different paths + if 'SAT' in dest: + ang = self.getRegExpElement('S20SY02', 'MKAC', 'design_kick') + self.setRegExpElement('S20SY02', 'MKAC0.0', 'cory', ang[0]) + ang = self.getRegExpElement('S20SY02', 'MKDC', 'design_kick') + self.setRegExpElement('S20SY02', 'MKDC0.0', 'cory', ang[0]) + self.setRegExpElement('S20SY02', 'MBND', 'angle', 2) + else: + self.setRegExpElement('S20SY02', 'MK.C0.0', 'cory', 0) + self.setRegExpElement('S20SY02', 'MBND', 'angle', 0) + if 'SPO' in dest: + ang = self.getRegExpElement('S30CB15', 'MKAC', 'design_kick') + if len(ang)> 0: + self.setRegExpElement('S30CB15', 'MKAC0.0', 'cory', ang[0]) + ang = self.getRegExpElement('S30CB15', 'MKDC', 'design_kick') + self.setRegExpElement('S30CB15', 'MKDC0.0', 'cory', ang[0]) + print('### setting angle for Porthos') + self.setRegExpElement('S30SY01', 'MBND', 'angle', -2.55) + else: + self.setRegExpElement('S30CB15', 'MK.C0.0', 'cory', 0) + self.setRegExpElement('S30SY01', 'MBND', 'angle', 0) + if 'S10' in dest: + self.setRegExpElement('S10DI01', 'MBND', 'angle', 20) + else: + self.setRegExpElement('S10DI01', 'MBND', 'angle', 0) + + + def setBranch(self, destination,start='SINLH01',end=None): + destname = destination.upper() + # default branch is SwissFEL Injector + dest = destname + if 'ARAMIS' in destname: + dest = 'SARBD02' + elif 'ATHOS' in destname: + dest = 'SATBD02' + elif 'PORTHOS' in destname : + dest = 'SPOSY02' + if end is None: + end = 'SPOSY02' + elif 'INJECTOR' in destname: + dest = 'S10BD01' + + self.setBranchingMagnets(dest) + sec = self.getSection(dest) + path = sec.mapping + self.branch = self.BeamPath(path) + if end is None: + self.branch.setRange(start) + else: + self.branch.setRange(start,end) + + def getBranchElements(self,filter='.*'): + order=LineOrder(filter) + self.exportBranch(order) + return order.order + + + def isUpstream(self,a,b): + # Tell you if a is upstream of b. + try: + A=self.ElementDB[a] + Amap=A.mapping + except: + try: + A=self.IntegratedCorrDB[a] + Amap=A.mapping + except: + print('Your input element name is not in Element data base...A:',a,flush=True) + return None + + try: + B=self.ElementDB[b] + Bmap=B.mapping + except: + try: + B=self.IntegratedCorrDB[b] + Bmap=B.mapping + except: + print('Your input element name is not in Element data base...B:',b,flush=True) + return None + + if Amap[0]==Bmap[0]: # Two elements are in the same part of the machine, easy case + if Amap[2]=0: + parent=self.PartsList[myself][1] + if Amap[0]==parent: + branch=self.PartsList[myself][2] + if branch==-1: + return True # The part is connected to the paprent at the end, so Elem A is upstream of B + elif Amap[1]<=branch: + return True # The part is connected to the parent at a branch and Elem A is upstream of the branch, and thus upstream of B + elif Amap[1]>branch: + return False # The part is connected to the parent at a branch but Elem A is downstream of the branch and thus not upstream of B + myself=parent + + return False # The part which includes Elem A is not upstream of Elem B + + def mapping(self): + + foundMarker={} + List1=[] + for i in range(0,len(self.PartsList)): + p=self.PartsList[i] + found=0 + #line=self.BeamPath(i) + j=1 + k=0 + for subsec in p[0].Element: + for ele in subsec.Element: + if ele.Tag=='MKBR': + key=str(i)+'_'+str(j) + found=found+1 + foundMarker[key]=ele.Name + j=j+1 + k=k+1 + ele.__dict__.update({'mapping':[i,j,k]}) + #if ele.__dict__.has_key('cor') or ele.__dict__.has_key('corx') or ele.__dict__.has_key('cor'): + if 'cor' in ele.__dict__ or 'corx' in ele.__dict__ or 'cor' in ele.__dict__: + if ('MQUA' in ele.Name) or ('MBND' in ele.Name) or ('UIND' in ele.Name): + #if ele.__dict__.has_key('cor') or ele.__dict__.has_key('corx'): + if 'cor' in ele.__dict__ or 'corx' in ele.__dict__: + cname=ele.Name.replace('MQUA','MCRX').replace('MBND','MCRX').replace('UIND','MCRX') + cele=self.IntegratedCorrDB[cname] + cele.__dict__.update({'mapping':[i,j,k]}) + #if ele.__dict__.has_key('cory'): + if 'cory' in ele.__dict__: + cname=ele.Name.replace('MQUA','MCRY').replace('MBND','MCRY').replace('UIND','MCRY') + cele=self.IntegratedCorrDB[cname] + cele.__dict__.update({'mapping':[i,j,k]}) + + section=self.SectionDB[subsec.Name] + section.__dict__.update({'mapping':i}) + List1.append([p[1],p[2],found]) + + self.foundMarker=foundMarker + + return List1 + + + def BeamPath(self, PathIndx, Branch=0): + + #Branch indicates that for a branching marker is searched and then stopped at the correct occurence (e.g. Branch=2 stops at the second marker) + if PathIndx<0 or PathIndx>len(self.PartsList)-1: + return None + if not isinstance(PathIndx,int): + return None + + path=self.PartsList[PathIndx] + parentID=path[1] # reference to parent beamline (0 means root and has no parent) + branchID=path[2] # branching point: the nth occurence of a marker in the parent beam line. If negative -1 than line is appended to parent line + line=LineContainer('') # clear current line + + if (branchID!=0): # get line from parents first + line=self.BeamPath(parentID,branchID) + + foundMarker=0 # to count makers in a given line + for sec in path[0].Element: # loob over all subsection in a line + for ele in sec.Element: # check if the element of the subsection is a marker and increase counter if so + #if isinstance(ele,Dipole): + # ele.angle=ele.design_angle + if isinstance(ele,Marker) and 'MKBR' in ele.Name: + foundMarker=foundMarker+1 + if Branch>0 and foundMarker==Branch: # check whether the branching point is found + found=0 # find dipole before marker to sert to desing_angle + idx=sec.Element.index(ele)-1 + print('### branching marker:', ele.Name) + while found==0 and idx >-1: + if isinstance(sec.Element[idx],Dipole): + print('### setting branching angle for',sec.Element[idx].Name,'for branch', Branch) + sec.Element[idx].angle=sec.Element[idx].design_angle; + found=1; + idx=idx-1; + secnew=LineContainer(sec.Name) # create a new line container to hold reduced section + for i in range(0,sec.Element.index(ele)): + ref=sec.Ref[sec.Element[i]] + secnew.append(sec.Element[i],ref['pos'],ref['ref']) # add element and its reference + line.append(secnew,0,'relative') + return line + else: + # check for dipoles to be set straight if not a branching point. + found=0; # search backwards from the marker + idx=sec.Element.index(ele)-1 + while found==0 and idx >-1: + if isinstance(sec.Element[idx],Dipole): + sec.Element[idx].angle=0; + found=1; + idx=idx-1; + line.append(sec,0,'relative') # append subsection + return line + + def writeFacility(self,app=None,simple=0,EM=None): + # Handle the entire facility at one time - experimental + # Particular application is the energy manager + + if not app: + 'no application, no action' + return + + for i in range(0,len(self.PartsList)): + p=self.PartsList[i] + if 'demandMapID' in dir(app): + if app.demandMapID(): # app is asking which part is under writing + if i==len(self.PartsList)-1: + app.MapIndx=-i # Give a negative index for the last part + else: + app.MapIndx=i + if simple: + line=p[0] # this might be wrong if a branching has happen, because it takes the full cenn and not the adjusted on eup to branch + else: + line=self.BeamPath(self.PartsList.index(p)) + print('Writing Beam line for: ', p[3],flush=True) + # if no energy manager is specified than use the internal one as default + if not EM: + EM=self + line.writeLattice(app,EM) + + + diff --git a/onlinemodel/core/facilityaccess.py b/onlinemodel/core/facilityaccess.py new file mode 100644 index 0000000..7226c9f --- /dev/null +++ b/onlinemodel/core/facilityaccess.py @@ -0,0 +1,337 @@ +import re +import numpy as np +from onlinemodel.core.apptemplate import ApplicationTemplate + +class LineOrder(ApplicationTemplate): + def __init__(self,RE='.*'): + self.order = [] + self.Re=re.compile(RE) + + def register(self,name): + name=name.replace('.','-') + if not name in self.order: + if self.Re.match(name): + self.order.append(name) + + def writeQuadrupole(self,ele): + self.register(ele.Name) + if 'corx' in ele.__dict__: + self.register(ele.Name.replace('MQUA','MCRX')) + if 'cory' in ele.__dict__: + self.register(ele.Name.replace('MQUA','MCRY')) + + def writeBend(self,ele): + self.register(ele.Name) + if 'cor' in ele.__dict__: + if ele.Tilt == 0: + self.register(ele.Name.replace('MBND','MCRX')) + else: + self.register(ele.Name.replace('MBND','MCRY')) + + + def writeCorrector(self,ele): + if 'corx' in ele.__dict__: + self.register(ele.Name.replace('MCOR','MCRX')) + if 'cory' in ele.__dict__: + self.register(ele.Name.replace('MCOR','MCRY')) + + def writeSextupole(self,ele): + self.register(ele.Name) + + def writeRF(self,ele): + self.register(ele.Name) + + def writeUndulator(self,ele): + self.register(ele.Name) + + def writeDiagnostic(self,ele): + self.register(ele.Name) + + def writeSolenoid(self,ele): + self.register(ele.Name) + + + + + + +class FacilityAccess: + def listElement(self, word=None, NoE=0): + # Return a list of elements that includes word + # Wild card (*) can be used + # NoE stands for Name or Element, switching the output + # 0: list of Element names, 1: list of element instances. + if word==None: + return + elif word=='*': + allRequested=1 + elif len(word)==0: + return [] + else: + word=word.upper() + word=word.split('*') + allRequested=0 + + ElementList=[] + for p in self.PartsList: + for sec in p[0].Element: + for ele in sec.Element: + requested=1 + for w in word: + if w not in ele.Name: + requested=allRequested + if requested: + if NoE: + ElementList.append(ele) + else: + if ele.Type!='Corrector': + ElementList.append(ele.Name) + if ('cor' in ele.__dict__) or ('corx' in ele.__dict__): + ElementList.append(ele.Name.replace(ele.Name[8:12],'MCRX')) + if ('cory' in ele.__dict__): + if (ele.Name[8:12]=='MKAC') or (ele.Name[8:12]=='MKDC'): + ElementList.append(ele.Name) + else: + ElementList.append(ele.Name.replace(ele.Name[8:12],'MCRY')) + + return ElementList + + def listSection(self, word=None, NoE=0): + # Return a list of elements that includes word + # Wild card (*) can be used + # NoE stands for Name or Element, switching the output + # 0: list of Element names, 1: list of line instances. + if word==None: + return + elif word=='*': + allRequested=1 + elif len(word)==0: + return [] + else: + word=upper(word) + word=word.split('*') + allRequested=0 + + SectionList=[] + for p in self.PartsList: + for sec in p[0].Element: + requested=1 + for w in word: + if w not in sec.Name: + requested=allRequested + if requested: + if NoE: + SectionList.append(sec) + else: + SectionList.append(sec.Name) + + return SectionList + + + + + def getElement(self,name): + name=name.upper() + if name in self.ElementDB: + return self.ElementDB[name] + else: + return None + + def getSection(self,name): + name=name.upper() + if name in self.SectionDB: + return self.SectionDB[name] + else: + return None + + + + def setGroup(self,dn,field,values): + # dn is a list of element names. It can be obtained from OMEpicsChannel + for i in range(0,len(dn)): + n=dn[i] + ele=self.ElementDB[n] + if field in dir(ele): + ele.__dict__.update({field:values[i]}) + else: + print('Error due to unregistered Channel:',ele,field) + + + def getGroup(self,dn,field): + values=[] + for i in range(0,len(dn)): + n=dn[i] + ele=self.ElementDB[n] + if field in dir(ele): + values.append(ele[field]) + else: + print('Error due to unregistered Channel:',ele,field) + + return values + + + + def setRegExpElement(self,section,device,field,value): + myRe=re.compile(section+'.*'+device) + for key in self.ElementDB.keys(): + if (myRe.match(key)): + ele=self.ElementDB[key] + if field in ele.__dict__.keys(): + if (field=='angle'): + ele.__dict__[field]=np.sign(ele['design_angle'])*abs(value) # this is needed to set collectively a bunch compressor angle without flipping signs + # print(ele.Name,ele.__dict__[field]) + else: + ele.__dict__[field]=value + + + + def getRegExpElement(self,section,device,field): + val=[] + myRe=re.compile(section+'.*'+device) + for key in self.ElementDB.keys(): + if (myRe.match(key)): + ele=self.ElementDB[key] + if field in ele.__dict__.keys(): + val.append(ele[field]) + return val + + +#------------------------- +# routines from former Energy Manager + + def EnergyAt(self,ele): + if isinstance(ele,str): + name=ele + else: + name=ele.Name + if name in self.EM.Energy.keys(): + return self.EM.Energy[name] + else: + print(name+' is not found in Energy manager of Module Facility!') + return [0.0,0.0] + + def forceEnergyAt(self,elename,Ei): + + # force an update of the energy profile + self.writeFacility(self.EM) + + if Ei<0: + # The initial energy should be positive + return + try: + # Check if elename exists + E0=self.EM.Energy[elename][0] + dE=Ei-E0 +# dE=Ei +# print('energy change',elename,Ei,E0) +# self.EM.Energy[elename][0]=Ei + except: + print('Element %s not found for forced energy update' % elename) + return False + + # First compute the energy with the present gradient and phase + # Note that the energy at the elements upstream of the given initial element is not touched. + + for k in self.EM.Energy.keys(): + if len(self.EM.Energy[k])==2: + self.EM.Energy[k]=[self.EM.Energy[k][0]+dE,self.EM.Energy[k][1]] + elif len(self.EM.Energy[k])==3: + self.EM.Energy[k]=[self.EM.Energy[k][0]+dE,self.EM.Energy[k][1],self.EM.Energy[k][2]] + return True + + def updateModel(self,dict_in, usek1l = True, referenceEnergy = {} ): + # routine to update the online model from a dictionary following the snapshop names as keys + # to do: 1) not enforce updates if the values are not defined in the input dictionary + # 2) Adjust gradient with reference points of beam energy (dipoles) + # 3) update field strength with dipole current and beam energy + # 4) udjust magnet strengths with magnet current, beam energy and calibration + + + # step 1: RF settings + RF = ['SINSB%2.2d' %i for i in range(1,5)] + ['SINXB01'] + ['S10CB%2.2d' % i for i in range(1,10)]+['S20CB%2.2d' % i for i in range(1,5)]+['S30CB%2.2d' % i for i in range(1,14)]+['SATCB01'] + nCav=[1,1,2,2,2] + [4 for i in range(9+4+13+1)] # linac 1+2+3+SATCB01 + #Linac 1 + for i,rf0 in enumerate(RF): + rfvol = 0 + rfphase = 0 + key = '%s-RMSM:SM-SET' % rf0 + if key in dict_in.keys() and int(dict_in[key]) == 9: + key = '%s-RSYS:REQUIRED-OP' % rf0 + if key in dict_in.keys() and int(dict_in[key]) == 1: + key = '%s-RSYS:SET-ACC-VOLT' % rf0 + if key in dict_in.keys(): + rfvol = dict_in[key]/nCav[i] + key = '%s-RSYS:SET-BEAM-PHASE' % rf0 + if key in dict_in.keys(): + rfphase = dict_in[key] + Lens = self.getRegExpElement(rf0,'RACC','Length') + rfvol=rfvol*1e6/Lens[0] + self.setRegExpElement(rf0,'RACC','Gradient',rfvol) + self.setRegExpElement(rf0,'RACC','Phase',rfphase) + # calculate energy profile + self.forceEnergyAt('SINLH01.DBAM010',155.*1e6) + + # step 2: Undulator + UND = ['SARUN%2.2d' % i for i in range(3,16)]+['SATUN%2.2d' %i for i in range(6,14)]+['SATUN%2.2d' %i for i in range(15,23)] + for und0 in UND: + key = '%s-UIND030:K_SET' % und0 + if key in dict_in.keys(): + K = dict_in[key] + self.setRegExpElement(und0,'UIND','K',K) + if 'SATUN' in und0: + key = key.replace('K_SET','POL-SET') + if key in dict_in.keys(): + pol = int(dict_in[key]) # 0 - LH, 1+2 - LV, 3+4 - C + fit = 0.04701805*K*K+ 0.20133565*K+ 0.83336356 # git formula from Apple-X radia model + if pol == 0: + kx = fit + ky = 1 - kx + elif pol == 1 or pol ==2: + ky = fit + kx = 1- ky + else: + kx = 0.5 + ky = 0.5 + self.setRegExpElement(und0,'UIND','kx',kx) + self.setRegExpElement(und0,'UIND','ky',ky) + key = 'SINLH02-UIND230:X-SP' + K = 0 + if key in dict_in.keys(): + x = dict_in[key] + if x > -40.: # undulator is in the beamline + key = 'SINLH02-UIND230:K-SP' + if key in dict_in.keys(): + K = dict_in[key] + self.setRegExpElement('SINLH02','UIND','K',K) + MOD = ['SATDI01','SATUN04'] + idx = [50,30] + for i,mod0 in enumerate(MOD): + for j in range(9): + key1 = '%s-UMOD%3.3d:K-D-SET' % (mod0,j+idx[i]) + key2 = '%s-UMOD%3.3d:K-U-SET' % (mod0,j+idx[i]) + if key1 in dict_in.keys() and key2 in dict_in.keys(): + K = 0.5*(dict_in[key1]+dict_in[key2]) + self.setRegExpElement(mod0,key[8:15],'K',K) + + # step 3 - Magnets + if usek1l: + for key in dict_in.keys(): + if ':K2L-SET' in key: + magnet = key.split(':')[0].replace('-','.') + ele = self.getElement(magnet) + if not ele is None: + ele.k2 = dict_in[key]/ele.Length + if ':K1L-SET' in key: + magnet = key.split(':')[0].replace('-','.') + ele = self.getElement(magnet) + if not ele is None: + ele.k1 = dict_in[key]/ele.Length + + BEND = ['SINLH02-MBND100','SINBC02-MBND100','S10BC02-MBND100','SATMA01-MBND300','SATUN05-MBND100','SATUN14-MBND100'] + for bend0 in BEND: + key = '%s:K0L-SET' % bend0 + if key in dict_in.keys(): + ang = dict_in[key]*180./np.pi + self.setRegExpElement(bend0[0:7],'MBND','angle',ang) + + # step 4 - Dechirper -- to be defined diff --git a/propertytypes.py b/propertytypes.py index f3d4661..93ec309 100644 --- a/propertytypes.py +++ b/propertytypes.py @@ -1,5 +1,3 @@ -from elements import Corrector - # Editor Names Name = ['Name','Prefix','Tag','Index'] Position=['Offset','Relative Offset', 'Length','Reserved Length','Reference Position','Displacement in X','Displacement in Y','Tilt'] diff --git a/protoliste.py b/protoliste.py new file mode 100644 index 0000000..b0505cf --- /dev/null +++ b/protoliste.py @@ -0,0 +1,472 @@ +from onlinemodel.core.apptemplate import ApplicationTemplate +import math +import numpy as np + +class ProtoListe(ApplicationTemplate): +# any calculation for the explicit layout of the machine + + def __init__(self,phase = 0): + #self.SF = Facility(1,phase) # generate instance of facility + self.phase = phase + choice = ['Current', 'Planned', 'Future'] + self.phaseName=choice[phase] + self.info={} + self.order=[] + self.dipole={} + + + + def generateLayout(self,line=None): + print('--------\nGenerating Layout for Phase: %s\n' % self.phaseName) + self.order.clear() + self.info.clear() + self.dipole.clear() + self.setRef(-0.1, -0.1) + line.writeLattice(self) + #self.SF.writeFacility(self) + + def setRef(self,sin=0,zin=0,xin=0,yin=1.2,exin=0,eyin=0,ezin=1): + self.s=sin + self.z=zin + self.x=xin + self.y=yin + # vector of direction + self.ex=exin + self.ey=eyin + self.ez=ezin + + self.vertx=0 + self.verty=0 + self.vertz=0 + self.verts=0 + + # some stuff to make it easier to generate the holy list. + self.p1 = np.array([sin,xin,yin,zin]) # starting point + self.pv = np.array([0,0,0,0]) # center point/vertex point + self.p2 = np.array([0,0,0,0]) # end point + self.ev1 = np.array([1,exin,eyin,ezin]) + self.ev2 = np.array([1,exin,eyin,ezin]) + self.ev = np.array([0,0,0]) + + + def advance(self, L, angle=0, tilt=0): # this is independent on the new coordinate system for the holy list + self.vertx = 0 + self.verty = 0 + self.vertz = 0 + self.verts = self.s + self.s = self.s + L + if angle == 0: + self.x = self.x + L * self.ex + self.y = self.y + L * self.ey + self.z = self.z + L * self.ez + else: + # save initial coordinates + xtmp = self.x + ytmp = self.y + ztmp = self.z + # radius + R = L / angle + # rotation vector - default around y axis + Rproj = self.ey # R dot e = Ry*ey=ey + Rx = -Rproj * self.ex + Ry = 1 - Rproj * self.ey + Rz = -Rproj * self.ez + # Orthonormal vector to direction of propagation and rotation axis + Ox = Ry * self.ez - Rz * self.ey + Oy = Rz * self.ex - Rx * self.ez + Oz = Rx * self.ey - Ry * self.ex + Rnorm = math.sqrt(Rx * Rx + Ry * Ry + Rz * Rz) + Onorm = math.sqrt(Ox * Ox + Oy * Oy + Oz * Oz) + Rx = Rx / Rnorm + Ry = Ry / Rnorm + Rz = Rz / Rnorm + Ox = Ox / Onorm + Oy = Oy / Onorm + Oz = Oz / Onorm + # Rot Vector is the axis of rotation and Ort vector pointing to the origin of the rotation + Ortx = Ox * math.cos(tilt) - Rx * math.sin(tilt) + Orty = Oy * math.cos(tilt) - Ry * math.sin(tilt) + Ortz = Oz * math.cos(tilt) - Rz * math.sin(tilt) + # the origin of the rotation + X0 = self.x + R * Ortx + Y0 = self.y + R * Orty + Z0 = self.z + R * Ortz + # rotation + + self.vertx = X0 - R * (Ortx * math.cos(angle * 0.5) - self.ex * math.sin(angle * 0.5)) / math.cos(angle * 0.5) + self.verty = Y0 - R * (Orty * math.cos(angle * 0.5) - self.ey * math.sin(angle * 0.5)) / math.cos(angle * 0.5) + self.vertz = Z0 - R * (Ortz * math.cos(angle * 0.5) - self.ez * math.sin(angle * 0.5)) / math.cos(angle * 0.5) + self.x = X0 - R * (Ortx * math.cos(angle) - self.ex * math.sin(angle)) + self.y = Y0 - R * (Orty * math.cos(angle) - self.ey * math.sin(angle)) + self.z = Z0 - R * (Ortz * math.cos(angle) - self.ez * math.sin(angle)) + self.ex = (self.ex * math.cos(angle) + Ortx * math.sin(angle)) + self.ey = (self.ey * math.cos(angle) + Orty * math.sin(angle)) + self.ez = (self.ez * math.cos(angle) + Ortz * math.sin(angle)) + # calculate vertex + ds = (self.vertx - xtmp) * (self.vertx - xtmp) + (self.verty - ytmp) * (self.verty - ytmp) + ( + self.vertz - ztmp) * (self.vertz - ztmp) + self.verts = self.verts + math.sqrt(ds) + + + def swap(self): + xtmp=self.x + ytmp=self.y + ztmp=self.z + stmp=self.s + self.x=self.vertx + self.y=self.verty + self.z=self.vertz + self.s=self.verts + self.vertx=xtmp + self.verty=ytmp + self.vertz=ztmp + self.verts=stmp + + + def getDomain(self,name): + domain = 'INJECTOR' + if name[1:3] == '10': + domain = 'LINAC1' + if name[1:3] == '20': + domain = 'LINAC2' + if name[1:3] == '30': + domain = 'LINAC3' + if name[1:3] == 'AR': + domain = 'ARAMIS-U15' + if name[1:3] == 'AT': + domain = 'ATHOS-U40' + if name[1:3] == 'PO': + domain = 'PORTHOS' + return domain + + def write(self,elements): + + name=elements[0] + reference=elements[2] + elelength=elements[1] + group=elements[3] + prefix = name[0:7] + suffix = name[8:12] + if len(name) < 15: + index=0 + else: + index = int(name[12:15]) + baugruppe = elements[5] + roll = elements[6] + psch = elements[7] + + align = 'Entrance' + namePL=name.replace('-','.') + if namePL in self.dipole.keys(): + align = self.dipole[namePL]['Alignment'] + + domain = self.getDomain(name) + + #if (name=='SINEG01.MSOL010'): ############### very first element - very ugly implementation + # self.setRef(-0.085,-0.085) + # self.pv=self.p1 + # self.p1 = self.pv - 0.5 * elelength*self.ev1 + # self.p2 = self.pv + 0.5 * elelength*self.ev1 + + + angh=np.arctan2(self.ev[0],self.ev[2])*180/np.pi + rh = np.sqrt(self.ev[2]*self.ev[2]+self.ev[0]*self.ev[0]) + angv=np.arctan2(self.ev[1],rh)*180/np.pi + + + pos_s=self.s + pos_z=self.z + + if name in self.order: + return + + if name[8:9]=='S': + group='Safety' + + # add some ad hoc shifts of elements + ds = np.array([0,0,0,0]) + + + if name[0:12]=='SINBC02.DCOL' or name[0:12]=='S10BC02.DCOL' or name[0:12]=='S10BC02.DSFH': + reference='Center' + pos_s=pos_s+0.5*elelength + pos_z=pos_z+0.5*elelength + ds = +0.5*elelength * np.array([1, 0, 0, 1]) + + if name[0:12]=='SINBC02.DSCR' or name[0:12]=='S10BC02.DSCR': + elelength=0.28 + pos_s=pos_s-0.09 + pos_z=pos_z-0.09 + ds = -0.09 * np.array([1,0,0,1]) + + + + info={'Domain': domain, 'Prefix': prefix, 'Suffix': suffix, 'Index': index, + 's (m)': pos_s, 'z (m)': pos_z, 'x (m)': self.x, 'y (m)': self.y, 'L (m)': elelength, + 'Reference': reference, 'Group': group, 'Baugruppe': baugruppe, 'Variante': 'A', 'Roll': roll, + 'PS-Ch.': psch,'P1':self.p1+ds,'PV':self.pv+ds,'P2':self.p2+ds,'EV1':self.ev1,'EV2':self.ev2, + 'Yaw': angh, 'Pitch': angv, 'Alignment': align} + + self.info[name] = info + self.order.append(name) + +#---------------------------------------- +# specific elements + + def isType(self,name): + if (name.find('holylist')>-1): + return 1 + else: + return 0 + + def writeLine(self,ele,seq): + if not ele.Name in self.info.keys() and len(ele.Name) > 1: + print('Empty Cell encountered:', ele.Name) + domain = self.getDomain(ele.Name) + align = 'Entrance' + angh = np.arctan2(self.ev[0], self.ev[2]) * 180 / np.pi + rh = np.sqrt(self.ev[2] * self.ev[2] + self.ev[0] * self.ev[0]) + angv = np.arctan2(self.ev[1], rh) * 180 / np.pi + info = {'Domain': domain, 'Prefix': ele.Name, 'Suffix': '', 'Index': '', + 's (m)': self.s, 'z (m)': self.z, 'x (m)': self.x, 'y (m)': self.y, 'L (m)': 0, + 'Reference': 'Start', 'Group': 'Section', 'Baugruppe': 'Marker', 'Variante': 'A', 'Roll': 0, + 'PS-Ch.': 0, 'P1': self.p1, 'PV': self.p1, 'P2': self.p1, 'EV1': self.ev1, + 'EV2': self.ev1, 'Yaw': angh, 'Pitch': angv, 'Alignment': align} + self.info[ele.Name] = info + self.order.append(ele.Name) + return + + def writeDrift(self,s): + self.p1 = self.p1 + s * self.ev1 # no change in vector - should replace advance in the future + self.advance(s,0) + + def writeAlignment(self,ele): + # print info of the element and advances the layout position over the reserved length + if 'dx' in ele.__dict__: + self.x=self.x+ele.dx + self.p1[1] += ele.dx + if 'dy' in ele.__dict__: + self.y=self.y+ele.dy + self.p1[2] += ele.dy + self.advance(ele.getResLength(),0) + + + def writeMarker(self,ele): + + if ele.Group == 'Marker': + if 'MREF' in ele.Name: + self.write((ele.Name, ele.LengthRes, 'Start', ele.Group, ele.Tag, ele.Baugruppe, 0, 0)) + if ele.Group=='Start': + self.write((ele.Name, ele.LengthRes, 'Start', ele.Group, ele.Tag, 'Marker', 0, 0)) + return + + def writeVacuum(self,ele): # all components, which have a well defined vacuum vessel. + self.p2 = self.p1 + ele.LengthRes*self.ev1 + self.pv = self.p1 + self.ev2 = self.ev1 + self.ev = self.ev1[1:4] + name = ele.Prefix+'-'+ele.Name + self.write((name,ele.LengthRes,'Start',ele.Group,ele.Tag,ele.Baugruppe,0,0)) + self.p1 = self.p2 # advance position + self.advance(ele.LengthRes,0) + + def writeRF(self,ele): + self.writeVacuum(ele) + + def writeDiagnostic(self,ele): + self.writeVacuum(ele) + + def writeUndulator(self,ele): + self.writeMagnets(ele) + + def writeCorrector(self,ele): + self.writeMagnets(ele) + + def writeQuadrupole(self,ele): + self.writeMagnets(ele) + + def writeSextupole(self,ele): + self.writeMagnets(ele) + + def writeSolenoid(self,ele): + self.writeMagnets(ele) + + + + def writeMagnets(self,ele): + channels=0 + if 'k1' in ele.__dict__: + channels=channels+1 + if 'k2' in ele.__dict__: + channels=channels+1 + if 'ks' in ele.__dict__: + channels=channels+1 + if 'corx' in ele.__dict__: + channels=channels+1 + if 'cory' in ele.__dict__: + channels=channels+1 + if 'cor' in ele.__dict__: + channels=channels+1 + if ele.Baugruppe=='QFU': + channels=channels-1 + if ele.Baugruppe=='QFUE': + channels=channels-1 + if ele.Baugruppe == 'QFDump': + channels = channels - 1 + if ele.Baugruppe=='UE38-DELAY': + channels=channels+1 + # print info of the element and advances the layout position over the reserved length + self.advance(ele.sRef+ele.Length*0.5,0) # propagation to the center point + self.p1 = self.p1 + ele.sRef*self.ev1 + self.pv = self.p1 + 0.5 * ele.Length*self.ev1 + self.p2 = self.p1 + ele.Length*self.ev1 + self.ev2 = self.ev1 + self.ev = self.ev1[1:4] + name = ele.Prefix+'-'+ele.Name + self.write((name,ele.Length,'Center',ele.Group,ele.Tag,ele.Baugruppe,ele.Tilt*90/math.asin(1),channels)) + self.p1 = self.p2 + (ele.LengthRes-ele.sRef-ele.Length)*self.ev1 + self.advance(ele.LengthRes-ele.sRef-ele.Length*0.5,0) # residual space including half the magnet length + + + + + + def writeBend(self,ele): + # print info of the element and advances the layout position over the reserved length + channels=1 + name = ele.Prefix+'-'+ele.Name + if ele.Baugruppe[0:3]=='AFP': + channels=0 + if 'cor' in ele.__dict__: + channels=channels+1 + self.advance(ele.sRef,0) + self.p1 = self.p1 + ele.sRef*self.ev1 + if (ele.design_angle==0): + self.pv = self.p1 + 0.5*ele.getLength()*self.ev1 + self.p2 = self.p1 + ele.getLength()*self.ev1 + self.ev2 = self.ev1 + self.ev= self.ev1[1:4] + self.advance(0.5*ele.getLength(),0) + if ele.Name not in self.dipole.keys(): + self.dipole[name] = {'Angle': ele.design_angle, 'Alignment': 'Entrance', 'Radius': '', + 'Transverse Wander': 0, 'Shift': 0, 'Length': ele.getLength()} + self.write((name, ele.getLength(), 'Vertex', ele.Group, ele.Tag, ele.Baugruppe, + ele.Tilt * 90 / math.asin(1), channels)) + self.p1 = self.p2+(ele.getResLength()-ele.sRef-ele.getLength())*self.ev1 + self.advance(ele.getResLength()-ele.sRef-0.5*ele.getLength(), 0) + else: + # print('Bending magnet', ele.Name,ele.angle,ele.design_angle,ele.sRef) + xtmp=self.x + ytmp=self.y + ztmp=self.z + stmp=self.s + extmp=self.ex + eytmp=self.ey + eztmp=self.ez + angletmp=ele.angle + ele.angle=ele.design_angle + angrad = np.pi* ele.angle / 180. + + # ignoring tilt: + R = ele.getLength()/angrad + M = np.zeros((4,4)) + + M[0,0]=1 + M[2,2]= 1 + M[1, 1] = np.cos(angrad) + M[3, 3] = np.cos(angrad) + M[1, 3] = np.sin(angrad) + M[3, 1] = -np.sin(angrad) + self.ev2 = np.matmul(M,self.ev1) + if not ele.Tilt == 0: + vec1 = np.array(self.ev1[1:]) + vec2 = np.array(self.ev2[1:]) + proj = np.sum(vec1*vec2) + vecper = vec2 - proj*vec1 + vecpar = vec2 - vecper + veccro = np.cross(vec1,vecper) + self.ev2[1:] = vecpar+np.cos(ele.Tilt)*vecpar-np.sin(ele.Tilt)*veccro + + Ld = R * np.sin(angrad*0.5) + T = Ld/np.cos(angrad*0.5) + self.pv = self.p1 + T * self.ev1 + self.p2 = self.pv + T * self.ev2 + self.ev = (1-ele.e1) * self.ev1[1:4] + (1-ele.e2) * self.ev2[1:4] + + # calculate shifts of dipoles for good field region + R = ele.getLength() / angrad; + dR = R * (1 - math.cos(angrad)) # beam wander + shift = 0.5 + if (ele.e1 == 0.5): + dR = R * (1 - math.cos(angrad * 0.5)) + shift = -0.5 + vec1 = np.array(self.ev1[1:]) + vec2 = np.array(self.ev2[1:]) + vnorm = np.cross(vec1, vec2) + vrot = (1 - ele.e1) * vec1 + (1 - ele.e2) * vec2 + vshift = np.cross(vnorm, vrot) + vlen = np.sqrt(np.sum(vshift*vshift)) + # ############################### + # disable shift in holy list!!!! + # + vlen = 0 + + if vlen > 0: + vshift = vshift /vlen + eshift = np.sign(angrad)*np.array([0,vshift[0],vshift[1],vshift[2]])*shift*dR + self.p1 += eshift + self.pv += eshift + self.p2 += eshift + + # advance according to design angle + self.advance(ele.getLength(), angrad, ele.Tilt) + self.swap() + + if ele.Name not in self.dipole.keys(): + if ele.e2 == 1: + align = 'Entrance' + elif ele.e1 ==1: + align = 'Exit' + else: + align = 'Middle' + self.dipole[name] = {'Angle': ele.design_angle, 'Alignment': align, 'Radius': np.abs(R), + 'Transverse Wander': np.abs(dR), 'Shift': shift*dR, 'Length': ele.getLength()} + + self.write((name, ele.Length, 'Vertex', ele.Group, ele.Tag, ele.Baugruppe, ele.Tilt * 90 / math.asin(1), channels)) + self.swap() + + if vlen > 0: + self.p1 -= eshift + self.pv -= eshift + self.p2 -= eshift + + # restore position before dipole and then advance according to given angle (this matters for branching point dipoles + self.x=xtmp + self.y=ytmp + self.z=ztmp + self.s=stmp + self.ex=extmp + self.ey=eytmp + self.ez=eztmp + + + + + + ele.angle=angletmp + if ele.angle == 0: + self.p1 = self.p1 + ele.getLength() *self.ev1 + else: + ps = self.p1[0]+ele.getLength() + self.p1 = self.p2 + self.p1[0]=ps + self.ev1 = self.ev2 + angrad=np.pi*ele.angle/180. + Lpath=ele.getLength() + self.advance(Lpath,angrad,ele.Tilt) + self.advance(ele.getResLength()-ele.sRef-Lpath,0) + self.p1 = self.p1 + (ele.getResLength()-ele.sRef-Lpath)*self.ev1 + + + +