Description of atcorrectdispersion

0001 function [rcor,inCOD,qs,ss]=atcorrectdispersion(...
0002     rerr,...
0003     indBPM,...
0004     indQCor,...
0005     indSCor,...
0006     inCOD,...
0007     neigSteerer,...
0008     correctflags,...
0009     scalefactor,...
0010     ModelRM,...
0011     refdispersion,...
0012     correctorslimit,...
0013     printouttext)
0014 % function [...
0015 %    rcor,...           1) corrected lattice
0016 %    inCOD,...          2) initial COD (dpp is stored here)
0017 %    qs,ss...           3) required normal and skew quad. strengths (total)
0018 %    ]=atcorrectdispersion(...
0019 %     rerr,...          1) AT lattice to correct
0020 %     indBPM,...        2) Nbx1 bpm indexes
0021 %     indHCor,...       3) Nhx1 hor. cor indexes
0022 %     indVCor,...       4) Nvx1 ver. cor indexes
0023 %     inCOD,...         5) 6x1 initial COD guess
0024 %     neigSteerer,...   6) 2xNiter eigenvectors for correction H and V at
0025 %                          each iteration (default: [Nh/2 Nv/2])
0026 %     correctflags,...  7) correct [dpp mean0](default: [true true])
0027 %     scalefactor,...   8) scale factor to correction (default: 1.0)
0028 %     ModelRM,...       9) ModelRM.Orb(H/V)Cor = 4x1 cell of orbit response mat.
0029 %                          ModelRM.Orb(H/V)DPP = 6x1 array of orbit
0030 %                          response to dpp
0031 %                          if [] compute RM (default: [])
0032 %     refdispersion,...10) 2xNbpm reference dispersion to correct to
0033 %                           (default rerr dispersion)
0034 %     correctorslimit  11) 2x1 limit of steerers abs(steerer)<steererlimit
0035 %                           (default: [], no limits)
0036 %     printouttext     12) if 1 or true, display rms orbit
0037 %     )
0038 %
0039 % features impelemented:
0040 % limit correctors strengths
0041 % ddp correction
0042 % sum of steerers = 0
0043 % 6D dispersion with BPM errors
0044 % initial coordinates
0045 % correction to reference dispersions refdispersion
0046 % use atsetfieldvalues, atgetcells
0047 %
0048 %
0049 %see also: qemsvd_mod finddispersion6Err getresponsematrices
0050 
0051 
0052 
0053 % response matrix kicks
0054 kval=1e-5;
0055 delta=1e-3;
0056 
0057 alpha=mcf(rerr);
0058 indrfc=find(atgetcells(rerr,'Frequency'));
0059 f0=rerr{indrfc(1)}.Frequency;
0060 
0061 % default arguments
0062 if nargin<12
0063     printouttext=true;
0064 end
0065 if nargin<11
0066     correctorslimit=[];
0067 end
0068 
0069 if nargin<4
0070     if printouttext
0071         disp('get BPM and Correctors indexes'); end;
0072     indBPM=finc(atgetcells(rerr,'Class','Monitor'));
0073     indQCor=finc(atgetcells(rerr,'Class','Quadrupole'));
0074     indSCor=finc(atgetcells(rerr,'iscorS','S'));
0075 end
0076 
0077 if nargin<5
0078     inCOD=[0 0 0 0 0 0]';
0079 end
0080 
0081 if nargin<6
0082     neigSteerer=[length(indQCor) length(indSCor)]/2;
0083 end
0084 
0085 if nargin<7
0086     correctflags=[true true];
0087 end
0088 
0089 if nargin<8
0090     if printouttext
0091         disp(' --- scale set to 1.0'); end;
0092     scalefactor=1.0;
0093 end
0094 
0095 if nargin<9
0096     if printouttext, disp(' --- computing orbit Response matrix'); end;
0097     ModelRM=[];
0098 end
0099 
0100 if nargin<10
0101      if printouttext, disp(' --- reference dispersion = 0 V, rerr disp H'); end;
0102     refdispersion=zeros(size(indBPM),2);
0103     [l,~,~]=atlinopt(rerr,0,indBPM);
0104     refdispersion(1,:)=arrayfun(@(a)a.Dispersion(1),l);
0105 end
0106 
0107 if scalefactor<0 || scalefactor>1
0108     if printouttext
0109         disp(' --- scale factor out of range. Set to 1.0'); end;
0110     scalefactor=1.0;
0111 end
0112 
0113 
0114 if correctflags(1) % dpp correction
0115     rmsel=[9 10 11];
0116 else
0117     rmsel=[1 2];
0118 end
0119 
0120 % load or compute response matrix
0121 if isempty(ModelRM)
0122     % get orbit RM
0123     if printouttext
0124         disp('get orbit RM'); end;
0125     
0126     ModelRM=getresponsematrices(...
0127         rerr,...          %1 AT lattice
0128         indBPM,...      %2 bpm indexes in at lattice
0129         [],...     %3 h cor indexes
0130         [],...     %4 v cor indexes
0131         indSCor,...     %5 skew  cor indexes
0132         indQCor,...     %6 quad cor indexes
0133         [],...
0134         inCOD,...       %7 initial coordinates
0135         rmsel...        %8 specifiy rm to be computed
0136         );
0137     
0138     if ~correctflags(1)
0139         
0140         ModelRM.DispHDPP=[];
0141         ModelRM.DispVDPP=[];
0142     end
0143     
0144 end
0145 
0146 drmH=ModelRM.DispQCor;
0147 drmV=ModelRM.DispSCor;
0148 % kval=ModelRM.kval;
0149 dppH=ModelRM.DispHDPP;
0150 dppV=ModelRM.DispVDPP;
0151 % delta=ModelRM.delta;
0152 alpha=mcf(rerr);
0153 indrfc=find(atgetcells(rerr,'Frequency'));
0154      
0155 % get initial dispersion
0156 d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0157 dx0=d(1,:);
0158 dy0=d(3,:);
0159 
0160 qs0=atgetfieldvalues(rerr,indQCor,'PolynomB',{1,2});
0161 
0162 % iterate correction
0163 Niter=size(neigSteerer,1);
0164 for iter=1:Niter
0165     
0166     if printouttext
0167         disp(['Dispersion correction iter ' num2str(iter,'%d, ') ...
0168             'n-eig: ' num2str(neigSteerer(iter,:),'%d, ')]);
0169     end
0170     
0171     % initial corrector strengths
0172     corq0=atgetfieldvalues(rerr,indQCor,'PolynomB',{1,2});
0173     cors0=atgetfieldvalues(rerr,indSCor,'PolynomA',{1,2});
0174     
0175     % get current orbit
0176     d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0177     dx=d(1,:);
0178     dy=d(3,:);
0179     
0180     % subtract reference orbit
0181     dx=dx-refdispersion(1,:);
0182     dy=dy-refdispersion(2,:);
0183     
0184     % build RMs
0185     if correctflags(1) && correctflags(2) % dpp and mean0
0186         RMH=[ [drmH{1};ones(size(indQCor))] [dppH';0] ];
0187         RMV=[ [drmV{3};ones(size(indSCor))] [dppV';0] ];
0188     elseif correctflags(1) && ~correctflags(2)% dpp no mean 0
0189         RMH=[ drmH{1} dppH' ];
0190         RMV=[ drmV{3} dppV' ];
0191     elseif ~correctflags(1) && correctflags(2) % mean0 no dpp
0192         RMH=[drmH{1};ones(size(indQCor))];
0193         RMV=[drmV{3};ones(size(indSCor))];
0194     elseif ~correctflags(1) && ~correctflags(2) % no dpp no mean0
0195         RMH=drmH{1};
0196         RMV=drmV{3};
0197     end
0198     
0199     % compute correction
0200     if correctflags(2) % mean 0
0201         dch=qemsvd_mod(RMH,[dx';0],neigSteerer(1));
0202         dcv=qemsvd_mod(RMV,[dy';0],neigSteerer(2));
0203     else % no constraint on correctors mean
0204         dch=qemsvd_mod(RMH,dx',neigSteerer(1));
0205         dcv=qemsvd_mod(RMV,dy',neigSteerer(2));
0206     end
0207     
0208     
0209     % get total correctors values and apply scaling
0210     if correctflags(1)
0211         qs=corq0-dch(1:end-1)*scalefactor;
0212         ss=cors0-dcv(1:end-1)*scalefactor;
0213         % energy deviation
0214         dd=-dch(end);
0215     else
0216         qs=corq0-dch*scalefactor;
0217         ss=cors0-dcv*scalefactor;
0218     end
0219     
0220     % limit correctors strengths
0221     if ~isempty(correctorslimit)
0222         qs(abs(qs)>correctorslimit(1))=correctorslimit(1);
0223         ss(abs(ss)>correctorslimit(2))=correctorslimit(2);
0224     end
0225     
0226     % apply correction in lattice
0227     rcor=atsetfieldvalues(rerr,indQCor,'PolynomB',{1,2},qs);
0228     rcor=atsetfieldvalues(rcor,indSCor,'PolynomA',{1,2},ss);
0229     
0230     if correctflags(1)
0231        
0232         rcor=atsetfieldvalues(rcor,indrfc,'Frequency',f0-alpha*(dd)*f0);
0233         
0234         if printouttext
0235             disp(['Delta RF : ' num2str(-alpha*(dd)*f0) ' Hz']);
0236         end
0237     end
0238     
0239     % lattice start point for next iteration
0240     rerr=rcor;
0241 end
0242 
0243 % get current orbit
0244 d=finddispersion6Err(rcor,indBPM,indrfc,alpha,delta,inCOD);
0245 dxc=d(1,:);
0246 dyc=d(3,:);
0247 
0248 if printouttext
0249     % display results
0250     disp(['before' '    ' '-->' '    ' 'after'])
0251     disp(['dX: ' num2str(std(dx0-refdispersion(1,:))*1e3,'%3.3f') ' -> ' num2str(std(dxc-refdispersion(1,:))*1e3,'%3.3f') 'mm'])
0252     disp(['dY: ' num2str(std(dy0-refdispersion(2,:))*1e3,'%3.3f') ' -> ' num2str(std(dyc-refdispersion(2,:))*1e3,'%3.3f') 'mm'])
0253     disp(['    ' 'min' '    ' 'mean' '    ' 'max'])
0254     disp(['hs:'  num2str([min(qs-qs0) mean(qs-qs0) max(qs-qs0)],' %2.4f ') ' 1/m2'])
0255     disp(['vs:'  num2str([min(ss) mean(ss) max(ss)],' %2.4f ') ' 1/m2'])
0256     disp(['dpp: ' num2str(inCOD(5))])
0257 
0258 end
atcorrectdispersion

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
