atdispersionfreesteering

 function [...
    rcor,...           1) corrected lattice
    inCOD,...          2) initial COD (dpp is stored here)
    hs,vs...           3) required steerers strengths (total)
    ]=atdispersionfreesteering(...
     rerr,...          1) AT lattice to correct
     indBPM,...        2) Nbx1 bpm indexes
     indHCor,...       3) Nhx1 hor. cor indexes
     indVCor,...       4) Nvx1 ver. cor indexes
     inCOD,...         5) 6x1 initial COD guess
     neigSteerer,...   6) 2xNiter eigenvectors for correction H and V at
                          each iteration (default: [Nh/2 Nv/2])
     correctflags,...  7) correct [dpp mean0](default: [true true])
     scalefactor,...   8) scale factor to correction (default: 1.0)
     wdisp,...         9) weight dispersion*wdisp and orbit*(1-wdisp)
                          (default: 0.7)
     ModelRM,...       10) ModelRM.Orb(H/V)Cor = 6x1 cell of orbit response mat.
                           ModelRM.Orb(H/V)DPP = 6x1 array of orbit
                           ModelRM.Disp(H/V)Cor = 6x1 cell of dispersion response mat.
                           ModelRM.Disp(H/V)DPP = 6x1 array of dispersion
                           if [] compute RM (default: [])
     reforbit,...      11) 2xNbpm reference orbit to correct to (default 0*2xNb)
     refdispersion,... 12) 2xNbpm reference orbit to correct to (default 0*2xNb)
     steererlimit      13) 2x1 limit of steerers abs(steerer)<steererlimit
                           (default: [], no limits)
     printouttext      14) if 1 or true, display rms orbit
     )

 features impelemented:
 limit correctors strengths
 ddp correction
 sum of steerers = 0
 6D orbit with BPM errors
 initial coordinate
 correction to reference orbit refx refy
 use atsetfieldvalues, atgetcells


http://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.3.121001

see also: qemsvd_mod findorbit6Err getresponsematrices

0001 function [rcor,inCOD,hs,vs]=atdispersionfreesteering(...
0002     rerr,...
0003     indBPM,...
0004     indHCor,...
0005     indVCor,...
0006     inCOD,...
0007     neigSteerer,...
0008     correctflags,...
0009     scalefactor,...
0010     wdisp,...
0011     ModelRM,...
0012     reforbit,...
0013     refdispersion,...
0014     steererlimit,...
0015     printouttext)
0016 % function [...
0017 %    rcor,...           1) corrected lattice
0018 %    inCOD,...          2) initial COD (dpp is stored here)
0019 %    hs,vs...           3) required steerers strengths (total)
0020 %    ]=atdispersionfreesteering(...
0021 %     rerr,...          1) AT lattice to correct
0022 %     indBPM,...        2) Nbx1 bpm indexes
0023 %     indHCor,...       3) Nhx1 hor. cor indexes
0024 %     indVCor,...       4) Nvx1 ver. cor indexes
0025 %     inCOD,...         5) 6x1 initial COD guess
0026 %     neigSteerer,...   6) 2xNiter eigenvectors for correction H and V at
0027 %                          each iteration (default: [Nh/2 Nv/2])
0028 %     correctflags,...  7) correct [dpp mean0](default: [true true])
0029 %     scalefactor,...   8) scale factor to correction (default: 1.0)
0030 %     wdisp,...         9) weight dispersion*wdisp and orbit*(1-wdisp)
0031 %                          (default: 0.7)
0032 %     ModelRM,...       10) ModelRM.Orb(H/V)Cor = 6x1 cell of orbit response mat.
0033 %                           ModelRM.Orb(H/V)DPP = 6x1 array of orbit
0034 %                           ModelRM.Disp(H/V)Cor = 6x1 cell of dispersion response mat.
0035 %                           ModelRM.Disp(H/V)DPP = 6x1 array of dispersion
0036 %                           if [] compute RM (default: [])
0037 %     reforbit,...      11) 2xNbpm reference orbit to correct to (default 0*2xNb)
0038 %     refdispersion,... 12) 2xNbpm reference orbit to correct to (default 0*2xNb)
0039 %     steererlimit      13) 2x1 limit of steerers abs(steerer)<steererlimit
0040 %                           (default: [], no limits)
0041 %     printouttext      14) if 1 or true, display rms orbit
0042 %     )
0043 %
0044 % features impelemented:
0045 % limit correctors strengths
0046 % ddp correction
0047 % sum of steerers = 0
0048 % 6D orbit with BPM errors
0049 % initial coordinate
0050 % correction to reference orbit refx refy
0051 % use atsetfieldvalues, atgetcells
0052 %
0053 %
0054 %http://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.3.121001
0055 %
0056 %see also: qemsvd_mod findorbit6Err getresponsematrices
0057 
0058 
0059 
0060 % response matrix kicks
0061 kval=1e-5;
0062 delta=1e-3;
0063 
0064 % default arguments
0065 if nargin<14
0066     printouttext=true;
0067 end
0068 if nargin<13
0069     steererlimit=[];
0070 end
0071 
0072 if nargin<4
0073     if printouttext
0074         disp('get BPM and Correctors indexes'); end;
0075     indBPM=finc(atgetcells(rerr,'Class','Monitor'));
0076     indHCor=finc(atgetcells(rerr,'iscorH','H'));
0077     indVCor=finc(atgetcells(rerr,'iscorV','V'));
0078 end
0079 
0080 if nargin<5
0081     inCOD=[0 0 0 0 0 0]';
0082 end
0083 
0084 if nargin<6
0085     neigSteerer=[length(indHCor) length(indVCor)]/2;
0086 end
0087 
0088 if nargin<7
0089     correctflags=[true true];
0090 end
0091 
0092 if nargin<8
0093     if printouttext
0094         disp(' --- scale set to 1.0'); end;
0095     scalefactor=1.0;
0096 end
0097 
0098 if nargin<9
0099     if printouttext, disp(' --- alpha=0.7'); end;
0100     wdisp=0.8;
0101 end
0102 
0103 if nargin<10
0104     if printouttext, disp(' --- computing orbit Response matrix'); end;
0105     ModelRM=[];
0106 end
0107 
0108 if nargin<11
0109     if printouttext, disp(' --- reference orbit = 0'); end;
0110     reforbit=zeros(size(indBPM),2);
0111 end
0112 
0113 if nargin<12
0114     if printouttext, disp(' --- reference dispersion = 0 V, rerr disp H'); end;
0115     refdispersion=zeros(size(indBPM),2);
0116     [l,~,~]=atlinopt(rerr,0,indBPM);
0117     refdispersion(1,:)=arrayfun(@(a)a.Dispersion(1),l);
0118 end
0119 
0120 
0121 if scalefactor<0 || scalefactor>1
0122     if printouttext
0123         disp(' --- scale factor out of range. Set to 1.0'); end;
0124     scalefactor=1.0;
0125 end
0126 
0127 
0128 if correctflags(1) % dpp correction
0129     rmsel=[1 2 3 7 8 9];
0130 else
0131     rmsel=[1 2 7 8];
0132 end
0133 
0134 % load or compute response matrix
0135 if isempty(ModelRM)
0136     % get orbit RM
0137     if printouttext
0138         disp('get RM'); end;
0139     
0140     ModelRM=getresponsematrices(...
0141         rerr,...          %1 AT lattice
0142         indBPM,...      %2 bpm indexes in at lattice
0143         indHCor,...     %3 h cor indexes
0144         indVCor,...     %4 v cor indexes
0145         [],...     %5 skew cor indexes
0146         [],...     %6 quad cor indexes
0147         [],...
0148         inCOD,...       %7 initial coordinates
0149         rmsel...        %8 specifiy rm to be computed
0150         );
0151     
0152     if ~correctflags(1)
0153         
0154         ModelRM.OrbHDPP=[];
0155         ModelRM.OrbVDPP=[];
0156         ModelRM.DispHDPP=[];
0157         ModelRM.DispVDPP=[];
0158     end
0159     
0160 end
0161 
0162 % load RM computed by getresponsematrices
0163 
0164 ormH=ModelRM.OrbHCor;
0165 ormV=ModelRM.OrbVCor;
0166 drmH=ModelRM.DispHCor;
0167 drmV=ModelRM.DispVCor;
0168 % kval=ModelRM.kval;
0169 dppH=ModelRM.OrbHDPP;
0170 dppV=ModelRM.OrbVDPP;
0171 dppHd=ModelRM.DispHDPP;
0172 dppVd=ModelRM.DispVDPP;
0173 % delta=ModelRM.delta;
0174 alpha=mcf(rerr);
0175 indrfc=find(atgetcells(rerr,'Frequency'));
0176             
0177 
0178 % get initial orbit
0179 o=findorbit6Err(rerr,indBPM,inCOD);
0180 ox0=o(1,:);
0181 oy0=o(3,:);
0182 d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0183 dx0=d(1,:);
0184 dy0=d(3,:);
0185 
0186 %rerr0=rerr;
0187 
0188 % iterate correction
0189 Niter=size(neigSteerer,1);
0190 for iter=1:Niter
0191     
0192     if printouttext
0193         disp(['Disp. Free Steering iter ' num2str(iter,'%d, ') ...
0194             ' n-eig: ' num2str(neigSteerer(iter,:),'%d, ') ...
0195             ' alpha: ' num2str(wdisp,'%2.2f ')]);
0196     end
0197     
0198     % initial corrector strengths
0199     corh0=atgetfieldvalues(rerr,indHCor,'PolynomB',{1,1});
0200     corv0=atgetfieldvalues(rerr,indVCor,'PolynomA',{1,1});
0201     
0202     % get current orbit
0203     o=findorbit6Err(rerr,indBPM,inCOD);
0204     ox=o(1,:);
0205     oy=o(3,:);
0206     d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0207     dx=d(1,:);
0208     dy=d(3,:);
0209     
0210     % subtract reference orbit
0211     ox=ox-reforbit(1,:);
0212     oy=oy-reforbit(2,:);
0213     % subtract reference dispersion
0214     dx=dx-refdispersion(1,:);
0215     dy=dy-refdispersion(2,:);
0216     
0217     % weigths between orbit and dispersion
0218     ox=ox*(1-wdisp);
0219     oy=oy*(1-wdisp);
0220     dx=dx*(wdisp);
0221     dy=dy*(wdisp);
0222     
0223     % build RMs
0224     if correctflags(1) && correctflags(2) % dpp and mean0
0225         RMH=[ [ormH{1}*(1-wdisp);drmH{1}*(wdisp);ones(size(indHCor))] [dppH'*(1-wdisp);dppHd'*(wdisp);0] ];
0226         RMV=[ [ormV{3}*(1-wdisp);drmV{3}*(wdisp);ones(size(indVCor))] [dppV'*(1-wdisp);dppVd'*(wdisp);0] ];
0227     elseif correctflags(1) && ~correctflags(2)% dpp no mean 0
0228         RMH=[ [ormH{1}*(1-wdisp);drmH{1}*(wdisp)] [dppH'*(1-wdisp);dppHd'*(wdisp)] ];
0229         RMV=[ [ormV{3}*(1-wdisp);drmV{3}*(wdisp)] [dppV'*(1-wdisp);dppVd'*(wdisp)] ];
0230     elseif ~correctflags(1) && correctflags(2) % mean0 no dpp
0231         RMH=[ormH{1}*(1-wdisp);drmH{1}*(wdisp);ones(size(indHCor))];
0232         RMV=[ormV{3}*(1-wdisp);drmV{3}*(wdisp);ones(size(indVCor))];
0233     elseif ~correctflags(1) && ~correctflags(2) % no dpp no mean0
0234         RMH=[ormH{1}*(1-wdisp);drmH{1}*(wdisp)];
0235         RMV=[ormV{3}*(1-wdisp);drmV{3}*(wdisp)];
0236     end
0237     
0238     % compute correction
0239     if correctflags(2) % mean 0
0240         dch=qemsvd_mod(RMH,[ox';dx';0],neigSteerer(1));
0241         dcv=qemsvd_mod(RMV,[oy';dy';0],neigSteerer(2));
0242     else % no constraint on correctors mean
0243         dch=qemsvd_mod(RMH,[ox';dx'],neigSteerer(1));
0244         dcv=qemsvd_mod(RMV,[oy';dy'],neigSteerer(2));
0245     end
0246     
0247     
0248     % get total correctors values and apply scaling
0249     if correctflags(1)
0250         hs=corh0-dch(1:end-1)*scalefactor;
0251         vs=corv0-dcv(1:end-1)*scalefactor;
0252         % energy deviation
0253         inCOD(5)=dch(end);
0254     else
0255         hs=corh0-dch*scalefactor;
0256         vs=corv0-dcv*scalefactor;
0257     end
0258     
0259     % limit steerers strengths
0260     if ~isempty(steererlimit)
0261         hs(abs(hs)>steererlimit(1))=steererlimit(1);
0262         vs(abs(vs)>steererlimit(2))=steererlimit(2);
0263     end
0264     
0265     % apply correction in lattice
0266     rcor=atsetfieldvalues(rerr,indHCor,'PolynomB',{1,1},hs);
0267     rcor=atsetfieldvalues(rcor,indVCor,'PolynomA',{1,1},vs);
0268     
0269     % lattice start point for next iteration
0270     rerr=rcor;
0271 end
0272 
0273 % get current orbit
0274 o=findorbit6Err(rcor,indBPM,inCOD);
0275 oxc=o(1,:);
0276 oyc=o(3,:);
0277 d=finddispersion6Err(rcor,indBPM,indrfc,alpha,delta,inCOD);
0278 dxc=d(1,:);
0279 dyc=d(3,:);
0280     
0281 Lh=atgetfieldvalues(rcor,indHCor,'Length');
0282 Lv=atgetfieldvalues(rcor,indVCor,'Length');
0283 hsL=hs.*Lh;
0284 vsL=vs.*Lv;
0285 
0286 if printouttext
0287     % display results
0288     disp(['      before' '    ' '-->' '    ' 'after'])
0289     disp(['oX: ' num2str(std(ox0-reforbit(1,:))*1e6,'%3.3f') ' -> ' num2str(std(oxc-reforbit(1,:))*1e6,'%3.3f') 'um']);
0290     disp(['oY: ' num2str(std(oy0-reforbit(2,:))*1e6,'%3.3f') ' -> ' num2str(std(oyc-reforbit(2,:))*1e6,'%3.3f') 'um']);
0291     disp(['dX: ' num2str(std(dx0-refdispersion(1,:))*1e3,'%3.3f') ' -> ' num2str(std(dxc-refdispersion(1,:))*1e3,'%3.3f') 'mm'])
0292     disp(['dY: ' num2str(std(dy0-refdispersion(2,:))*1e3,'%3.3f') ' -> ' num2str(std(dyc-refdispersion(2,:))*1e3,'%3.3f') 'mm'])
0293     disp(['    ' 'min' '    ' 'mean' '    ' 'max'])
0294     disp(['hs:'  num2str([min(hsL) mean(hsL) max(hsL)]*1e3,' %2.2f ') ' mrad'])
0295     disp(['vs:'  num2str([min(vsL) mean(vsL) max(vsL)]*1e3,' %2.2f ') ' mrad'])
0296     disp(['dpp: ' num2str(inCOD(5))])
0297 end
0298 
0299 
0300 end