forked from ThorstenHellert/SC
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SCsynchPhaseCorrection.m
291 lines (251 loc) · 8.44 KB
/
SCsynchPhaseCorrection.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
function [deltaPhi,ERROR] = SCsynchPhaseCorrection(SC,varargin)
% SCsynchPhaseCorrection
% ======================
%
% NAME
% ----
% SCsynchPhaseCorrection - Calculates a beam based correction to the rf cavity phase
%
% SYNOPSIS
% --------
% `[deltaPhi, ERROR] = SCsynchPhaseCorrection(SC [, options])`
%
%
% DESCRIPTION
% -----------
%
% Changes the cavity phase within the phase interval [-pi,pi] stepwise and evaluates the mean
% turn-by-turn horizontal BPM deviation. A sinusoidal function is fitted to the data and the zero
% crossing is identified. It is assumed that for sufficiently small number of turns the synchrotron
% motion is negligible and injection at the synchronous phase will result in zero turn-by-turn
% energy variation. Thus, the horizontal turn-by-turn BPM readings should in first approximation not
% differ. The number of evaluated turns should be significantly smaller than half a sunchrotron
% period. Note that if more than one cavity is specified the same phase steps are applied to all
% cavities. Results might be compromised.
%
%
% INPUTS
% ------
% `SC`::
% The SC base structure
%
%
% OPTIONS
% -------
% The following options can be given as name/value-pairs:
%
% `'cavOrd'` (`SC.ORD.Cavity`)::
% Ordinate of evaluated cavity
% `'nSteps'` (`30`)::
% Number of phase steps to be evaluated
% `'nTurns'` (`30`)::
% Number of turns to be evaluated
% `'BPMords'` (`SC.ORD.BPM`)::
% BPM ordinates where the dispersive offset change is evaluated
% `'plotResults'` (`0`)::
% If true, results are plotted.
% `'plotProgress'` (`0`)::
% If true, each phase step is plotted.
% `'verbose'` (`0`)::
% If true, debug information is printed.
%
%
% RETURN VALUES
% -------------
% `deltaPhi`::
% Phase correction step to be added to cavity field 'TimeLag'.
% `ERROR`::
% Error value.
%
% ERRORS
% ------
% `0`::
% All good.
% `1`::
% Horizontal TBT BPM deviation shows no zero crossing
% `2`::
% Sinusoidal fit function shows no zero crossing
%
%
% SEE ALSO
% --------
% *SCsetCavs2SetPoints*, *SCgetBPMreading*, *SCsynchEnergyCorrection*
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input check
% Parse input
p = inputParser;
addOptional(p,'cavOrd',SC.ORD.Cavity);
addOptional(p,'nSteps',15);
addOptional(p,'nTurns',20);
addOptional(p,'BPMords',SC.ORD.BPM);
addOptional(p,'plotResults',0);
addOptional(p,'plotProgress',0);
addOptional(p,'verbose',0);
parse(p,varargin{:});
par = p.Results;
inputCheck(SC,par);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initialization
% Initialize output
ERROR = 0;
deltaPhi = 0;
% Initialize TbT BPM shift
BPMshift = nan(1,par.nSteps);
% RF wavelength [m]
lambda = 299792458/SC.RING{par.cavOrd}.Frequency;
% Define phase scan range
lambdaTestVec = 1/2 * lambda * linspace(-1,1,par.nSteps);
% Adjust number of turns
SC.INJ.nTurns = par.nTurns;
if par.verbose
fprintf('Calibrate RF phase with: \n %d Particles \n %d Turns \n %d Shots \n %d Phase steps.\n\n',SC.INJ.nParticles,SC.INJ.nTurns,SC.INJ.nShots,par.nSteps)
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Main script
% Loop over different phase offsets
for nL = 1:length(lambdaTestVec)
% Change phase
tmpSC = SCsetCavs2SetPoints(SC,par.cavOrd,'TimeLag',lambdaTestVec(nL),'add');
% Calculate turn-by-turn horizontal trajectory shift
[BPMshift(nL),TBTdE] = getTbTEnergyShift(tmpSC,par);
% Plot current step
if par.plotProgress
plotProgress(TBTdE,BPMshift,lambdaTestVec,nL)
end
end
% Check if zero crossing is reached in data
if ~(max(BPMshift)>0 && min(BPMshift)<0)
fprintf('Zero crossing not within data set.\n')
ERROR = 1;
return
end
% Define sinusoidal function
sinFun = @(par,s) par(1)*sin(2*pi*(par(4)*s + par(2))) +par(3);
% Define merit function
fomFun = @(par) sum( (sinFun(par,lambdaTestVec) - BPMshift).^2 );
% Loop over different start point guesses
for startPhase=[-pi,-pi/2,-pi/4,0]
% Run fminsearch to determine sinusoidal parameters
sol = fminsearch(fomFun, [max(BPMshift)-mean(BPMshift); startPhase; mean(BPMshift) ; 0.5/max(lambdaTestVec)]);
% Save solution in sinusoidal function
solFun = @(x) sinFun(sol,x);
% Generate densely populated lambda-vector
xp = linspace(lambdaTestVec(1),lambdaTestVec(end),100);
% Check if zero crossing is reached in fitted function
if ~(max(solFun(xp))>0 && min(solFun(xp))<0)
fprintf('Zero crossing not within fit function, trying different start point guess.\n')
else
break
end
end
% Check if zero crossing is reached in fitted function
if ~(max(solFun(xp))>0 && min(solFun(xp))<0)
fprintf('Zero crossing not within fit function\n')
ERROR = 2;
% return
end
% Identify maximum value
[~,maxValInd] = max(solFun(xp));
% Find zero crossing on the left hand side of maximum value
deltaPhi = fzero(solFun,xp(maxValInd)-abs(lambdaTestVec(1))/2);
% Shift phase into [-pi,pi]-intervall
if deltaPhi > lambda/2
deltaPhi = deltaPhi - lambda;
elseif deltaPhi < -lambda/2
deltaPhi = deltaPhi + lambda;
end
% Plot results
if par.plotResults
plotFunction()
end
% Check for NaN results
if isnan(deltaPhi)
ERROR = 3;
fprintf('SCrfCommissioning: ERROR (NaN phase)\n')
return
end
% Print results
if par.verbose
% Initial phase
XCO = findorbit6(SC.RING);
% Apply phase correction
tmpSC = SCsetCavs2SetPoints(SC,par.cavOrd,'TimeLag',deltaPhi,'add');
% Calculate corrected closed orbit
XCOfinal = findorbit6(tmpSC.RING);
% Calculate figures of merrit
initial = rem((XCO(6) -SC.INJ.Z0(6))/lambda*360,360);
final = rem((XCOfinal(6)-SC.INJ.Z0(6))/lambda*360,360);
% Print result
fprintf('Phase correction step: %.3m\n',deltaPhi)
fprintf('>> Static phase error corrected from %.0fdeg to %.1fdeg\n',initial,final)
end
% Plot final results
function plotFunction
figure(87);clf;hold on
% Plot TBT BPM shift
plot((lambdaTestVec+SC.INJ.Z0(6))/lambda*360,1E6*BPMshift,'o')
% Plot fitted sinosoidal functions
plot((xp+SC.INJ.Z0(6))/lambda*360,1E6* solFun(xp))
% Plot initial phase
plot((SC.INJ.Z0(6))/lambda*2*180,1E6* (sol(1)*sin(2*pi*(sol(4)*0 + sol(2))) +sol(3)),'rD','MarkerSize',16)
% Plot final phase
plot((deltaPhi)/lambda*2*180,0,'kX','MarkerSize',16)
set(gca,'XLim',180*[-1 1],'box','on');
legend({'Measurement','Fit','Initial phase','Phase correction'})
xlabel('RF phase [$^\circ$]');ylabel('$<\Delta x>$ [$\mu$m/turn]');
set(findall(gcf,'-property','FontSize'),'FontSize',18);
set(findall(gcf,'-property','Interpreter'),'Interpreter','latex');
set(findall(gcf,'-property','TickLabelInterpreter'),'TickLabelInterpreter','latex');
set(gcf,'color','w');drawnow
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary functions
% Input check
function inputCheck(SC,par)
if ~strcmp(SC.INJ.trackMode,'TBT')
error('Trackmode should be turn-by-turn (''SC.INJ.trackmode=TBT'').')
end
if par.nSteps<2 || par.nTurns<2
error('Number of steps and number of turns must be larger than 2.')
end
if ~isfield(SC.RING{par.cavOrd},'Frequency')
error('This is not a cavity (ord: %d)',par.cavOrd)
end
if ~any(strcmp(SC.RING{par.cavOrd}.PassMethod,{'CavityPass','RFCavityPass'}))
warning('Cavity (ord: %d) seemed to be switched off.',par.cavOrd)
end
end
% Calculate mean turn-by-turn horizontal BPM shift
function [BPMshift,dE] = getTbTEnergyShift(SC,par)
% Calculate beam reading
B = SCgetBPMreading(SC,par.BPMords);
% Reshape horizontal BPM readings turn-by-turn
BB = reshape(B(1,:),[],SC.INJ.nTurns);
% Calculate mean turn-by-turn difference
dE = mean( BB - repmat(BB(:,1),1,SC.INJ.nTurns) );
% Calculate mean turn-by-turn energy shift
x = (1:(SC.INJ.nTurns-1))';
y = dE(2:end)';
% Exclude nan BPM readings
x(isnan(y),:) = [];
y(isnan(y)) = [];
% Perform least squares fit
BPMshift= x\y;
end
% Plot current step
function plotProgress(TBTdE,BPMshift,fTestVec,nE)
figure(2);clf
subplot(2,1,1);hold on
plot(TBTdE,'o');
plot([1:length(TBTdE)] * BPMshift(nE),'--')
xlabel('Number of turns');ylabel('$<\Delta x_\mathrm{TBT}>$ [m]');
subplot(2,1,2);
plot(fTestVec(1:nE),BPMshift(1:nE),'o')
xlabel('$\Delta \phi$ [m]');ylabel('$<\Delta x>$ [m/turn]');
set(findall(gcf,'-property','FontSize'),'FontSize',18);
set(findall(gcf,'-property','Interpreter'),'Interpreter','latex');
set(findall(gcf,'-property','TickLabelInterpreter'),'TickLabelInterpreter','latex');
set(gcf,'color','w');drawnow
drawnow
end