source: MML/trunk/mml/setorbitbump.m @ 4

Last change on this file since 4 was 4, checked in by zhangj, 11 years ago

Initial import--MML version from SOLEIL@2013
File size: 17.4 KB
Line 
1function [OCS, OCS0, V, S, ErrorFlag] = setorbitbump(varargin)
2%SETORBITBUMP - Local bump program (uses setorbit)
3%  [OCS, OCS0, V, S] = setorbitbump(BPMFamily, BPMDeviceList, GoalOrbit, CMFamily, CMIncrementList, NIter, SVDIndex, BPMWeight)
4%  [OCS, OCS0, V, S] = setorbitbump(OCS, GoalOrbit, CMIncrementList, NIter, SVDIndex, BPMWeight)
5%
6%  INPUTS
7%  1. BPMFamily - BPM family
8%  2. BPMDeviceList - BPM device list
9%  3. GoalOrbit - Goal orbit position (Referenced to an absolute or incremental orbit
10%                                      depending on 'Incremental' or 'Absolute' flag)
11%  4. CMFamily - Corrector magnet family
12%  5. CMIncrementList - Magnet list upstream and downstream of the BPMs
13%                       For instance, [-2 1 3] uses the 2nd correctors upstream and
14%                       1st and 3rd corrector downstream from the BPMs in the bump.
15%          or
16%     CMIncrementList can be an ordinary device list used for specifying correctors.
17%  6. NIter - Number of iterations  {Default: 1}
18%            (set NIter to 0 if no BPMs are functioning, see setorbit)
19%  7. SVDIndex - vector, maximum number, 'All', or
20%                base on a threshold of min/max singular value {Default: see setorbit}
21%  8. BPMWeight - Weight applied to the BPMs inside the bump (GoalOrbit).  BPM weighting
22%                 should only be required when the corrector are in a region of dispersion
23%                 and the RF frequency is not used in the orbit correction.
24%                {Default BPMWeight is to weight the leakage twice as much as the GoalOrbit}
25%  9. Optional flags (any flags used in setorbit can also be used here)
26%     'Incremental' {Default} or 'Absolute' - Type of orbit change
27%     'Display' - Display bump characteristics  (no display is the default)
28%     'SetPV'   or 'SetSP'   - Make the setpoint change {Default}
29%     'NoSetPV' or 'NoSetSP' - Don't set the magnets, just return the OCS
30%     'FitRF' or 'SetRF' - Flag to include the RF frequency as part of orbit correction.
31%     'LeakageCorrection' - Option to try to clean up the leakage after the bump is applied.
32%                           Note: 1. NIter for leakage correction is the same for the bump
33%                                 2. It's just like setting the weight to zero for a second set of iterations.
34%                                 3. For the horizontal plane, usually one should include the RF frequency
35%                                    along with leakage correction.  Otherwise, the goal orbit may not be attained.
36%     'Tolerance', Tol - Quit correction if the std(error) < Tol {Tol: 0}
37%     'RampSteps', NSteps - Number of steps to ramp in the correction {NSteps: 1}
38%                                    This is used to avoid a large transient between setpoint changes.
39%
40%  OUTPUTS (same as setorbit)
41%  1. OCS    - Orbit correction structure (OCS) usable by setorbit
42%  2. OCS0   - Starting OCS
43%  3. V      - Corrector magnet singular vectors
44%  4. Svalue - Singular values (vector)
45%
46%  NOTES
47%  1. setorbitbump creates an OCS structure and uses setorbit to actually change the orbit.
48%
49%  EXAMPLES
50%  1. 4 magnet horizontal incremental bump at BPM(6,6) and BPM(7,1) including the RF frequency
51%     OCS = setorbitbump('BPMx',[6 6;7 1],[1;-1],'HCM',[-2 -1 1 2],'FitRF', 'Display');
52%     On 12-16-2003, [-2 -1 1 2] corresponded to HCM(6,3), HCM(6,4), HCM(7,1), HCM(7,3)
53%     Display plots information before applying the bump.
54%
55%  2. If the response matrix is not perfect, a clean (minimal rms orbit error) bump can be created
56%     by increasing the iteration (note: in regions of dispersion, all bumps will have leakage unless
57%     the RF frequency is included). 
58%     [OCS, OCS0] = setorbitbump('BPMx',[6 6;7 1],[1;-1],'HCM',[-2 -1 1 2], 4, 'Display');
59%     DeltaHCM = OCS.CM.Data - OCS0.CM.Data;
60%     BPMxDeviceList = OCS.BPM.DeviceList;
61%     HCMDeviceList  = OCS.CM.DeviceList;
62%
63%
64%  See also setorbit, orbitcorrectionmethods, setorbitgui, rmdisp, plotcm, setorbitbumpgui
65%
66%  Written by Greg Portmann
67
68
69%%%%%%%%%%%%
70% Defaults %
71%%%%%%%%%%%%
72OCS.BPM.FamilyName = [];
73OCS.BPM.DeviceList = [];
74GoalOrbit = [];
75OCS.CM.FamilyName = [];
76OCS.CM.DeviceList = [];
77CMIncrementList = [];
78OCS.NIter = 1;
79BPMWeight = [];
80OCS.Incremental = 1;
81DisplayFlag = 0;
82SetSPFlag = 1;
83NoLeakageFlag = 0;
84UnitsFlag = '';
85ModeFlag = '';
86
87L = getfamilydata('Circumference');
88
89
90%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91% Input parsing and checking %
92%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93Flags = {};
94for i = length(varargin):-1:1
95    if isstruct(varargin{i})
96        % Ignor structures
97    elseif iscell(varargin{i})
98        % Ignor cells
99    elseif strcmpi(varargin{i},'Golden')
100        % Use the golden orbit
101        %GoalOrbit = 'Golden';
102        %varargin(i) = [];
103    elseif strcmpi(varargin{i},'Offset')
104        % Use the offset orbit
105        %GoalOrbit = 'Offset';
106        %varargin(i) = [];
107    elseif strcmpi(varargin{i},'Display')
108        Flags = [Flags varargin(i)];
109        DisplayFlag = 1;
110        varargin(i) = [];
111    elseif strcmpi(varargin{i},'NoDisplay')
112        Flags = [Flags varargin(i)];
113        DisplayFlag = 0;
114        varargin(i) = [];
115    elseif strcmpi(varargin{i},'Inc') || strcmpi(varargin{i},'Incremental')
116        OCS.Incremental = 1;
117        varargin(i) = [];
118    elseif strcmpi(varargin{i},'Abs') || strcmpi(varargin{i},'Absolute')
119        OCS.Incremental = 0;
120        varargin(i) = [];
121    elseif strcmpi(varargin{i},'simulator') || strcmpi(varargin{i},'model')
122        ModeFlag = 'SIMULATOR';
123        varargin(i) = [];
124    elseif strcmpi(varargin{i},'Online')
125        ModeFlag = 'Online';
126        varargin(i) = [];
127    elseif strcmpi(varargin{i},'Manual')
128        ModeFlag = 'Manual';
129        varargin(i) = [];
130    elseif strcmpi(varargin{i},'physics')
131        UnitsFlag = 'Physics';
132        varargin(i) = [];
133    elseif strcmpi(varargin{i},'hardware')
134        UnitsFlag = 'Hardware';
135        varargin(i) = [];
136    elseif any(strcmpi(varargin{i},{'SetPV','SetSP'}))
137        SetSPFlag = 1;
138        varargin(i) = [];
139    elseif any(strcmpi(varargin{i},{'NoSetPV','NoSetSP'}))
140        SetSPFlag = 0;
141        Flags = [Flags varargin(i)];
142        varargin(i) = [];
143    elseif strcmpi(varargin{i},'LeakageCorrection')       
144        NoLeakageFlag = 1;
145        varargin(i) = [];
146       
147    elseif strcmpi(varargin{i},'archive') || strcmpi(varargin{i},'noarchive')
148        % Just remove
149        varargin(i) = [];
150    elseif strcmpi(varargin{i},'struct') || strcmpi(varargin{i},'numeric')
151        % Just remove
152        varargin(i) = [];
153
154    elseif strcmpi(varargin{i},'CorrectorGain')
155        Flags = [Flags varargin(i) varargin(i+1)];
156        varargin(i+1) = [];
157        varargin(i)   = [];
158
159    elseif strcmpi(varargin{i},'Tolerance')
160        Flags = [Flags varargin(i) varargin(i+1)];
161        varargin(i+1) = [];
162        varargin(i)   = [];
163
164    elseif strcmpi(varargin{i},'RampSteps')
165        Flags = [Flags varargin(i) varargin(i+1)];
166        varargin(i+1) = [];
167        varargin(i)   = [];
168
169    elseif strcmpi(varargin{i},'ModelResp')
170        Flags = [Flags varargin(i)];
171        varargin(i) = [];
172    elseif strcmpi(varargin{i},'ModelDisp')
173        Flags = [Flags varargin(i)];
174        varargin(i) = [];
175    elseif strcmpi(varargin{i},'MeasDisp')
176        Flags = [Flags varargin(i)];
177        varargin(i) = [];
178    elseif strcmpi(varargin{i},'GoldenDisp')
179        Flags = [Flags varargin(i)];
180        varargin(i) = [];
181       
182    elseif any(strcmpi(varargin{i},{'FitRF','SetRF','RFCorrector'}))
183        Flags = [Flags varargin(i)];
184        varargin(i) = [];
185    elseif any(strcmpi(varargin{i},{'FitRFHCM0','FitRFEnergy','SetRFHCM0','SetRFEnergy'}))
186        Flags = [Flags varargin(i)];
187        varargin(i) = [];
188    elseif any(strcmpi(varargin{i},{'FitRFDeltaHCM0','SetRFDeltaHCM0'}))
189        Flags = [Flags varargin(i)];
190        varargin(i) = [];
191    elseif any(strcmpi(varargin{i},{'FitRFDispersionOrbit','SetRFDispersionOrbit'}))
192        Flags = [Flags varargin(i)];
193        varargin(i) = [];
194    elseif any(strcmpi(varargin{i},{'FitRFDispersionHCM','SetRFDispersionHCM'}))
195        Flags = [Flags varargin(i)];
196        varargin(i) = [];
197
198    end
199end
200
201
202if length(varargin) >= 1
203    if isstruct(varargin{1})
204        OCS = varargin{1};
205        varargin(1) = [];
206        if length(varargin) < 2
207            error('An OCS plus at least 2 inputs are required');
208        end
209        if length(varargin) >= 1
210            GoalOrbit = varargin{1};
211            varargin(1) = [];
212        end
213        if length(varargin) >= 1
214            CMIncrementList = varargin{1};
215            varargin(1) = [];
216        end
217    else
218        OCS.BPM.FamilyName = varargin{1};
219        varargin(1) = [];
220        if length(varargin) >= 1
221            OCS.BPM.DeviceList = varargin{1};
222            varargin(1) = [];
223        end
224        if length(varargin) >= 1
225            GoalOrbit = varargin{1};
226            varargin(1) = [];
227        end
228        if length(varargin) >= 1
229            OCS.CM.FamilyName = varargin{1};
230            varargin(1) = [];
231        end
232        if length(varargin) >= 1
233            CMIncrementList = varargin{1};
234            varargin(1) = [];
235        end
236    end
237   
238    % Get NIter
239    if length(varargin) >= 1
240        if isnumeric(varargin{1})
241            OCS.NIter = varargin{1};
242            varargin(1) = [];
243        end
244    end
245   
246    % Get SVDIndex (can be a fractional number, vector, or 'All')
247    if length(varargin) >= 1
248        %if isnumeric(varargin{1})
249            OCS.SVDIndex = varargin{1};
250            varargin(1) = [];
251        %end
252    end
253
254    % Get BPMWeight
255    if length(varargin) >= 1
256        if isnumeric(varargin{1})
257            BPMWeight = varargin{1};
258            varargin(1) = [];
259        end
260    end
261end
262
263
264% Pass the extra inputs on
265if length(varargin) >= 1
266    Flags = [Flags varargin];
267end
268
269
270% Fill the empty fields
271if isempty(OCS.BPM.FamilyName)
272    i = menu('Select a plane','Horizontal','Vertical');
273    if i == 1
274        OCS.BPM.FamilyName = gethbpmfamily;
275        OCS.CM.FamilyName  = gethcmfamily;
276    elseif i == 2
277        OCS.BPM.FamilyName = getvbpmfamily;
278        OCS.CM.FamilyName  = getvcmfamily;
279    end
280end
281if isempty(OCS.BPM.DeviceList)
282    OCS.BPM.DeviceList = editlist(family2dev(OCS.BPM.FamilyName), OCS.BPM.FamilyName, 0);
283    % Getting a list can messes up the order for bumps going from the last sector to the first
284    for i = 2:length(OCS.BPM.DeviceList)
285        if getspos(OCS.BPM.FamilyName,OCS.BPM.DeviceList(i,:))-getspos(OCS.BPM.FamilyName,OCS.BPM.DeviceList(i-1,:)) > L/2
286            break;
287        end
288    end
289    if i < length(OCS.BPM.DeviceList)
290        OCS.BPM.DeviceList = [OCS.BPM.DeviceList(i:end,:); OCS.BPM.DeviceList(1:i-1,:)];
291    end
292end
293if isempty(OCS.CM.FamilyName)
294    if strcmpi(OCS.BPM.FamilyName, gethbpmfamily)
295        OCS.CM.FamilyName = gethcmfamily;
296    else
297        OCS.CM.FamilyName = getvcmfamily;
298    end
299end
300
301if isempty(UnitsFlag)
302    UnitsFlag = getunits(OCS.BPM.FamilyName);
303end
304
305if isempty(ModeFlag)
306    ModeFlag = getmode(OCS.BPM.FamilyName);
307end
308
309if NoLeakageFlag && OCS.NIter==0
310    error('Leakage cleanup flag cannot be used with zero iterations.');
311end
312%%%%%%%%%%%%%%%%%%%%%%%%
313% End of input parsing %
314%%%%%%%%%%%%%%%%%%%%%%%%
315
316
317% Possibly convert the GoalOrbit
318if any(isnan(GoalOrbit))
319    error('Goal orbit has a NaN');
320end
321if isempty(GoalOrbit)
322    % Input the goal orbit
323    for i = 1:size(OCS.BPM.DeviceList,1)
324        labels{i} = sprintf('Input the goal orbit at %s(%d,%d)',OCS.BPM.FamilyName, OCS.BPM.DeviceList(i,1), OCS.BPM.DeviceList(i,2));
325        StartingPoint{i} = '0';
326    end
327    answer = inputdlg(labels,'Local Bump',1,StartingPoint);
328    if isempty(answer)
329        fprintf('   Local bump cancelled');
330        return
331    end
332    for i = 1:size(OCS.BPM.DeviceList,1)
333        GoalOrbit(i,1) = str2num(answer{i});
334    end
335end
336if ischar(GoalOrbit)
337    if strcmpi(GoalOrbit, 'Golden')
338        GoalOrbit = getgolden(OCS.BPM.FamilyName, OCS.BPM.DeviceList);
339    elseif strcmpi(GoalOrbit, 'Offset')
340        GoalOrbit = getoffset(OCS.BPM.FamilyName, OCS.BPM.DeviceList);
341    else
342        error('Goal unknown');
343    end
344end
345GoalOrbit = GoalOrbit(:);
346
347% Check the length
348if length(GoalOrbit) ~= size(OCS.BPM.DeviceList,1)
349    error('Length of the GoalOrbit must equal the number of devices in BPMList');
350end
351
352
353% CMIncrementList must be a vector, no zeros, no repeats
354CheckCM = 0;
355if isempty(CMIncrementList)
356    CheckCM = 1;
357    Ncm = size(OCS.BPM.DeviceList,1)+2;
358    if rem(Ncm,2) == 0
359        NcmDiv2 = Ncm / 2;
360        CMIncrementList = [-(Ncm/2):-1 1:(Ncm/2)];
361    else
362        NcmDiv2 = Ncm / 2;
363        CMIncrementList = [-floor(Ncm/2):-1 1:(floor(Ncm/2)+1)];
364    end
365end
366
367% Check for a device list input
368if ~(size(CMIncrementList,2) == 2 && size(CMIncrementList,1) > 1)
369    CMIncrementList = CMIncrementList(:);
370    CMIncrementList = sort(CMIncrementList);
371    CMIncrementList(find(CMIncrementList==0)) = [];
372    CMIncrementList(find(diff(CMIncrementList)==0)) = [];
373end
374
375% Get BPM positions
376BPMListTotal = family2dev(OCS.BPM.FamilyName, 1);
377BPMsposTotal = getspos(OCS.BPM.FamilyName, BPMListTotal);
378BPMspos      = getspos(OCS.BPM.FamilyName, OCS.BPM.DeviceList);
379
380
381% Stack 3 rings so you you don't have to worry about the L to 0 transition
382CMListTotal = family2dev(OCS.CM.FamilyName, 1);
383CMsposTotal = getspos(OCS.CM.FamilyName, CMListTotal);
384CMListTotal = [CMListTotal; CMListTotal; CMListTotal];
385CMsposTotal = [CMsposTotal-L; CMsposTotal; CMsposTotal+L];
386
387
388% Find the correctors
389% Check for a device list input
390if size(CMIncrementList,2) == 2 && size(CMIncrementList,1) > 1
391    % DeviceList
392    OCS.CM.DeviceList = CMIncrementList;
393else
394    for i = 1:length(CMIncrementList)
395        if CMIncrementList(i) <= 0
396            j = find(CMsposTotal <= BPMspos(1));
397            OCS.CM.DeviceList(i,:) = CMListTotal(j(end)+CMIncrementList(i)+1,:);
398        else
399            j = find(CMsposTotal >= BPMspos(end));
400            OCS.CM.DeviceList(i,:) = CMListTotal(j(1)+CMIncrementList(i)-1,:);
401        end
402    end
403end
404
405if CheckCM
406    DevTotal = family2dev(OCS.CM.FamilyName);
407    checkvec = zeros(size(DevTotal,1),1);
408    checkvec(findrowindex(OCS.CM.DeviceList, DevTotal)) = 1;
409    OCS.CM.DeviceList = editlist(DevTotal, OCS.CM.FamilyName, checkvec);
410
411    % Getting a list can messes up the order for bumps going from the last sector to the first
412    for i = 2:length(OCS.CM.DeviceList)
413        if getspos(OCS.CM.FamilyName,OCS.CM.DeviceList(i,:)) - getspos(OCS.CM.FamilyName,OCS.CM.DeviceList(i-1,:)) > L/2
414            break;
415        end
416    end
417    if i < length(OCS.CM.DeviceList)
418        OCS.CM.DeviceList = [OCS.CM.DeviceList(i:end,:); OCS.CM.DeviceList(1:i-1,:)];
419    end
420end
421
422
423% Find all BPMs outside the bump (leakage control BPMs)
424CMspos = getspos(OCS.CM.FamilyName, OCS.CM.DeviceList);
425if CMspos(1) > CMspos(end)
426    j1 = intersect(find(BPMsposTotal < CMspos(1)), find(BPMsposTotal > CMspos(end)));
427    j2 = [];
428else
429    j1 = find(BPMsposTotal < CMspos(1));
430    j2 = find(BPMsposTotal > CMspos(end));
431end
432if isempty(j1) && isempty(j2)
433    error('Cound not find any leakage control BPMs');
434end
435OCS.BPM.DeviceList   = [BPMListTotal(j1,:); OCS.BPM.DeviceList; BPMListTotal(j2,:);];
436BPMDeviceListLeakage = [BPMListTotal(j1,:);                     BPMListTotal(j2,:);];
437
438
439if OCS.Incremental
440    % Incremental
441    OCS.GoalOrbit = [zeros(length(j1),1); GoalOrbit; zeros(length(j2),1)];
442
443    if OCS.NIter == 0
444        LeakageGoalOrbit = zeros(size(BPMDeviceListLeakage,1),1);
445    else
446        if isempty(j1)
447            LeakageOrbit1 = [];
448        else
449            LeakageOrbit1 = getam(OCS.BPM.FamilyName,BPMListTotal(j1,:), UnitsFlag, ModeFlag);
450        end
451        if isempty(j2)
452            LeakageOrbit2 = [];
453        else
454            LeakageOrbit2 = getam(OCS.BPM.FamilyName,BPMListTotal(j2,:), UnitsFlag, ModeFlag);
455        end
456        LeakageGoalOrbit = [LeakageOrbit1; LeakageOrbit2];
457    end
458elseif ~OCS.Incremental
459    % Absolute
460    if isempty(j1)
461        LeakageOrbit1 = [];
462    else
463        LeakageOrbit1 = getam(OCS.BPM.FamilyName,BPMListTotal(j1,:), UnitsFlag, ModeFlag);
464    end
465    if isempty(j2)
466        LeakageOrbit2 = [];
467    else
468        LeakageOrbit2 = getam(OCS.BPM.FamilyName,BPMListTotal(j2,:), UnitsFlag, ModeFlag);
469    end
470    OCS.GoalOrbit    = [LeakageOrbit1; GoalOrbit; LeakageOrbit2];
471    LeakageGoalOrbit = [LeakageOrbit1; LeakageOrbit2];
472else
473    error('Absolute or incremental input unknown.');
474end
475
476
477% Add a weight
478if isempty(BPMWeight)
479    % Default BPMWeight is to weight the leakage twice as much as the bump goal
480    BPMWeight = .5 * (length(j1) + length(j2)) / length(GoalOrbit);
481    OCS.BPMWeight = [ones(length(j1),1); BPMWeight .* ones(length(GoalOrbit),1); ones(length(j2),1)];
482else
483    OCS.BPMWeight = [ones(length(j1),1); BPMWeight .* ones(length(GoalOrbit),1); ones(length(j2),1)];
484end
485
486
487% Just to fill the data structure fields
488OCS.BPM = family2datastruct(OCS.BPM.FamilyName, 'Monitor',  OCS.BPM.DeviceList, UnitsFlag, ModeFlag);
489OCS.CM  = family2datastruct(OCS.CM.FamilyName,  'Setpoint', OCS.CM.DeviceList,  UnitsFlag, ModeFlag);
490
491
492% Bumps stats
493if DisplayFlag
494    % Corrector strength (meters/radian) and (mm/amp)
495    % Warn on when to add RF frequency (generated dispersion is greater
496    % than mm per amp or % of mm/amp bump)   
497end
498
499
500% Set the bump
501[OCS, OCS0, V, S, ErrorFlag] = setorbit(OCS, UnitsFlag, ModeFlag, Flags{:});
502
503
504% Clean up the leakage
505if NoLeakageFlag
506    if DisplayFlag
507        fprintf('   Correcting the leakage\n');
508    end
509    OCS1 = OCS;
510    OCS1.GoalOrbit = LeakageGoalOrbit;
511    OCS1.BPM = family2datastruct(OCS.BPM.FamilyName, 'Monitor', BPMDeviceListLeakage);
512    OCS1.BPMWeight = [];
513    OCS1.Incremental = 0;
514    OCS1 = setorbit(OCS1, UnitsFlag, ModeFlag);
515end
516
Note: See TracBrowser for help on using the repository browser.