source: MML/trunk/applications/loco/locoplot.m @ 4

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

Initial import--MML version from SOLEIL@2013
File size: 18.6 KB
Line 
1function ElementsInput = locoplot(FileName, IterationNumber, PlotString, PlaneString, ElementsInput)
2%LOCOPLOT - Subroutines for plotting LOCO output
3%  locoplot(FileName, IterationNumber, PlotString, PlaneString, Elements)
4%  locoplot({BPMData, CMData, LocoMeasData, LocoModel, FitParameters, LocoFlags, RINGData}, IterationNumber, PlotString, PlaneString, Elements)
5%
6%  INPUTS
7%  1. FileName = data file name
8%  2. IterationNumber = 0, 1, 2, etc
9%       or
10%  1. BPMData
11%  2. CMData
12%  3. LocoMeasData
13%  4. LocoModel
14%  5. FitParameters
15%  6. LocoFlags
16%  7. RINGData
17%
18%  What data to plot:
19%     PlotString = 'Model', 'Measured', 'Model-Measured', , 'Model-Measured+EnergyShifts'
20%
21%  For 3-D surf plots of the entire response matrix:
22%     PlaneString = 'XX', 'YY', 'XY', 'YX', 'All'
23%
24%  For 2-D row (BPM) or column (CM) plots:
25%     PlaneString = 'HBPM', 'VBPM', 'HCM', 'VCM'
26%     Elements = element number within the response matrix (based on HBPMIndex,VBPMIndex,HCMIndex,VCMIndex)
27%
28%  Written by Greg Portmann
29
30if ~nargin == 5
31    error('Requires 5 inputs (see help locoplot).');
32end
33
34if isempty(FileName)
35    return;
36end
37
38if iscell(FileName)
39    BPMData       = FileName{1};
40    CMData        = FileName{2};
41    LocoMeasData  = FileName{3};
42    LocoModel     = FileName{4};
43    FitParameters = FileName{5};
44    LocoFlags     = FileName{6};
45    RINGData      = FileName{7};
46elseif isstr(FileName)   
47    try
48        load(FileName);
49    catch
50        fprintf('   LOCOPLOT:  File does not exist or is not a *.mat file type.\n'); return;
51    end   
52else
53    error('Input problem');
54end
55
56legend off
57
58if length(BPMData) > 1
59    IterationNumber = IterationNumber + 1;
60    if IterationNumber > length(BPMData)
61        fprintf('   LOCOPLOT:  The data file only has %d iterations.  Hence, the input InterationNumber cannot be %d.\n', length(BPMData)-1, IterationNumber-1);
62        return
63    end
64   
65    BPMData = BPMData(IterationNumber);
66    CMData = CMData(IterationNumber);
67    LocoModel = LocoModel(IterationNumber);
68    FitParameters = FitParameters(IterationNumber);
69    LocoFlags = LocoFlags(IterationNumber);
70end
71
72% if isempty(LocoModel.M)
73%     fprintf('   LOCOPLOT:  Model data is empty.\n');
74%     return;
75% end
76if isempty(LocoMeasData.M)
77    fprintf('   LOCOPLOT:  Measured data is empty. \n');
78    return;
79end
80
81Mmodel = LocoModel.M;
82
83Outliers = LocoModel.OutlierIndex;
84
85Mmeas = LocoMeasData.M;
86Mmeas = Mmeas([BPMData.HBPMGoodDataIndex(:)' length(BPMData.HBPMIndex)+BPMData.VBPMGoodDataIndex(:)'], [CMData.HCMGoodDataIndex(:)' length(CMData.HCMIndex)+CMData.VCMGoodDataIndex(:)']);
87
88NHBPM = length(BPMData.HBPMGoodDataIndex);
89NVBPM = length(BPMData.VBPMGoodDataIndex);
90NBPM  = NHBPM + NVBPM;
91NHCM = length(CMData.HCMGoodDataIndex);
92NVCM = length(CMData.VCMGoodDataIndex);
93
94% Mark the outliers with NaN
95if isempty(Mmodel)
96    Mmodel = NaN * Mmeas;
97else
98    Mmodel = Mmodel(:);
99    Mmodel(Outliers) = NaN;
100    Mmodel = reshape(Mmodel,  NBPM, NHCM+NVCM);
101end
102
103if strcmp(lower(PlotString),'model')
104    data = Mmodel;
105elseif strcmp(lower(PlotString),'measured')
106    data = Mmeas;
107elseif strcmp(lower(PlotString),'model-measured')
108    data = Mmodel-Mmeas;
109elseif strcmp(lower(PlotString),'model-measured+energyshifts')
110    %[Model, Meas] = locodata(FileName, IterationNumber-1, 'ResponseMatrix', 'Model+EnergyShift', 'ResponseMatrix', 'Meas');
111   
112    % Add the dispersion term (measured) to the model response matrix
113    HCMEnergyShift = CMData.HCMEnergyShift(CMData.HCMGoodDataIndex);
114    VCMEnergyShift = CMData.VCMEnergyShift(CMData.VCMGoodDataIndex);
115   
116    % Set the lattice model
117    for j = 1:length(FitParameters.Params)
118        RINGData = locosetlatticeparam(RINGData, FitParameters.Params{j}, FitParameters.Values(j));
119    end
120    AlphaMCF = locomcf(RINGData);
121    EtaXmcf = -AlphaMCF * LocoMeasData.RF * LocoMeasData.Eta(BPMData.HBPMGoodDataIndex) / LocoMeasData.DeltaRF;
122    EtaYmcf = -AlphaMCF * LocoMeasData.RF * LocoMeasData.Eta(length(BPMData.HBPMIndex)+BPMData.VBPMGoodDataIndex) / LocoMeasData.DeltaRF;
123   
124    for i = 1:length(HCMEnergyShift)
125        Mmodel(:,i) = Mmodel(:,i) + HCMEnergyShift(i) * [EtaXmcf; EtaYmcf];
126    end
127   
128    for i = 1:length(VCMEnergyShift)
129        Mmodel(:,NHCM+i) = Mmodel(:,NHCM+i) + VCMEnergyShift(i) * [EtaXmcf; EtaYmcf];
130    end
131    data = Mmodel-Mmeas;
132else
133    error('   PlotString unknown type.');
134end
135
136PlaneString = upper(PlaneString);
137
138if strcmp(PlaneString,'XX') | strcmp(PlaneString,'YY') | strcmp(PlaneString,'YX') | strcmp(PlaneString,'XY') | strcmp(PlaneString,'ALL')
139    % 3-D Plots   
140    if strcmp(PlaneString,'XX')
141        % XX
142        [X,Y] = meshgrid(BPMData.HBPMGoodDataIndex, CMData.HCMGoodDataIndex);
143        surf(X, Y, data(1:NHBPM,1:NHCM)');
144        ylabel('HCM #');
145        xlabel('HBPM #');
146       
147    elseif strcmp(PlaneString,'YY')
148        % YY
149        [X,Y] = meshgrid(BPMData.VBPMGoodDataIndex, CMData.VCMGoodDataIndex);
150        surf(X, Y, data(NHBPM+1:NHBPM+NVBPM,NHCM+1:NHCM+NVCM)');
151        ylabel('VCM #');
152        xlabel('VBPM #');
153    elseif strcmp(PlaneString,'XY')
154        % XY
155        [X,Y] = meshgrid(BPMData.HBPMGoodDataIndex, CMData.VCMGoodDataIndex);
156        surf(X, Y, data(1:NHBPM,NHCM+1:NHCM+NVCM)');
157        ylabel('VCM #');
158        xlabel('HBPM #');
159       
160    elseif strcmp(PlaneString,'YX')
161        % YX
162        [X,Y] = meshgrid(BPMData.VBPMGoodDataIndex, CMData.HCMGoodDataIndex);
163        surf(X, Y, data(NHBPM+1:NHBPM+NVBPM,1:NHCM)');
164        ylabel('HCM #');
165        xlabel('VBPM #');
166    elseif strcmp(PlaneString,'ALL')
167        % All
168        surf(data');
169        %ylabel('CM''s');
170        %xlabel('BPM''s');
171       
172        % Change label to BPM numbers
173        [Ny,Nx] = size(data');
174        Ticks = 1:ceil(Nx/10):Nx;
175        set(gca, 'XTick', Ticks);     
176        TickNumber = [BPMData.HBPMGoodDataIndex(:)' BPMData.VBPMGoodDataIndex(:)'];
177        TickNumber = TickNumber(Ticks);
178        set(gca, 'XTickLabel', num2cell(TickNumber));         
179        xlabel('HBPM# and VBPM#');
180
181        % Change label to HCM and VCM numbers
182        Ticks = 1:ceil(Ny/10):Ny;
183        set(gca, 'YTick', Ticks);     
184        TickNumber = [CMData.HCMGoodDataIndex(:)' CMData.VCMGoodDataIndex(:)'];
185        TickNumber = TickNumber(Ticks);
186        set(gca, 'YTickLabel', num2cell(TickNumber));         
187        ylabel('HCM# and VCM#');
188    else
189        fprintf('   LOCOPLOT:  PlaneString unknown. \n');
190        return;       
191    end
192   
193   
194    if strcmp(lower(PlotString),'model')
195        title('Model Response Matrix');
196        zlabel('[mm]');
197    elseif strcmp(lower(PlotString),'measured')
198        title('Measured Response Matrix');
199        zlabel('[mm]');
200    elseif strcmp(lower(PlotString),'model-measured')
201        title('Model - Measured Response Matrix');
202        zlabel('Error [mm]');
203    elseif strcmp(lower(PlotString),'model-measured+energyshifts')
204        title('Model-Measured+EnergyShifts Response Matrix');
205        zlabel('Error [mm]');
206    end
207   
208    % For some reason axis tight has a problem with 3-D plots with only NaN's
209    if ~all(isnan(data))
210        axis tight
211    end
212end
213
214
215% 2-D plots
216
217if strcmp(PlaneString,'HBPM')
218    BPMData.BPMIndex = BPMData.BPMIndex(:);
219    if ~exist('ElementsInput','var')
220        ElementsInput = locoeditlist([BPMData.HBPMGoodDataIndex(:) BPMData.BPMIndex(BPMData.HBPMIndex(BPMData.HBPMGoodDataIndex))], PlaneString, zeros(length(BPMData.HBPMGoodDataIndex),1));
221        if isempty(ElementsInput)
222            return
223        end
224        ElementsInput = ElementsInput(:,1)';
225    end
226    if isempty(ElementsInput)
227        ElementsInput = locoeditlist([BPMData.HBPMGoodDataIndex(:) BPMData.BPMIndex(BPMData.HBPMIndex(BPMData.HBPMGoodDataIndex))], PlaneString, zeros(length(BPMData.HBPMGoodDataIndex),1));
228        if isempty(ElementsInput)
229            return
230        end
231        ElementsInput = ElementsInput(:,1)';
232    end
233   
234    % Convert to elements in the HBPMGoodDataIndex vector
235    [tmp1, tmp2, Elements] = intersect(ElementsInput, BPMData.HBPMGoodDataIndex);
236   
237    if ~isempty(Elements)
238        plot(data(Elements,:)');
239       
240        % Change label to HCM and VCM numbers
241        N = NHCM + NVCM;
242        Ticks = 1:ceil(N/10):N;
243        set(gca, 'XTick', Ticks);     
244        TickNumber = [CMData.HCMGoodDataIndex(:)' CMData.VCMGoodDataIndex(:)'];
245        TickNumber = TickNumber(Ticks);
246        set(gca, 'XTickLabel', num2cell(TickNumber));         
247        xlabel('Horizontal and Vertical Corrector Magnets');
248       
249        if strcmp(lower(PlotString),'model')
250            title('Model Response Matrix');
251            if length(Elements) == 1
252                ylabel(sprintf('Horizontal BPM(%d) [mm]', BPMData.HBPMGoodDataIndex(Elements)));
253            else
254                ylabel('Horizontal BPM [mm]');
255            end
256        elseif strcmp(lower(PlotString),'measured')
257            title('Measured Response Matrix');
258            if length(Elements) == 1
259                ylabel(sprintf('Horizontal BPM(%d) [mm]', BPMData.HBPMGoodDataIndex(Elements)));
260            else
261                ylabel('Horizontal BPM [mm]');
262            end
263        elseif strcmp(lower(PlotString),'model-measured')
264            title('Model - Measured Response Matrix');
265            if length(Elements) == 1
266                ylabel(sprintf('Horizontal Error BPM(%d) [mm]', BPMData.HBPMGoodDataIndex(Elements)));
267            else
268                ylabel('Horizontal BPM Error [mm]');
269            end
270        elseif strcmp(lower(PlotString),'model-measured+energyshifts')
271            title('Model-Measured+EnergyShifts Response Matrix');
272            if length(Elements) == 1
273                ylabel(sprintf('Horizontal Error BPM(%d) [mm]', BPMData.HBPMGoodDataIndex(Elements)));
274            else
275                ylabel('Horizontal BPM Error [mm]');
276            end
277        end
278       
279        % Only put a legend up to size 10
280        if length(Elements)>1 & length(Elements)<=10
281            for i = 1:length(Elements)
282                legendstr{i} = sprintf('BPM(%d)',BPMData.HBPMGoodDataIndex(Elements(i)));
283            end
284            legend(legendstr,0);
285        end
286        axis tight
287    end
288end
289
290
291if strcmp(PlaneString,'VBPM')
292    BPMData.BPMIndex = BPMData.BPMIndex(:);
293    if ~exist('ElementsInput','var')
294        ElementsInput = locoeditlist([BPMData.VBPMGoodDataIndex(:) BPMData.BPMIndex(BPMData.VBPMIndex(BPMData.VBPMGoodDataIndex))], PlaneString, zeros(length(BPMData.VBPMGoodDataIndex),1));
295        if isempty(ElementsInput)
296            return
297        end
298        ElementsInput = ElementsInput(:,1)';
299    end
300    if isempty(ElementsInput)
301        ElementsInput = locoeditlist([BPMData.VBPMGoodDataIndex(:) BPMData.BPMIndex(BPMData.VBPMIndex(BPMData.VBPMGoodDataIndex))], PlaneString, zeros(length(BPMData.VBPMGoodDataIndex),1));
302        if isempty(ElementsInput)
303            return
304        end
305        ElementsInput = ElementsInput(:,1)';
306    end
307   
308    % Convert to elements in the HBPMGoodDataIndex vector
309    [tmp1, tmp2, Elements] = intersect(ElementsInput, BPMData.VBPMGoodDataIndex);
310   
311    if ~isempty(Elements)
312        plot(data(NHBPM+Elements,:)');
313       
314        % Change label to HCM and VCM numbers
315        N = NHCM + NVCM;
316        Ticks = 1:ceil(N/10):N;
317        set(gca, 'XTick', Ticks);     
318        TickNumber = [CMData.HCMGoodDataIndex(:)' CMData.VCMGoodDataIndex(:)'];
319        TickNumber = TickNumber(Ticks);
320        set(gca, 'XTickLabel', num2cell(TickNumber));         
321        xlabel('Horizontal and Vertical Corrector Magnets');
322       
323        if strcmp(lower(PlotString),'model')
324            title('Model Response Matrix');
325            if length(Elements) == 1
326                ylabel(sprintf('Vertical BPM(%d) [mm]', BPMData.VBPMGoodDataIndex(Elements)));
327            else
328                ylabel('Vertical BPM [mm]');
329            end
330        elseif strcmp(lower(PlotString),'measured')
331            title('Measured Response Matrix');
332            if length(Elements) == 1
333                ylabel(sprintf('Vertical BPM(%d) [mm]', BPMData.VBPMGoodDataIndex(Elements)));
334            else
335                ylabel('Vertical BPM [mm]');
336            end
337        elseif strcmp(lower(PlotString),'model-measured')
338            title('Model - Measured Response Matrix');
339            if length(Elements) == 1
340                ylabel(sprintf('Vertical Error BPM(%d) [mm]', BPMData.VBPMGoodDataIndex(Elements)));
341            else
342                ylabel('Vertical BPM Error [mm]');
343            end
344        elseif strcmp(lower(PlotString),'model-measured+energyshifts')
345            title('Model-Measured+EnergyShifts Response Matrix');
346            if length(Elements) == 1
347                ylabel(sprintf('Vertical Error BPM(%d) [mm]', BPMData.VBPMGoodDataIndex(Elements)));
348            else
349                ylabel('Vertical BPM Error [mm]');
350            end
351        end
352       
353        % Only put a legend up to size 10
354        if length(Elements)>1 & length(Elements)<=10
355            for i = 1:length(Elements)
356                legendstr{i} = sprintf('BPM(%d)',BPMData.VBPMGoodDataIndex(Elements(i)));
357            end
358            legend(legendstr,0);
359        end
360        axis tight
361    end
362end
363
364
365if strcmp(PlaneString,'HCM')   
366    CMData.HCMGoodDataIndex = CMData.HCMGoodDataIndex(:)';
367    CMData.HCMIndex = CMData.HCMIndex(:);
368    if ~exist('ElementsInput','var')
369        ElementsInput = locoeditlist([CMData.HCMGoodDataIndex(:) CMData.HCMIndex(CMData.HCMGoodDataIndex)], PlaneString, zeros(length(CMData.HCMGoodDataIndex),1));
370        if isempty(ElementsInput)
371            return
372        end
373        ElementsInput = ElementsInput(:,1)';
374    end
375    if isempty(ElementsInput)
376        ElementsInput = locoeditlist([CMData.HCMGoodDataIndex(:) CMData.HCMIndex(CMData.HCMGoodDataIndex)], PlaneString, zeros(length(CMData.HCMGoodDataIndex),1));
377        if isempty(ElementsInput)
378            return
379        end
380        ElementsInput = ElementsInput(:,1)';
381    end
382
383    % Convert to elements in the HCMGoodDataIndex vector
384    [tmp1, tmp2, Elements] = intersect(ElementsInput, CMData.HCMGoodDataIndex);
385   
386    if ~isempty(Elements)
387        plot(data(:,Elements));
388
389        % Change label to BPM numbers
390        N = NBPM;
391        Ticks = 1:ceil(N/10):N;
392        set(gca, 'XTick', Ticks);     
393        TickNumber = [BPMData.HBPMGoodDataIndex(:)' BPMData.VBPMGoodDataIndex(:)'];
394        TickNumber = TickNumber(Ticks);
395        set(gca, 'XTickLabel', num2cell(TickNumber));         
396        xlabel('Horizontal and Vertical BPM Numbers');
397               
398        if strcmp(lower(PlotString),'model')
399            title('Model Response Matrix');
400            if length(ElementsInput) == 1
401                ylabel(sprintf('HCM(%d) [mm]', ElementsInput));
402            else
403                ylabel('HCM [mm]');
404            end
405        elseif strcmp(lower(PlotString),'measured')
406            title('Measured Response Matrix');
407            if length(ElementsInput) == 1
408                ylabel(sprintf('HCM(%d) [mm]', ElementsInput));
409            else
410                ylabel('HCM [mm]');
411            end
412        elseif strcmp(lower(PlotString),'model-measured')
413            title('Model - Measured Response Matrix');
414            if length(ElementsInput) == 1
415                ylabel(sprintf('HCM(%d) Error [mm]', ElementsInput));
416            else
417                ylabel('HCM Error [mm]');
418            end
419        elseif strcmp(lower(PlotString),'model-measured+energyshifts')
420            title('Model-Measured+EnergyShifts Response Matrix');
421            if length(ElementsInput) == 1
422                ylabel(sprintf('HCM(%d) Error [mm]', ElementsInput));
423            else
424                ylabel('HCM Error [mm]');
425            end
426        end
427       
428        % Only put a legend up to size 10
429        if length(ElementsInput)>1 & length(ElementsInput)<=10
430            for i = 1:length(ElementsInput)
431                legendstr{i} = sprintf('HCM(%d)',ElementsInput(i));
432            end
433            legend(legendstr,0);
434        end
435        axis tight
436    end
437end
438
439
440
441if strcmp(PlaneString,'VCM')
442    CMData.VCMGoodDataIndex = CMData.VCMGoodDataIndex(:)';
443    CMData.VCMIndex = CMData.VCMIndex(:);
444    if ~exist('ElementsInput','var')
445        ElementsInput = locoeditlist([CMData.VCMGoodDataIndex(:) CMData.VCMIndex(CMData.VCMGoodDataIndex)], PlaneString, zeros(length(CMData.VCMGoodDataIndex),1));
446        if isempty(ElementsInput)
447            return
448        end
449        ElementsInput = ElementsInput(:,1)';
450    end
451    if isempty(ElementsInput)
452        ElementsInput = locoeditlist([CMData.VCMGoodDataIndex(:) CMData.VCMIndex(CMData.VCMGoodDataIndex)], PlaneString, zeros(length(CMData.VCMGoodDataIndex),1));
453        if isempty(ElementsInput)
454            return
455        end
456        ElementsInput = ElementsInput(:,1)';
457    end
458           
459    % Convert to elements in the VCMGoodDataIndex vector
460    [tmp1, tmp2, Elements] = intersect(ElementsInput, CMData.VCMGoodDataIndex);
461   
462    if ~isempty(Elements)
463        plot(data(:,NHCM+Elements));
464       
465        % Change label to BPM numbers
466        N = NBPM;
467        Ticks = 1:ceil(N/10):N;
468        set(gca, 'XTick', Ticks);     
469        TickNumber = [BPMData.HBPMGoodDataIndex(:)' BPMData.VBPMGoodDataIndex(:)'];
470        TickNumber = TickNumber(Ticks);
471        set(gca, 'XTickLabel', num2cell(TickNumber));         
472        xlabel('Horizontal and Vertical BPM Numbers');
473               
474        if strcmp(lower(PlotString),'model')
475            title('Model Response Matrix');
476            if length(Elements) == 1
477                ylabel(sprintf('VCM(%d) [mm]', CMData.VCMGoodDataIndex(Elements)));
478            else
479                ylabel('VCM [mm]');
480            end
481        elseif strcmp(lower(PlotString),'measured')
482            title('Measured Response Matrix');
483            if length(Elements) == 1
484                ylabel(sprintf('VCM(%d) [mm]', CMData.VCMGoodDataIndex(Elements)));
485            else
486                ylabel('VCM [mm]');
487            end
488        elseif strcmp(lower(PlotString),'model-measured')
489            title('Model - Measured Response Matrix');
490            if length(Elements) == 1
491                ylabel(sprintf('VCM(%d) Error [mm]', CMData.VCMGoodDataIndex(Elements)));
492            else
493                ylabel('VCM Error [mm]');
494            end
495        elseif strcmp(lower(PlotString),'model-measured+energyshifts')
496            title('Model-Measured+EnergyShifts Response Matrix');
497            if length(Elements) == 1
498                ylabel(sprintf('VCM(%d) Error [mm]', CMData.VCMGoodDataIndex(Elements)));
499            else
500                ylabel('VCM Error [mm]');
501            end
502        end
503       
504        % Only put a legend up to size 10
505        if length(Elements)>1 & length(Elements)<=10
506            for i = 1:length(Elements)
507                legendstr{i} = sprintf('VCM(%d)',CMData.VCMGoodDataIndex(Elements(i)));
508            end
509            legend(legendstr,0);
510        end
511        axis tight
512    end
513end
514
515if ~exist('ElementsInput','var') & nargout > 0
516    ElementsInput = [];
517end
Note: See TracBrowser for help on using the repository browser.