source: MML/trunk/mml/meascmhysteresis.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: 10.9 KB
Line 
1function meascmhysteresis(BPMList, CMFamily, CMList, MaxChange, NSteps)
2%MEASCMHYSTERESIS - Measure corrector magnet hysteresis
3%
4%  meascmhysteresis(BPMList, MaxChange, NSteps)
5%
6%  INPUTS
7%  1. BPMList   - BPM list
8%  2. CMFamily  - Corrector family
9%  3. CMList    - Corrector device list
10%  4. MaxChange - Maximum vertical change from the starting point
11%  5. NSteps    - Number of step from starting point to the maximum
12%
13%  NOTES
14%  1. This function looks for a bump coefficient file corresponding to BPM list. 
15%     For instance, VerticalBumpCoef-7_8-6_1 is the vertical file for [7 8;6 1].
16%     If you don't want to use this file, delete it and a new new will be created.
17%     All files are stored in <DataRoot><aperturescan>
18
19%
20%  Written by Gregory J. Portmann
21%  Adapted by Laurent S. Nadolski
22
23BPMxFamily = gethbpmfamily;
24BPMyFamily = getvbpmfamily;
25HCMFamily  = gethcmfamily;
26VCMFamily  = getvcmfamily;
27
28
29%%%%%%%%%%%%%%%%%%
30% Input checking %
31%%%%%%%%%%%%%%%%%%
32if nargin < 1
33    BPMList = [7 6; 8 1];
34end
35if nargin < 2
36    CMFamily = VCMFamily;
37end
38if nargin < 3
39    if size(BPMList,1) == 1
40        CMList = [4 4];
41    else
42        CMList = [-4 -3 -2 -1 1 2 3 4];
43    end
44end
45if nargin < 4
46    MaxChange = 1.5;
47end
48if nargin < 5
49    NSteps = 5;     % Steps from 0 to MaxChange
50end
51
52
53% Corrector starting points
54HCM0 = getsp(HCMFamily);
55VCM0 = getsp(VCMFamily);
56
57
58% Offset
59Xoffset = getoffset(BPMxFamily, BPMList);
60Yoffset = getoffset(BPMyFamily, BPMList);
61
62% Gain
63Xgain = getphysdata(BPMxFamily, 'Gain', BPMList);
64Ygain = getphysdata(BPMyFamily, 'Gain', BPMList);
65
66
67BPMxIndex = findrowindex(BPMList, getlist(BPMxFamily));
68BPMyIndex = findrowindex(BPMList, getlist(BPMyFamily));
69
70BPMxList = getlist(BPMxFamily);
71BPMyList = getlist(BPMyFamily);
72
73FileNumberString =  sprintf('-%d_%d', BPMList');
74DirectoryName = [fullfile(getfamilydata('Directory','DataRoot'),'aperturescan'),filesep];
75
76
77%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78% Compute corrector coefficients %
79%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80BPMWeight = 40;
81SVDIndex = 1e-3;
82
83if strcmpi(CMFamily, HCMFamily)
84    HorizontalBumpFile = sprintf('%sHorizontalBumpCoef%s', DirectoryName, FileNumberString);
85    if size(BPMList,1) == 1 && isempty(CMList)
86        % For single BPM just use one corrector
87        % Pick the corrector based on the response matrix
88        R = getbpmresp('Struct');
89        [i, iNotFound] = findrowindex(BPMList, R(1,1).Monitor.DeviceList);
90        m = R(1,1).Data(i,:);
91        [MaxValue, j] = max(abs(m));
92        HCMDeviceList = R(1,1).Actuator.DeviceList(j,:);
93        %DeltaHCM = (1/m(j)) ./ Xgain;  % Amps/mm  (at the BPM BPMList)
94        DeltaHCM = (1/max(abs(m))) ./ Xgain;  % Amps/mm  (maximum at many BPM)
95       
96    elseif size(CMList,1) == 1 && size(CMList,2) <= 2
97        % For single corrector
98        R = getbpmresp('Struct');
99        [i, iNotFound] = findrowindex(BPMList, R(1,1).Monitor.DeviceList);
100        [j, jNotFound] = findrowindex(CMList,  R(1,1).Actuator.DeviceList);
101        m = R(1,1).Data(i,j);
102        DeltaHCM = (1/m) ./ Xgain;  % Amps/mm
103       
104    elseif exist([HorizontalBumpFile,'.mat'], 'file')
105        fprintf('   Loading horizontal corrector magnet bump coefficients from %s\n', ...
106            HorizontalBumpFile);
107        load(HorizontalBumpFile);
108       
109    else
110        % Set the hysteresis
111        fprintf('   Finding horizontal corrector magnet bump coefficients\n');
112        DeltaX = [.2;.2] ./ Xgain;
113        x0 = getam(BPMxFamily, BPMList);
114        %[HOCS, RF, HOCS0] = setorbitbump(BPMxFamily, BPMList, DeltaX, HCMFamily, CMList, 1, SVDIndex, BPMWeight, 'Inc', 'NoDisplay');
115        HOCS0 = setorbitbump(BPMxFamily, BPMList, DeltaX, HCMFamily, CMList, 1, ...
116            BPMWeight, 'Inc', 'NoDisplay');
117        x1 = getam(BPMxFamily, BPMList);
118       
119        % Get a clean bump
120        figure(1);
121        clf reset
122        x2 = getam(BPMxFamily, BPMList);
123        %[HOCS, RF, HOCS0] = setorbitbump(BPMxFamily, BPMList, DeltaX, HCMFamily, CMList, 1, SVDIndex, BPMWeight, 'Inc', 'NoDisplay');
124        HOCS = setorbitbump(BPMxFamily, BPMList, DeltaX, HCMFamily, CMList, 1, ...
125            BPMWeight, 'Inc', 'NoDisplay');
126        x3 = getam(BPMxFamily, BPMList);
127        drawnow;
128        DeltaHCM = (HOCS.CM.Data - HOCS0.CM.Data)/.2;  % amps/mm  (real units)
129        BPMDeviceList = HOCS.BPM.DeviceList;
130        HCMDeviceList = HOCS.CM.DeviceList;
131        setsp(HCMFamily, HCM0);
132        save(HorizontalBumpFile, 'DeltaHCM', 'HCMDeviceList', 'BPMDeviceList', 'HOCS', 'HOCS0');
133    end
134   
135    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136    % Linearity and hystersis check %
137    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138   
139    % Starting corrector values
140    hcm0 = getsp(HCMFamily, HCMDeviceList);
141   
142    % Set hysteresis in an upward direction
143    setsp(HCMFamily, hcm0 - MaxChange * DeltaHCM, HCMDeviceList, -2);
144    setsp(HCMFamily, hcm0, HCMDeviceList);
145   
146    % Build loop
147    Loop = [linspace(0,MaxChange,NSteps+1) linspace(MaxChange-MaxChange/NSteps,-MaxChange,2*NSteps) linspace(-MaxChange+MaxChange/NSteps,0,NSteps)];
148    LoopTotal = [Loop Loop(2:end)];
149   
150    i = 0;
151    for y = LoopTotal
152        i = i + 1;
153        setsp(HCMFamily, hcm0 + y * DeltaHCM, HCMDeviceList, -2);
154        DCCT(i) = getdcct;
155        HCMsp(:,i) = getsp(HCMFamily, HCMDeviceList);
156        HCMam(:,i) = getam(HCMFamily, HCMDeviceList);
157        %BPMx(:,i) = raw2real(BPMxFamily, getx);
158        %BPMy(:,i) = raw2real(BPMyFamily, getz);
159        BPMx(:,i) = getx;
160        BPMy(:,i) = getz;
161    end
162   
163    % Reset correctors
164    setsp(HCMFamily, HCM0);
165   
166elseif strcmpi(CMFamily, VCMFamily)
167    VerticalBumpFile = sprintf('VerticalBumpCoef%s', FileNumberString);
168    if size(BPMList,1) == 1 && isempty(CMList)
169        % For single BPM just use one corrector
170        % Pick the corrector based on the response matrix
171        R = getbpmresp('Struct');
172        [i, iNotFound] = findrowindex(BPMList, R(2,2).Monitor.DeviceList);
173        m = R(2,2).Data(i,:);
174        [MaxValue, j] = max(abs(m));
175        VCMDeviceList = R(2,2).Actuator.DeviceList(j,:);
176        %DeltaVCM = (1/m(j)) ./ Ygain;  % Amps/mm  (at the BPM BPMList)
177        DeltaVCM = (1/max(abs(m))) ./ Ygain;  % Amps/mm  (maximum at many BPM)
178       
179    elseif size(CMList,1) == 1 && size(CMList,2) <= 2
180        % For single corrector
181        R = getbpmresp('Struct');
182        [i, iNotFound] = findrowindex(BPMList, R(2,2).Monitor.DeviceList);
183        [j, jNotFound] = findrowindex(CMList,  R(2,2).Actuator.DeviceList);
184        m = R(2,2).Data(i,j);
185        DeltaVCM = (1/m) ./ Ygain;  % Amps/mm
186        VCMDeviceList = CMList;
187       
188    elseif exist([VerticalBumpFile,'.mat'], 'file')
189        fprintf('   Loading vertical corrector magnet bump coefficients from %s\n', VerticalBumpFile);
190        load(VerticalBumpFile);
191    else
192        % Set the hysteresis
193        fprintf('   Finding vertical corrector magnet bump coefficients\n');
194        DeltaY = [.2;.2] ./ Ygain;
195        y0 = getam(BPMyFamily, BPMList);
196        %[VOCS, RF, VOCS0] = setorbitbump(BPMyFamily, BPMList, DeltaY, VCMFamily, CMList, 1, SVDIndex, BPMWeight, 'Inc', 'NoDisplay');
197        % in one iteration       
198        VOCS0 = setorbitbump(BPMyFamily, BPMList, DeltaY, VCMFamily, CMList, ...
199            1, BPMWeight, 'Inc', 'NoDisplay');
200        y1 = getam(BPMyFamily, BPMList);
201       
202        % Get a clean bump
203        %vcm0= getsp(VCMFamily, VOCS.CM.DeviceList);
204        figure(1);
205        clf reset
206        y2 = getam(BPMyFamily, BPMList);
207        VOCS = setorbitbump(BPMyFamily, BPMList, DeltaY, VCMFamily, CMList, ...
208            5,  BPMWeight, 'Inc', 'NoDisplay');
209%        [VOCS, RF, VOCS0] = setorbitbump(BPMyFamily, BPMList, DeltaY, VCMFamily, CMList, 5, SVDIndex, BPMWeight, 'Inc', 'NoDisplay');
210        y3 = getam(BPMyFamily, BPMList);
211        drawnow;
212        %vcm1= getsp(VCMFamily, VOCS.CM.DeviceList);
213        DeltaVCM = (VOCS.CM.Data - VOCS0.CM.Data) / .2;  % amps/mm  (real units)
214        BPMDeviceList = VOCS.BPM.DeviceList;
215        VCMDeviceList = VOCS.CM.DeviceList;
216        setsp(VCMFamily, VCM0);
217        save(VerticalBumpFile, 'DeltaVCM', 'VCMDeviceList', 'BPMDeviceList', 'VOCS', 'VOCS0');
218    end
219   
220   
221    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222    % Linearity and hystersis check %
223    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
224   
225    % Starting corrector values
226    vcm0 = getsp(VCMFamily, VCMDeviceList);
227   
228    % Set hysteresis in an upward direction
229    %setsp(VCMFamily, vcm0 - MaxChange * DeltaVCM, VCMDeviceList, -2);
230    %setsp(VCMFamily, vcm0, VCMDeviceList);
231   
232    % Build loop
233    Loop = [linspace(0,MaxChange,NSteps+1) ...
234        linspace(MaxChange-MaxChange/NSteps,-MaxChange,2*NSteps) ...
235        linspace(-MaxChange+MaxChange/NSteps,0,NSteps)];
236    LoopTotal = [Loop Loop(2:end)];
237   
238    i = 0;
239    for y = LoopTotal
240        i = i + 1;
241        setsp(VCMFamily, vcm0 + y * DeltaVCM, VCMDeviceList, -2);
242        DCCT(i) = getdcct;
243        VCMsp(:,i) = getsp(VCMFamily, VCMDeviceList);
244        VCMam(:,i) = getam(VCMFamily, VCMDeviceList);
245        %BPMx(:,i) = raw2real(BPMxFamily, getx);
246        %BPMy(:,i) = raw2real(BPMyFamily, getz);
247        BPMx(:,i) = getx;
248        BPMy(:,i) = getz;
249       
250    end
251   
252    % Reset correctors
253    setsp(VCMFamily, VCM0);
254   
255else
256    error('Corrector magnet family unknown');
257end
258
259
260% Save data
261HysteresisFilename = sprintf('Hysteresis%s', FileNumberString);
262HysteresisFilename = appendtimestamp(HysteresisFilename);
263save(HysteresisFilename);
264fprintf('   Hysteresis data saved to %s.mat\n', HysteresisFilename);
265
266
267% Plot results
268clf reset
269LineType = '.-b';
270if size(BPMList,1) == 1   
271   
272    subplot(2,1,1);
273    plot(LoopTotal, BPMy(BPMyIndex(1),:),LineType);
274    xlabel(sprintf('Goal Orbit at %s(%d,%d) [mm]', BPMyFamily, BPMList(1,:)));
275    ylabel(sprintf('%s(%d,%d) [mm]', BPMyFamily, BPMList(1,:)));
276   
277    subplot(2,1,2);
278    plot(LoopTotal, BPMy(BPMyIndex(1),:)-LoopTotal,LineType);
279    xlabel(sprintf('Goal Orbit at %s(%d,%d) [mm]', BPMyFamily, BPMList(1,:)));
280    ylabel(sprintf('%s(%d,%d) Error [mm]', BPMyFamily, BPMList(1,:)));
281    title('Hysteresis Loops');
282   
283else
284   
285    subplot(2,2,1);
286    plot(LoopTotal, BPMy(BPMyIndex(1),:),LineType);
287    xlabel(sprintf('Goal Orbit at %s(%d,%d) [mm]', BPMyFamily, BPMList(1,:)));
288    ylabel(sprintf('%s(%d,%d) [mm]', BPMyFamily, BPMList(1,:)));
289    grid on;
290
291    subplot(2,2,3);
292    plot(LoopTotal, BPMy(BPMyIndex(1),:)-LoopTotal,LineType);
293    xlabel(sprintf('Goal Orbit at %s(%d,%d) [mm]', BPMyFamily, BPMList(1,:)));
294    ylabel(sprintf('%s(%d,%d) Error [mm]', BPMyFamily, BPMList(1,:)));
295    grid on;
296   
297    subplot(2,2,2);
298    plot(LoopTotal, BPMy(BPMyIndex(2),:),LineType);
299    xlabel(sprintf('Goal Orbit at %s(%d,%d) [mm]', BPMyFamily, BPMList(2,:)));
300    ylabel(sprintf('%s(%d,%d) [mm]', BPMyFamily, BPMList(2,:)));
301    grid on;
302   
303    subplot(2,2,4);
304    plot(LoopTotal, BPMy(BPMyIndex(2),:)-LoopTotal,LineType);
305    xlabel(sprintf('Goal Orbit at %s(%d,%d) [mm]', BPMyFamily, BPMList(2,:)));
306    ylabel(sprintf('%s(%d,%d) Error [mm]', BPMyFamily, BPMList(2,:)));
307    grid on;
308   
309    h = addlabel(.5, 1, 'Hysteresis Loops');
310    set(h, 'Fontsize', 12);
311    orient tall
312   
313end
314
Note: See TracBrowser for help on using the repository browser.