FitLocalErrors.m 6.89 KB
Newer Older
1
function [knl, ksl, serr] = FitLocalErrors(....
2
3
    ReferenceRM,...
    PerturbedRM,...
4
5
    errorlocation,...
    neig)
6
7
8
9
%FITLOCALERRORS fit the best error to reproduced the difference between 
% Reference and Perturebed RM at errorlocation s positions
% 
% INPUT:
10
11
12
13
% ReferenceRM   : folder where the RM data is loaded or RM itself
% PerturbedRM   : folder where the RM data is loaded or RM itself
% errorlocation : list of s positions of error sources (
% neig          : # of eigenvectors for svd (default: length(errorlocation))
14
15
%
% OUTPUT:
16
17
% integrated normal (knl) and skew (ksl) quadrupole gradients (1/m) and
% serr ordered locations of errors
18
19
20
%
% EXAMPLE:
% 
21
% >> [knl, ksl] = FitLocalErrors('/machfs/MDT/2018/Sep26/respAC14','/machfs/MDT/2018/Sep26/respAC22',[286.34 584.438]);
22
%
23
24
25
% WARNING: 
% errorlocation should not be inside dipoles! an
% 
26
%see also: qemextract atsbreak
27
28
29
if nargin<4
    neig=length(errorlocation); 
end
30
31
32
33
34
35
36
37
38
39
40

% load reference (normal and skew)
[qemb, qemres, semres] = qemextract(ReferenceRM);

% load perturbed (normal and skew)
[~, qemresPert, semresPert] = qemextract(PerturbedRM);

% introduce generic straight multipole element with zero strenghts at given
% s location
r = qemb(2).at;

41
42
43
44
45
46
C = findspos(r,length(r)+1);

if any(errorlocation>C) | any(errorlocation<0)
    error(['Some locations are below 0 or above length of lattice (' num2str(C) ')']);
end

47
48
49
50
[rerr,errpos]=atsbreak(r,errorlocation);
Lthin = 1e-6;

rerr = atsetfieldvalues(rerr,errpos,'PassMethod', 'StrMPoleSymplectic4Pass');%'ThinMPolePass';
51
rerr = atsetfieldvalues(rerr,errpos,'Class', 'Quadrupole');% NECESSARY for FIT function
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
rerr = atsetfieldvalues(rerr,errpos,'PolynomB',{1,1}, 0);
rerr = atsetfieldvalues(rerr,errpos,'PolynomA',{1,1},0);
rerr = atsetfieldvalues(rerr,errpos,'PolynomB',{1,2}, 0);
rerr = atsetfieldvalues(rerr,errpos,'PolynomA',{1,2},0);
rerr = atsetfieldvalues(rerr,errpos,'MaxOrder', 1);
rerr = atsetfieldvalues(rerr,errpos,'Length', Lthin);
rerr = atsetfieldvalues(rerr,errpos,'NumIntSteps', 1);

% redefine fields of qemres and semres for fit of rm difference with local error
qemres.qpidx = errpos; % index of location to fit
qemres.qpfit = true(size(errpos));
qemres.skewidx = errpos;
semres.skewkeep = true(size(errpos));

qemb(2).at = rerr;
qemres.at = rerr;
% recompute indexes (new elements introduced)

thinsext=atgetcells(qemres.at(:),'Class',@isthinsext,'ThinMultipole');
thicksext=atgetcells(qemres.at(:),'Class','Sextupole');
qemres.sextidx=reshape(find(thinsext|thicksext),1,[]);
qemres.steerhidx=qemres.sextidx(selcor(8));
qemres.steervidx=qemres.sextidx(selcor(9));
qemres.bpmidx=findcells(qemres.at(:),'Class','Monitor');
dipfold=sr.fold(findcells(qemres.at(:),'Class','Bend'));
qemres.dipidx=reshape(dipfold([1 4],:),1,64);
qemres.ql=atgetfieldvalues(qemres.at(qemres.qpidx),'Length'); % this vector is used by semquadrespAnalytic
qemres.ql(qemres.ql==0) = 1; % integrated strengths;

% get analytic rm derivatives
% hw=waitbar(0,'Fitting quadrupole strengths (analytic)...');
% dispfunc=@(i,itot) waitbar(i/itot,hw);
dispfunc=@(a,b)disp(['Fitting quadrupole strengths (analytic): ' num2str(a) '/' num2str(b)]);
beamdyn's avatar
beamdyn committed
85
okbpm(1:224) = true; okbpm=okbpm';% added by Andrea
86
87
88
89
90
91
92
93
94
95
96
97
98
okbpm(qemres.wrongbpms) = false;

shidx=qemres.steerhidx(qemres.hlist);
svidx=qemres.steervidx(qemres.vlist);
nhst=length(shidx);
nvst=length(svidx);

[~,o0]=findsyncorbit(qemb(2).at,qemres.ct,qemres.bpmidx);
dpp=o0(5);
[~,fractunes]=atlinopt(qemb(2).at,dpp,[],o0);

dresp=qemderivAnalytic(qemb(2).at,dpp,@setk,0.0001,qemres.qpidx,qemres.bpmidx,shidx,svidx,dispfunc);
drespskew=semquadrespAnalytic(qemb(2),qemres,semres);
beamdyn's avatar
beamdyn committed
99
100
% save -ascii pippo.txt dresp;
% return; 
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

% load RM
wrongbpms=qemres.wrongbpms; % find(bpmm.All_Status.read~=0)
[semres.resph,semres.respv]=semloadresp(qemres,semres,wrongbpms);
[qemres.resph,qemres.respv,qemres.frespx,semres.frespz]=qemcheckresp(qemres,semres,wrongbpms);
semres.respiax=load_sriaxresp(qemres.datadir,qemres.opticsdir,'h2v',semres.hlist);
qemres.respiax=load_sriaxresp(qemres.datadir,qemres.opticsdir,'v',qemres.vlist);

wrongbpms=qemresPert.wrongbpms; % find(bpmm.All_Status.read~=0)
qemresPert.hlist = qemres.hlist;
qemresPert.vlist = qemres.vlist;
semresPert.hlist = semres.hlist;
semresPert.vlist = semres.vlist;

[semresPert.resph,semresPert.respv]=semloadresp(qemresPert,semresPert,wrongbpms);
[qemresPert.resph,qemresPert.respv,qemresPert.frespx,semresPert.frespz]=qemcheckresp(qemresPert,semresPert,wrongbpms);
semresPert.respiax=load_sriaxresp(qemresPert.datadir,qemresPert.opticsdir,'h2v',semresPert.hlist);
qemresPert.respiax=load_sriaxresp(qemresPert.datadir,qemresPert.opticsdir,'v',qemresPert.vlist);

% get perturbed RM
[rh,rh2v,rv2h,rv,frh,frv]=qemdecode(qemb(2).at,qemresPert.ct,qemresPert,...
    qemresPert.resph,semresPert.respv,semresPert.resph,qemresPert.respv,...
    qemresPert.frespx,semresPert.frespz,...
    qemres.khrot,qemres.khgain,qemres.kvrot,qemres.kvgain,...
    qemres.brot,qemres.bhgain,qemres.bvgain);

% get reference RM
[rh0,rh2v0,rv2h0,rv0,frh0,frv0]=qemdecode(qemb(2).at,qemres.ct,qemres,...
    qemres.resph,semres.respv,semres.resph,qemres.respv,...
    qemres.frespx,semres.frespz,...
    qemres.khrot,qemres.khgain,qemres.kvrot,qemres.kvgain,...
    qemres.brot,qemres.bhgain,qemres.bvgain);

resp0=[rh0(:);rv0(:);frh0;qemres.fractunes'];
resp=[rh(:);rv(:);qemres.bhscale*frh;qemresPert.fractunes'];
resps0=[rh2v0(:);rv2h0(:);frv0];
resps=[rh2v(:);rv2h(:);semres.bvscale*frv];

% return fitted errors and residual of fit before and after
mode = struct;
okfit=true(size(qemres.qpidx));
if ~isfield(mode,'vnorm')
    vv=sqrt([0.2;1.1;1.95;2;2;1.05;0.64;0.9]);
    vv=sr.unfold(repmat(vv,1,32));
    mode.vnorm=vv(okfit)';
end

beamdyn's avatar
beamdyn committed
148
mode.nsets=4; 
149
mode.neigs=neig;
beamdyn's avatar
beamdyn committed
150
151
152
mode.tuneweight=0; 
mode.dispweight=0; 

153
154
kn = qemerrfit(nhst,nvst,resp-resp0,dresp,mode,okbpm);

beamdyn's avatar
beamdyn committed
155
mode.vnorm=sqrt([0.77;2;0.9;1.15;1.15;0.9;2;1]);
156
157
158
159
160
161
162

resps(isnan(resps))=0;
resps0(isnan(resps0))=0;

ks = semerrfit(length(semres.hlist),length(semres.vlist),resps-resps0,...
    drespskew,mode);

163
164
knl = kn *Lthin;
ksl = ks *Lthin;
165
166

% display fit residual
beamdyn's avatar
beamdyn committed
167
rm_in=(resp-resp0); rm_in(end-length(qemres.bpmidx)-2:end)=[];
168
dresp(end-length(qemres.bpmidx)-2:end,:)=[];
beamdyn's avatar
beamdyn committed
169
170
disp(['residual on- diagonal blocks: ' num2str(std2(rm_in)) '->'...
    num2str(std2(rm_in - dresp*kn))])
171
172
173
disp(['residual off-diagonal blocks: ' num2str(std2(resps-resps0)) '->'...
    num2str(std2(resps-resps0 - drespskew*ks))])

174
175
176
177
178
179
180
181
% display results
serr = findspos(rerr,sort(errpos)); % fitted gradient are sorted in lattice order
for il=1:length(serr)
    disp(['@ ' num2str(serr(il),'%6.2f') ' m : knl= '...
        num2str(knl(il),'%9.2e') ' 1/m, ksl= '...
        num2str(ksl(il),'%9.2e') ' 1/m']);
end

182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
end


% ancillary functions
function ok=isthinsext(elem,fieldname,fieldvalue)
if isfield(elem,fieldname) && strcmp(elem.(fieldname),fieldvalue)
    pb=find(elem.PolynomB ~= 0,1);
    ok=~isempty(pb) && (pb == 3);
else
    ok=false;
end
end

function elem=setk(elem,dk)
strength=elem.PolynomB(2)+dk;
elem.K=strength;
elem.PolynomB(2)=strength;
end