FitLocalErrors: additional input n eigvectors, error for bad positions,...

FitLocalErrors: additional input n eigvectors, error for bad positions, display of results, additional output of sorted positions

qem/FitLocalErrors/FitLocalErrors.m

-function [knl, ksl] = FitLocalErrors(....
+function [knl, ksl, serr] = FitLocalErrors(....
     ReferenceRM,...
     PerturbedRM,...
-    errorlocation)
+    errorlocation,...
+    neig)
 %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
+% 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))
 %
 % OUTPUT:
-% integrated normal (knl) and skew (ksl) quadrupole gradients (1/m)
+% integrated normal (knl) and skew (ksl) quadrupole gradients (1/m) and
+% serr ordered locations of errors
 %
 % EXAMPLE:
 % 
-% >> [knl, ksl] = FitLocalErrors('/machfs/MDT/2018/Sep17/respAC1','/machfs/MDT/2018/Mar25/resp1',[345.5,824.2])
+% >> [knl, ksl] = FitLocalErrors('/machfs/MDT/2018/Sep26/respAC14','/machfs/MDT/2018/Sep26/respAC22',[286.34 584.438]);
 %
+% WARNING: 
+% errorlocation should not be inside dipoles! an
+% 
 %see also: qemextract atsbreak
+if nargin<4
+    neig=length(errorlocation); 
+end
 
 % load reference (normal and skew)
 [qemb, qemres, semres] = qemextract(ReferenceRM);
@@ -29,6 +38,12 @@ function [knl, ksl] = FitLocalErrors(....
 % s location
 r = qemb(2).at;
 
+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
+
 [rerr,errpos]=atsbreak(r,errorlocation);
 Lthin = 1e-6;
 
@@ -131,7 +146,7 @@ if ~isfield(mode,'vnorm')
 end
 
 mode.nsets=4; 
-mode.neigs=length(errorlocation); 
+mode.neigs=neig;
 mode.tuneweight=0; 
 mode.dispweight=0; 
 
@@ -145,17 +160,25 @@ resps0(isnan(resps0))=0;
 ks = semerrfit(length(semres.hlist),length(semres.vlist),resps-resps0,...
     drespskew,mode);
 
-knl = kn ;
-ksl = ks ;
+knl = kn *Lthin;
+ksl = ks *Lthin;
 
 % display fit residual
 rm_in=(resp-resp0); rm_in(end-length(qemres.bpmidx)-2:end)=[];
-dresp(end-length(qemres.bpmidx)-2:end)=[];
+dresp(end-length(qemres.bpmidx)-2:end,:)=[];
 disp(['residual on- diagonal blocks: ' num2str(std2(rm_in)) '->'...
     num2str(std2(rm_in - dresp*kn))])
 disp(['residual off-diagonal blocks: ' num2str(std2(resps-resps0)) '->'...
     num2str(std2(resps-resps0 - drespskew*ks))])
 
+% 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
+
 end