fit difference between two RM using normal and skew quadrupole gradients at given s locations

qem/FitLocalErrors/FitLocalErrors.m

+function [knl, ksl] = FitLocalErrors(....
+    ReferenceRM,...
+    PerturbedRM,...
+    errorlocation)
+%FITLOCALERRORS fit the best error to reproduced the difference between 
+% Reference and Perturebed RM at errorlocation s positions
+% 
+% INPUT:
+% 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
+%
+% OUTPUT:
+% integrated normal (knl) and skew (ksl) quadrupole gradients (1/m)
+%
+% EXAMPLE:
+% 
+% >> [knl, ksl] = FitLocalErrors('/machfs/MDT/2018/Sep17/respAC1','/machfs/MDT/2018/Mar25/resp1',[345.5,824.2])
+%
+%see also: qemextract atsbreak
+
+% 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;
+
+[rerr,errpos]=atsbreak(r,errorlocation);
+Lthin = 1e-6;
+
+rerr = atsetfieldvalues(rerr,errpos,'PassMethod', 'StrMPoleSymplectic4Pass');%'ThinMPolePass';
+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)]);
+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);
+
+% 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
+
+kn = qemerrfit(nhst,nvst,resp-resp0,dresp,mode,okbpm);
+
+mode=struct('vnorm',sqrt([0.77;2;0.9;1.15;1.15;0.9;2;1]));
+
+resps(isnan(resps))=0;
+resps0(isnan(resps0))=0;
+
+ks = semerrfit(length(semres.hlist),length(semres.vlist),resps-resps0,...
+    drespskew,mode);
+
+knl = kn .*Lthin;
+ksl = ks .*Lthin;
+
+% display fit residual
+disp(['residual on- diagonal blocks: ' num2str(std2(resp-resp0)) '->'...
+    num2str(std2(resp-resp0 - dresp*kn))])
+disp(['residual off-diagonal blocks: ' num2str(std2(resps-resps0)) '->'...
+    num2str(std2(resps-resps0 - drespskew*ks))])
+
+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
+