From 4ed38feaacf54de55ec7494b76c1edd1030504f4 Mon Sep 17 00:00:00 2001 From: Shishlo Date: Fri, 26 Sep 2025 13:17:58 -0400 Subject: [PATCH 1/2] The method of setUsePhaseAtEntrance(usePhaseAtEntrance = True/False) was added to RF_Cavity class. It defines the place from where the cavity phase is defined. True means at the entrance of the 1st RF gap of the cavity, and False means the center of 1st gap. This parameter is used only for Axis Field RF gaps. --- .../py_linac/lattice/LinacAccLatticeLib.py | 28 +++++ .../lattice/LinacFieldOverlappingNodes.py | 109 ++++++++-------- py/orbit/py_linac/lattice/LinacRfGapNodes.py | 117 +++++++++--------- 3 files changed, 134 insertions(+), 120 deletions(-) diff --git a/py/orbit/py_linac/lattice/LinacAccLatticeLib.py b/py/orbit/py_linac/lattice/LinacAccLatticeLib.py index bdb2ebc0..b76455f8 100644 --- a/py/orbit/py_linac/lattice/LinacAccLatticeLib.py +++ b/py/orbit/py_linac/lattice/LinacAccLatticeLib.py @@ -375,6 +375,14 @@ def __init__(self, name="none"): self.addParam("designArrivalTime", 0.0) self.addParam("isDesignSetUp", False) self.addParam("pos", 0.0) + #---- This parameter is only for RF gap types with + #---- the continuous RF fields on the axis of the cavity. + #---- If usePhaseAtEntrance = False we will use the phase at the center + #---- of RF gap as for a standard cavity representations as a set + #---- of zero-length RF gaps. + #---- If usePhaseAtEntrance = False we use the phase at the entrance + #---- of the cavity. + self.usePhaseAtEntrance = False def setDesignSetUp(self, designOnOf): """Sets the design set up information (True,False).""" @@ -482,6 +490,26 @@ def getAvgGapPhase(self): def getAvgGapPhaseDeg(self): """Returns average phase in degrees for all RF gaps in the cavity""" return self.getAvgGapPhase() * 180.0 / math.pi + + def setUsePhaseAtEntrance(self,usePhaseAtEntrance): + """ + This parameter is only for RF gap types with the continuous RF fields + on the axis of the cavity. + If usePhaseAtEntrance = False we will use the phase at the center + of 1 st RF gap as for a standard cavity representations as a set + of zero-length RF gaps. + If usePhaseAtEntrance = False we use the phase at the entrance + of the cavity. + By default it is False. Switching to True will change cavity phase and + all longitudinal dynamics calculations. You should calculate and set + of the cavity phase before the bunch tracking. And, of course, you + should start with the trackDesignBunch(...) lattice method. + """ + self.usePhaseAtEntrance = usePhaseAtEntrance + + def getUsePhaseAtEntrance(self): + """ Returns the usePhaseAtEntrance parameter of the cavity """ + return self.usePhaseAtEntrance def removeAllGapNodes(self): """Remove all rf gaps from this cavity.""" diff --git a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py index c1e8f3cc..eae8cc49 100755 --- a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py +++ b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py @@ -21,6 +21,8 @@ # from linac import the RF gap classes from orbit.core.linac import RfGapThreePointTTF, RfGapThreePointTTF_slow +# import from orbit c++ utilities +from orbit.core.orbit_utils import Function class AxisField_and_Quad_RF_Gap(AbstractRF_Gap): """ @@ -61,8 +63,8 @@ def __init__(self, axis_field_rf_gap): self.z_step = 0.01 self.z_min = 0.0 self.z_max = 0.0 - # ---- gap_phase_vs_z_arr keeps [pos,phase] pairs after the tracking - self.gap_phase_vs_z_arr = [] + # ---- gap_pos_phase_func keeps phase=function(pos) after the tracking + self.gap_pos_phase_func = Function() # ---- The position of the particle during the run. # ---- It is used for the path length accounting. self.part_pos = 0.0 @@ -300,6 +302,9 @@ def track(self, paramsDict): designArrivalTime = rfCavity.getDesignArrivalTime() phase_shift = rfCavity.getPhase() - rfCavity.getDesignPhase() phase = rfCavity.getFirstGapEtnrancePhase() + phase_shift + usePhaseAtEntrance = rfCavity.getUsePhaseAtEntrance() + if(usePhaseAtEntrance): + phase = rfCavity.getPhase() # ---------------------------------------- phase = math.fmod( frequency * (arrival_time - designArrivalTime) * 2.0 * math.pi + phase, @@ -307,9 +312,8 @@ def track(self, paramsDict): ) if index == 0: self.part_pos = self.z_min - self.gap_phase_vs_z_arr = [ - [self.part_pos, phase], - ] + self.gap_pos_phase_func.clean() + self.gap_pos_phase_func.add(self.part_pos, phase) zm = self.part_pos z0 = zm + part_length / 2 zp = z0 + part_length / 2 @@ -338,7 +342,7 @@ def track(self, paramsDict): # call rf gap model to track the bunch time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -377,7 +381,7 @@ def track(self, paramsDict): self.part_pos += part_length / 2 time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -388,29 +392,35 @@ def track(self, paramsDict): # ---- this part is the debugging ---STOP--- # ---- Calculate the phase at the center if index == (nParts - 1): - pos_old = self.gap_phase_vs_z_arr[0][0] - phase_gap = self.gap_phase_vs_z_arr[0][1] - ind_min = -1 - for ind in range(1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - if math.fabs(pos) >= math.fabs(pos_old): - ind_min = ind - 1 - phase_gap = self.gap_phase_vs_z_arr[ind_min][1] - phase_gap = phaseNearTargetPhase(phase_gap, 0.0) - self.gap_phase_vs_z_arr[ind_min][1] = phase_gap - break - pos_old = pos - self.setGapPhase(phase_gap) - # ---- wrap all gap part's phases around the central one - if ind_min > 0: - for ind in range(ind_min - 1, -1, -1): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind + 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) - for ind in range(ind_min + 1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind - 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) + self._phase_gap_func_analysis() + + def _phase_gap_func_analysis(self): + """ + Performs analysis of the RF gap function phase vs. postion + and calculates the accelerating phase at the center of the gap. + """ + n_pos_points = self.gap_pos_phase_func.getSize() + phase_gap = 0. + if(n_pos_points != 0): + phase_gap_old = self.gap_pos_phase_func.y(0) + for pos_ind in range(n_pos_points): + phase_gap = self.gap_pos_phase_func.y(pos_ind) + phase_gap = phaseNearTargetPhase(phase_gap,phase_gap_old) + self.gap_pos_phase_func.updatePoint(pos_ind,phase_gap) + phase_gap_old = phase_gap + phase_gap = self.gap_pos_phase_func.getY(0.) + phase_gap_center = phaseNearTargetPhase(phase_gap,0.) + phase_gap_shift = phase_gap - phase_gap_center + for pos_ind in range(n_pos_points): + phase_gap = self.gap_pos_phase_func.y(pos_ind) + self.gap_pos_phase_func.updatePoint(pos_ind,phase_gap - phase_gap_shift) + self.setGapPhase(phase_gap_center) + + def getPhaseVsPositionFuncion(self): + """ + Retuns the RF gap function phase vs. postion. + """ + return self.gap_pos_phase_func def trackDesign(self, paramsDict): """ @@ -433,10 +443,14 @@ def trackDesign(self, paramsDict): arrival_time = syncPart.time() frequency = rfCavity.getFrequency() phase = rfCavity.getFirstGapEtnrancePhase() + usePhaseAtEntrance = rfCavity.getUsePhaseAtEntrance() + if(usePhaseAtEntrance): + phase = rfCavity.getPhase() # ---- calculate the entance phase if self.isFirstRFGap() and index == 0: rfCavity.setDesignArrivalTime(arrival_time) - phase = self.axis_field_rf_gap.calculate_first_part_phase(bunch) + if(not usePhaseAtEntrance): + phase = self.axis_field_rf_gap.calculate_first_part_phase(bunch) rfCavity.setFirstGapEtnrancePhase(phase) rfCavity.setFirstGapEtnranceDesignPhase(phase) rfCavity.setDesignSetUp(True) @@ -453,9 +467,8 @@ def trackDesign(self, paramsDict): # print "debug design name=",self.getName()," arr_time=",arrival_time," phase=",phase*180./math.pi," E0TL=",E0TL*1.0e+3," freq=",frequency if index == 0: self.part_pos = self.z_min - self.gap_phase_vs_z_arr = [ - [self.part_pos, phase], - ] + self.gap_pos_phase_func.clean() + self.gap_pos_phase_func.add(self.part_pos, phase) zm = self.part_pos z0 = zm + part_length / 2 zp = z0 + part_length / 2 @@ -468,7 +481,7 @@ def trackDesign(self, paramsDict): # call rf gap model to track the bunch time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -491,7 +504,7 @@ def trackDesign(self, paramsDict): self.part_pos += part_length / 2 time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -502,29 +515,7 @@ def trackDesign(self, paramsDict): # ---- this part is the debugging ---STOP--- # ---- Calculate the phase at the center if index == (nParts - 1): - pos_old = self.gap_phase_vs_z_arr[0][0] - phase_gap = self.gap_phase_vs_z_arr[0][1] - ind_min = -1 - for ind in range(1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - if math.fabs(pos) >= math.fabs(pos_old): - ind_min = ind - 1 - phase_gap = self.gap_phase_vs_z_arr[ind_min][1] - phase_gap = phaseNearTargetPhase(phase_gap, 0.0) - self.gap_phase_vs_z_arr[ind_min][1] = phase_gap - break - pos_old = pos - self.setGapPhase(phase_gap) - # ---- wrap all gap part's phases around the central one - if ind_min > 0: - for ind in range(ind_min - 1, -1, -1): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind + 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) - for ind in range(ind_min + 1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind - 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) + self._phase_gap_func_analysis() class OverlappingQuadsNode(BaseLinacNode): diff --git a/py/orbit/py_linac/lattice/LinacRfGapNodes.py b/py/orbit/py_linac/lattice/LinacRfGapNodes.py index d5b062e5..2d6e8f1b 100644 --- a/py/orbit/py_linac/lattice/LinacRfGapNodes.py +++ b/py/orbit/py_linac/lattice/LinacRfGapNodes.py @@ -268,6 +268,8 @@ def trackDesign(self, paramsDict): if self.__isFirstGap: rfCavity.setDesignArrivalTime(arrival_time) rfCavity.setDesignSetUp(True) + rfCavity.setFirstGapEtnrancePhase(rfCavity.getPhase()) + rfCavity.setFirstGapEtnranceDesignPhase(rfCavity.getPhase()) rfCavity._setDesignPhase(rfCavity.getPhase()) rfCavity._setDesignAmp(rfCavity.getAmp()) else: @@ -469,8 +471,8 @@ def __init__(self, baserf_gap): self.z_max = 0.0 self.z_tolerance = 0.000001 # in meters self.phase_tolerance = 0.001 # in degrees - # ---- gap_phase_vs_z_arr keeps [pos,phase] pairs after the tracking - self.gap_phase_vs_z_arr = [] + # ---- gap_pos_phase_func keeps phase=function(pos) after the tracking + self.gap_pos_phase_func = Function() # ---- The position of the particle during the run. # ---- It is used for the path length accounting. self.part_pos = 0.0 @@ -508,11 +510,15 @@ def getAxisFieldFunction(self): """ return self.axis_field_func - def setAxisFieldFunction(self, axis_field_func): + def setAxisFieldFunction(self, axis_field_func, z_step=0.01): """ It sets the axis field function. """ self.axis_field_func = axis_field_func + z_min = self.axis_field_func.getMinX() + z_max = self.axis_field_func.getMaxX() + self.z_step = z_step + self.setZ_Min_Max(z_min, z_max) def getZ_Step(self): """ @@ -611,13 +617,15 @@ def track(self, paramsDict): designArrivalTime = rfCavity.getDesignArrivalTime() phase_shift = rfCavity.getPhase() - rfCavity.getDesignPhase() phase = rfCavity.getFirstGapEtnrancePhase() + phase_shift + usePhaseAtEntrance = rfCavity.getUsePhaseAtEntrance() + if(usePhaseAtEntrance): + phase = rfCavity.getPhase() # ---------------------------------------- phase = math.fmod(frequency * (arrival_time - designArrivalTime) * 2.0 * math.pi + phase, 2.0 * math.pi) if index == 0: self.part_pos = self.z_min - self.gap_phase_vs_z_arr = [ - [self.part_pos, phase], - ] + self.gap_pos_phase_func.clean() + self.gap_pos_phase_func.add(self.part_pos, phase) zm = self.part_pos z0 = zm + part_length / 2 zp = z0 + part_length / 2 @@ -630,7 +638,7 @@ def track(self, paramsDict): # call rf gap model to track the bunch time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -645,7 +653,7 @@ def track(self, paramsDict): self.part_pos += part_length / 2 time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -656,29 +664,35 @@ def track(self, paramsDict): # ---- this part is the debugging ---STOP--- # ---- Calculate the phase at the center if index == (nParts - 1): - pos_old = self.gap_phase_vs_z_arr[0][0] - phase_gap = self.gap_phase_vs_z_arr[0][1] - ind_min = -1 - for ind in range(1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - if math.fabs(pos) >= math.fabs(pos_old): - ind_min = ind - 1 - phase_gap = self.gap_phase_vs_z_arr[ind_min][1] - phase_gap = phaseNearTargetPhase(phase_gap, 0.0) - self.gap_phase_vs_z_arr[ind_min][1] = phase_gap - break - pos_old = pos - self.setGapPhase(phase_gap) - # ---- wrap all gap part's phases around the central one - if ind_min > 0: - for ind in range(ind_min - 1, -1, -1): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind + 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) - for ind in range(ind_min + 1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind - 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) + self._phase_gap_func_analysis() + + def _phase_gap_func_analysis(self): + """ + Performs analysis of the RF gap function phase vs. postion + and calculates the accelerating phase at the center of the gap. + """ + n_pos_points = self.gap_pos_phase_func.getSize() + phase_gap = 0. + if(n_pos_points != 0): + phase_gap_old = self.gap_pos_phase_func.y(0) + for pos_ind in range(n_pos_points): + phase_gap = self.gap_pos_phase_func.y(pos_ind) + phase_gap = phaseNearTargetPhase(phase_gap,phase_gap_old) + self.gap_pos_phase_func.updatePoint(pos_ind,phase_gap) + phase_gap_old = phase_gap + phase_gap = self.gap_pos_phase_func.getY(0.) + phase_gap_center = phaseNearTargetPhase(phase_gap,0.) + phase_gap_shift = phase_gap - phase_gap_center + for pos_ind in range(n_pos_points): + phase_gap = self.gap_pos_phase_func.y(pos_ind) + self.gap_pos_phase_func.updatePoint(pos_ind,phase_gap - phase_gap_shift) + self.setGapPhase(phase_gap_center) + + def getPhaseVsPositionFuncion(self): + """ + Retuns the RF gap function phase vs. postion. + """ + return self.gap_pos_phase_func def trackDesign(self, paramsDict): """ @@ -701,10 +715,14 @@ def trackDesign(self, paramsDict): arrival_time = syncPart.time() frequency = rfCavity.getFrequency() phase = rfCavity.getFirstGapEtnrancePhase() + usePhaseAtEntrance = rfCavity.getUsePhaseAtEntrance() + if(usePhaseAtEntrance): + phase = rfCavity.getPhase() # ---- calculate the entance phase if self.isFirstRFGap() and index == 0: rfCavity.setDesignArrivalTime(arrival_time) - phase = self.__calculate_first_part_phase(bunch) + if(not usePhaseAtEntrance): + phase = self.__calculate_first_part_phase(bunch) rfCavity.setFirstGapEtnrancePhase(phase) rfCavity.setFirstGapEtnranceDesignPhase(phase) rfCavity.setDesignSetUp(True) @@ -717,9 +735,8 @@ def trackDesign(self, paramsDict): phase = math.fmod(frequency * (arrival_time - first_gap_arr_time) * 2.0 * math.pi + phase, 2.0 * math.pi) if index == 0: self.part_pos = self.z_min - self.gap_phase_vs_z_arr = [ - [self.part_pos, phase], - ] + self.gap_pos_phase_func.clean() + self.gap_pos_phase_func.add(self.part_pos, phase) # print "debug design name=",self.getName()," index=",index," pos=",self.part_pos," arr_time=",arrival_time," phase=",phase*180./math.pi," freq=",frequency zm = self.part_pos z0 = zm + part_length / 2 @@ -733,7 +750,7 @@ def trackDesign(self, paramsDict): # call rf gap model to track the bunch time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -748,7 +765,7 @@ def trackDesign(self, paramsDict): self.part_pos += part_length / 2 time_middle_gap = syncPart.time() - arrival_time delta_phase = math.fmod(2 * math.pi * time_middle_gap * frequency, 2.0 * math.pi) - self.gap_phase_vs_z_arr.append([self.part_pos, phase + delta_phase]) + self.gap_pos_phase_func.add(self.part_pos, phase + delta_phase) # ---- this part is the debugging ---START--- # eKin_out = syncPart.kinEnergy() # s = "debug pos[mm]= %7.2f "%(self.part_pos*1000.) @@ -759,29 +776,7 @@ def trackDesign(self, paramsDict): # ---- this part is the debugging ---STOP--- # ---- Calculate the phase at the center if index == (nParts - 1): - pos_old = self.gap_phase_vs_z_arr[0][0] - phase_gap = self.gap_phase_vs_z_arr[0][1] - ind_min = -1 - for ind in range(1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - if math.fabs(pos) >= math.fabs(pos_old): - ind_min = ind - 1 - phase_gap = self.gap_phase_vs_z_arr[ind_min][1] - phase_gap = phaseNearTargetPhase(phase_gap, 0.0) - self.gap_phase_vs_z_arr[ind_min][1] = phase_gap - break - pos_old = pos - self.setGapPhase(phase_gap) - # ---- wrap all gap part's phases around the central one - if ind_min > 0: - for ind in range(ind_min - 1, -1, -1): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind + 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) - for ind in range(ind_min + 1, len(self.gap_phase_vs_z_arr)): - [pos, phase_gap] = self.gap_phase_vs_z_arr[ind] - [pos, phase_gap1] = self.gap_phase_vs_z_arr[ind - 1] - self.gap_phase_vs_z_arr[ind][1] = phaseNearTargetPhase(phase_gap, phase_gap1) + self._phase_gap_func_analysis() def calculate_first_part_phase(self, bunch_in): """ @@ -862,5 +857,5 @@ def __calculate_first_part_phase(self, bunch_in): phase_start -= 0.8 * (phase_cavity_new - phase_cavity) # ---- undo the last change in the while loop phase_start += 0.8 * (phase_cavity_new - phase_cavity) - # print "debug phase_start=",phase_start*180./math.pi + #print ("debug phase_start=",phase_start*180./math.pi) return phase_start From bd8dc6e90daa7c68c62366561101fa4a5e0ba41d Mon Sep 17 00:00:00 2001 From: Shishlo Date: Fri, 26 Sep 2025 13:29:35 -0400 Subject: [PATCH 2/2] The test of AxisFieldRF_Gap model and benchmark with BaseRfGap RF gap model. --- .../rfgap_model_longitudinal_test.py | 329 ++++++++++++++++++ 1 file changed, 329 insertions(+) create mode 100755 examples/SNS_Linac/rf_quad_overlapping_fields/rfgap_model_longitudinal_test.py diff --git a/examples/SNS_Linac/rf_quad_overlapping_fields/rfgap_model_longitudinal_test.py b/examples/SNS_Linac/rf_quad_overlapping_fields/rfgap_model_longitudinal_test.py new file mode 100755 index 00000000..6660ead1 --- /dev/null +++ b/examples/SNS_Linac/rf_quad_overlapping_fields/rfgap_model_longitudinal_test.py @@ -0,0 +1,329 @@ +#-------------------------------------------------------------------- +# This script tests the functionality of AxisFieldRF_Gap class. +# It also performs benchmark this calss with the zero-length BaseRfGap() model. +# There are cases: +# ============================================================= +# Case 1: Cavity phase scan assuming the phase at z=0 position. +# Case 2: Cavity phase scan assuming the phase at the entrance +# of the cavity +# Case 3: Cavity phase scan with non-empty bunch assuming the phase +# at the entrance of the cavity +# Case 4: Thin (zero length) RF Gap and Axis Field Gap comparison +#--------------------------------------------------------------------- + +import math +import random +import time +import sys + +# from linac import the C++ RF gap classes +from orbit.core.linac import BaseRfGap, MatrixRfGap, RfGapTTF + +from orbit.core.bunch import Bunch +from orbit.py_linac.lattice import BaseRF_Gap, AxisFieldRF_Gap, RF_Cavity +from orbit.py_linac.lattice import Drift + +from orbit.core.orbit_utils import Function, SplineCH, GaussLegendreIntegrator, Polynomial +from orbit.utils.fitting import PolynomialFit + +from orbit.utils import phaseNearTargetPhase +from orbit.utils import phaseNearTargetPhaseDeg + +#-------------------------------------------- +# Auxiliary Functions +#-------------------------------------------- + +def getRectangularRfAxisField(z_min,z_max): + """ + The RF axial field step function normalized to the integral value of 1. + The center position will be assumed at z=0. + """ + value = 1./(z_max - z_min) + npoints = 10 + z_step = (z_max - z_min)/(npoints - 1) + func = Function() + for ind in range(npoints): + z = z_min + z_step*ind + func.add(z,value) + return func + +def transitTimeFactor(bunch,cav_frequency,length): + """ + Returns the TTF parameter for the thin RF Gap + with a rectangular axis field. + beta - v/c -relativistic factor + w = cav_frequency [Hz] + L2 = length/2 [m] + lambda = beta*c/w + TTF = lambda/L2*sin((L2/lambda) + """ + L2 = length/2 + beta = bunch.getSyncParticle().beta() + w = 2*math.pi*cav_frequency + lmbd = beta*2.997924e+8/w + TTF = (lmbd/L2)*math.sin(L2/lmbd) + return TTF +#-------------------------------------------- +# Script start +#-------------------------------------------- + +#---- The RF axial field step function normalized to the integral value of 1. +z_min = -0.1 +z_max = 0.1 +axis_field_func = getRectangularRfAxisField(z_min,z_max) + +#---- Integral E*L= 20 MeV , here E0L parameter is in GeV +E0L = 0.020 +#---- The RF gap object representing one RF gap +accNode = BaseRF_Gap("BaseRfGap") +three_point_gap = AxisFieldRF_Gap(accNode) +three_point_gap.setAsFirstRFGap(True) +three_point_gap.setAxisFieldFunction(axis_field_func, z_step = 0.01) +three_point_gap.setParam("E0L", E0L) +#---- mode = 0 - phase as it is, mode = 1 shift phase by +PI +three_point_gap.addParam("mode", 1) + +#---- The RF cavity. Here it has only one RF gap +cav_amp = 1.0 +cav_frequency = 704.42e6 +cav = RF_Cavity("AxisFieldCavity") +cav.setAmp(cav_amp) +cav.setFrequency(cav_frequency) +cav.setPosition((z_min+z_max)/2.) +cav.addRF_GapNode(three_point_gap) + +#---- At the entrance of ESS spoke. +#---- Kinetic energy and mass in GeV +mass = 0.93827208943 +Ekin_init = 0.6201 +bunch_init = Bunch() +bunch_init.mass(mass) +bunch_init.charge(1.0) +bunch_init.getSyncParticle().kinEnergy(Ekin_init) + +#---- Parameters of the cavity phase scans. +#---- Phases are in radians. +Nph = 72 +phase_step = 2*math.pi/Nph + +#-------------------------------------------------------------------- +# Case 1: Cavity phase scan assuming the phase at z=0 position. +# We have several types of parameters: +# cav_phase - the input parameter. It is a phase at +# the center of the gap +# cav_1st_gap_entr_phase0 - the phase at the entrance of the 1-st gap +# (we have only one gap) for the cav_phase = 0. +# cav_1st_gap_entr_phase - the phase at the entrance of the 1-st gap +# for the defined cav_phase +# gap_phase - the calculated phase in the center of the RF gap after fitting +# process to get the cav_phase (default accuracy is 0.001 deg) +# delta_phase = cav_1st_gap_entr_phase - cav_phase - cav_1st_gap_entr_phase0 +# this phase will be 0 at cav_phase = 0, and it characterizes +# how wrong we are from the realistic phase scan when we +# we define cavity phase at the entrance of the cavity +# delta_eKin_out - the bunch energy gain after acceleration +#-------------------------------------------------------------------- +three_point_gap.getRF_Cavity().setPhase(0.) +bunch = Bunch() +bunch_init.copyEmptyBunchTo(bunch) +three_point_gap.trackDesignBunch(bunch) +cav_1st_gap_entr_phase0 = three_point_gap.getRF_Cavity().getFirstGapEtnrancePhase() + +print ("=============================================================") +print ("Case 1: Cavity phase scan assuming the phase at z=0 position.") +print ("=============================================================") + +st = "CavPhase[deg] CavPhaseErr[deg] DeltaPhase[deg] 1stGapPhase[deg] DeltaE[MeV] " +print (st) + +for ind in range(Nph+1): + cav_phase = ind*phase_step + three_point_gap.getRF_Cavity().setPhase(cav_phase) + + bunch = Bunch() + bunch_init.copyEmptyBunchTo(bunch) + three_point_gap.trackDesignBunch(bunch) + + cav_1st_gap_entr_phase = phaseNearTargetPhase(three_point_gap.getRF_Cavity().getFirstGapEtnrancePhase(),0.) + delta_phase = phaseNearTargetPhase(cav_1st_gap_entr_phase - cav_phase - cav_1st_gap_entr_phase0,0.) + gap_phase = phaseNearTargetPhase(three_point_gap.getGapPhase(),cav_phase) + delta_eKin_out = bunch.getSyncParticle().kinEnergy() - Ekin_init + st = " %+8.2f "%(cav_phase*180./math.pi) + st += " %+9.6f "%((gap_phase-cav_phase)*180./math.pi) + st += " %+8.4f "%(delta_phase*180./math.pi) + st += " %+8.2f "%(cav_1st_gap_entr_phase*180./math.pi) + st += " %+10.6f "%(delta_eKin_out*1000.) + print (st) + + +#-------------------------------------------------------------------- +# Case 2: Cavity phase scan assuming the phase at the entrance +# of the cavity +# We have several types of parameters: +# cav_phase - the input parameter. It is a phase at +# the center of the gap +# cav_1st_gap_entr_phase0 - the phase at the entrance of the 1-st gap +# (we have only one gap) for the cav_phase = 0. +# cav_1st_gap_entr_phase - the phase at the entrance of the 1-st gap +# for the defined cav_phase +# gap_phase - the calculated phase in the center of the RF gap after fitting +# process to get the cav_phase (default accuracy is 0.001 deg) +# delta_phase = cav_1st_gap_entr_phase - cav_phase - cav_1st_gap_entr_phase0 +# this phase will be 0 at cav_phase = 0, and it characterizes +# how wrong we are from the realistic phase scan when we +# we define cavity phase at the entrance of the cavity +# delta_eKin_out - the bunch energy gain after acceleration +#-------------------------------------------------------------------- + +#---- Set the cavity property +three_point_gap.getRF_Cavity().setUsePhaseAtEntrance(True) + +three_point_gap.getRF_Cavity().setPhase(0.) +bunch = Bunch() +bunch_init.copyEmptyBunchTo(bunch) +three_point_gap.trackDesignBunch(bunch) +cav_1st_gap_entr_phase0 = three_point_gap.getRF_Cavity().getFirstGapEtnrancePhase() + +print ("=============================================================") +print ("Case 2: Cavity phase scan assuming the phase at the entrance.") +print ("=============================================================") + +st = "CavPhase[deg] CavPhaseErr[deg] DeltaPhase[deg] 1stGapPhase[deg] DeltaE[MeV] " +print (st) + +for ind in range(Nph+1): + cav_phase = ind*phase_step + three_point_gap.getRF_Cavity().setPhase(cav_phase) + + bunch = Bunch() + bunch_init.copyEmptyBunchTo(bunch) + three_point_gap.trackDesignBunch(bunch) + + cav_1st_gap_entr_phase = phaseNearTargetPhase(three_point_gap.getRF_Cavity().getFirstGapEtnrancePhase(),0.) + delta_phase = phaseNearTargetPhase(cav_1st_gap_entr_phase - cav_phase - cav_1st_gap_entr_phase0,0.) + gap_phase = phaseNearTargetPhase(three_point_gap.getGapPhase(),cav_phase) + delta_eKin_out = bunch.getSyncParticle().kinEnergy() - Ekin_init + st = " %+8.2f "%(cav_phase*180./math.pi) + st += " %+8.4f "%((gap_phase-cav_phase)*180./math.pi) + st += " %+8.2f "%(delta_phase*180./math.pi) + st += " %+8.2f "%(cav_1st_gap_entr_phase*180./math.pi) + st += " %+10.6f "%(delta_eKin_out*1000.) + print (st) + +#-------------------------------------------------------------------- +# Case 3: Cavity phase scan with non-empty bunch assuming the phase +# at the entrance of the cavity. +# We have several types of parameters: +# cav_phase - the input parameter. It is a phase at +# the center of the gap +# delta_eKin_out - the bunch energy gain after acceleration +# (x,xp,y,yp,z,dE) - 6D coordinates of the particle +#-------------------------------------------------------------------- + +#---- Set the cavity property +three_point_gap.getRF_Cavity().setUsePhaseAtEntrance(True) + +three_point_gap.getRF_Cavity().setPhase(0.) +bunch = Bunch() +bunch_init.copyEmptyBunchTo(bunch) +three_point_gap.trackDesignBunch(bunch) + +print ("=============================================================") +print ("Case 3: Cavity phase scan with one particle in the bunch.") +print ("=============================================================") + +st = "CavPhase[deg] DeltaE[MeV] x[mm] xp[mrad] y[mm] yp[mrad] z[mm] dE[MeV] " +print (st) + +for ind in range(Nph+1): + cav_phase = ind*phase_step + three_point_gap.getRF_Cavity().setPhase(cav_phase) + + bunch = Bunch() + bunch_init.copyBunchTo(bunch) + bunch.addParticle(0.01,0.,0.01,0.,0.,0.) + + three_point_gap.trackBunch(bunch) + + delta_eKin_out = bunch.getSyncParticle().kinEnergy() - Ekin_init + st = " %+8.2f %+8.4f "%(cav_phase*180./math.pi,delta_eKin_out*1000.) + st += " %8.3f %+8.3f "%(bunch.x(0)*1000.,bunch.xp(0)*1000.) + st += " %8.3f %+8.3f "%(bunch.y(0)*1000.,bunch.yp(0)*1000.) + st += " %8.3f %+8.3f "%(bunch.z(0)*1000.,bunch.dE(0)*1000.) + print (st) + +print ("====================================================================") +print ("Case 4: Thin (zero length) RF Gap and Axis Field Gap comparison") +print ("====================================================================") +TTF = transitTimeFactor(bunch_init,cav_frequency,z_max-z_min) +print ("Zero-length RF Gap: Transit Time Factor TTF=",TTF) + +base_rf_gap = BaseRF_Gap("BaseRfGap") +base_rf_gap.setParam("E0TL", E0L*TTF) +#---- setting RF gap tracker +base_rf_gap.setCppGapModel(BaseRfGap()) + +#---- The RF cavity for zero-length RF gap. Here it has only one RF gap +cav_amp = 1.0 +cav_frequency = 704.42e6 +cav_base = RF_Cavity("BaseGapCavity") +cav_base.setAmp(cav_amp) +cav_base.setFrequency(cav_frequency) +cav_base.setPosition((z_min+z_max)/2.) +cav_base.addRF_GapNode(base_rf_gap) + +#---- Tracking design for both types of RF gaps +base_rf_gap.getRF_Cavity().setPhase(0.) +bunch = Bunch() +bunch_init.copyEmptyBunchTo(bunch) +base_rf_gap.trackDesignBunch(bunch) + +three_point_gap.getRF_Cavity().setPhase(0.) +bunch = Bunch() +bunch_init.copyEmptyBunchTo(bunch) +three_point_gap.trackDesignBunch(bunch) + +#---- drifts before and after zero-length RF gap +drift_before = Drift("drift_before") +drift_before.setLength((z_max-z_min)/2) +drift_after = Drift("drift_after") +drift_after.setLength((z_max-z_min)/2) + +print ("==== 0 - for zero-length and 1 - for axis field RF gap ====") +st = "CavPhase[deg] deltaE_0[MeV] deltaE_1[MeV] delta_0-1[MeV] delta_phase_0-1[deg]" +print (st) + +for ind in range(Nph+1): + cav_phase = ind*phase_step + base_rf_gap.getRF_Cavity().setPhase(cav_phase) + three_point_gap.getRF_Cavity().setPhase(cav_phase) + + bunch = Bunch() + bunch_init.copyBunchTo(bunch) + drift_before.trackBunch(bunch) + base_rf_gap.trackBunch(bunch) + drift_after.trackBunch(bunch) + delta_eKin_out_0 = bunch.getSyncParticle().kinEnergy() - Ekin_init + exit_phase0 = (bunch.getSyncParticle().time()*cav_frequency)*180./math.pi + exit_phase0 = phaseNearTargetPhaseDeg(exit_phase0,0.) + + bunch = Bunch() + bunch_init.copyBunchTo(bunch) + three_point_gap.trackBunch(bunch) + delta_eKin_out_1 = bunch.getSyncParticle().kinEnergy() - Ekin_init + exit_phase1 = (bunch.getSyncParticle().time()*cav_frequency)*180./math.pi + exit_phase1 = phaseNearTargetPhaseDeg(exit_phase1,exit_phase0) + + delta_phase = exit_phase0 - exit_phase1 + + st = " %+8.2f "%(cav_phase*180./math.pi) + st += " %+8.4f "%(delta_eKin_out_0*1000.) + st += " %+8.4f "%(delta_eKin_out_1*1000.) + st += " %+8.4f "%((delta_eKin_out_0 - delta_eKin_out_1)*1000.) + st += " %+8.4f "%delta_phase + + print (st) + + +print ("Stop.") +sys.exit(0) \ No newline at end of file