Description of atRDTdispersioncorrection

0001 function [rcor,inCOD,qs,ss]=atRDTdispersioncorrection(...
0002     rerr,...
0003     r0,...
0004     indBPM,...
0005     indQCor,...
0006     indSCor,...
0007     inCOD,...
0008     neigSteerer,...
0009     correctflags,...
0010     scalefactor,...
0011     wdisp,...
0012     ModelRM,...
0013     steererlimit,...
0014     printouttext)
0015 % function [...
0016 %    rcor,...            1) corrected lattice
0017 %    inCOD,...           2) initial COD (dpp is stored here)
0018 %    hs,vs...            3) required steerers strengths (total)
0019 %    ]=atdispersionfreesteering(...
0020 %     rerr,...           1) AT lattice to correct
0021 %     r0, ...            2) 2xNbpm reference rdt to correct to
0022 %     indBPM,...         3) Nbx1 bpm indexes       (default: monitor)
0023 %     indQCor,...        4) Nqx1 quad. cor indexes (default: quadrupole)
0024 %     indSCor,...        5) Nsx1 skew. cor indexes (default: sextupole)
0025 %     inCOD,...          6) 6x1 initial COD guess  (default: 6x1 zero)
0026 %     neigSteerer,...    7) 2xNiter eigenvectors for correction H and V at
0027 %                          each iteration (default: [Nh/2 Nv/2])
0028 %     correctflags,...   8) correct [ mean0](default: [ true])
0029 %     scalefactor,...    9) scale factor to correction (default: 1.0)
0030 %     [wdisph wtunes wdispv],...
0031 %                       10) dispersion and tune weight:
0032 %                           dispersionH*wdisph and orbith*(1-wdisph-wtune)
0033 %                           dispersionV*wdispv and orbith*(1-wdispv)
0034 %                           (default: 0.7 0.1 0.7)
0035 %     ModelRM,...       11) ModelRM.Disp(N/S)Quad = 6x1 cell of dispersion
0036 %                           response mat. if [] compute RM (default: [])
0037 %                           (default 0*2xNb, or from r0 if reftune is r0)
0038 %     steererlimit      12) 2x1 limit of steerers abs(steerer)<steererlimit
0039 %                           (default: [], no limits)
0040 %     printouttext      13) if 1 or true, display rms orbit
0041 %     )
0042 %
0043 % features impelemented:
0044 % - limit correctors strengths
0045 % - ddp correction
0046 % - sum of steerers = 0
0047 % - 6D orbit with BPM errors
0048 % - initial coordinate
0049 % - correction to reference rdt tune dispersion from r0 lattice
0050 % - retrival of normal and skew quadrupole components also from alignment
0051 %   errors and rotations
0052 % - use atsetfieldvalues, atgetcells
0053 %
0054 %
0055 % http://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.14.034002
0056 %
0057 %see also: qemsvd_mod findorbit6Err getresponsematrices
0058 
0059 
0060 
0061 % response matrix kicks
0062 %kval=1e-5;
0063 delta=1e-3;
0064 
0065 % default arguments
0066 if nargin<13
0067     printouttext=true;
0068 end
0069 if nargin<12
0070     steererlimit=[];
0071 end
0072 
0073 if nargin<5
0074     if printouttext
0075         disp('get BPM and Correctors indexes'); end;
0076     indBPM=finc(atgetcells(rerr,'Class','Monitor'));
0077     indQCor=finc(atgetcells(rerr,'Class','Quadrupole'));
0078     indSCor=finc(atgetcells(rerr,'Class','Sextupole'));
0079 end
0080 
0081 if nargin<6
0082     inCOD=[0 0 0 0 0 0]';
0083 end
0084 
0085 if nargin<7
0086     neigSteerer=[length(indQCor) length(indSCor)]/2;
0087 end
0088 
0089 if nargin<8
0090     correctflags=true;
0091 end
0092 
0093 if nargin<9
0094     if printouttext
0095         disp(' --- scale set to 1.0'); end;
0096     scalefactor=1.0;
0097 end
0098 
0099 if nargin<10
0100     if printouttext, disp(' --- wdisph=0.7 wtune=0.1 wdispv=0.7'); end;
0101     wdisp=[.7 .1 .7];
0102 end
0103 
0104 if nargin<11
0105     if printouttext, disp(' --- computing orbit Response matrix'); end;
0106     ModelRM=[];
0107 end
0108 
0109 
0110 if scalefactor<0 || scalefactor>1
0111     if printouttext
0112         disp(' --- scale factor out of range. Set to 1.0'); end;
0113     scalefactor=1.0;
0114 end
0115 
0116 
0117 % load or compute response matrix
0118 if isempty(ModelRM)
0119     % get orbit RM
0120     if printouttext
0121         disp('get RM'); end;
0122     
0123     
0124         ModelRM=getresponsematrices(...
0125             rerr,...          %1 AT lattice
0126             indBPM,...      %2 bpm indexes in at lattice
0127             [],...     %3 h cor indexes
0128             [],...     %4 v cor indexes
0129             indSCor,...     %5 skew cor indexes
0130             indQCor,...     %6 quad cor indexes
0131             [],...
0132             inCOD,...       %7 initial coordinates
0133             [10 11 12]...        %8 specifiy rm to be computed
0134             );
0135         
0136     
0137 end
0138 
0139 % load RM computed by getresponsematrices
0140 
0141 drmQ=ModelRM.DispQCor;
0142 drmS=ModelRM.DispSCor;
0143 
0144 
0145 tuneQ=[ModelRM.TuneQCor{1};ModelRM.TuneQCor{2}];
0146 
0147 % quad RDT RM
0148 [~,~,ind]=EquivalentGradientsFromAlignments6D(rerr,inCOD);
0149 indAllQuad=ind;
0150 indAllSkew=ind;
0151 
0152 %indAllQuad=[indQCor indSCor];
0153 %indAllSkew=[indQCor indSCor];
0154 
0155 [respqx,respqz]=qemrdtresp_mod(rerr,indBPM,indAllQuad);    % RDT response matrix assumes K=1
0156 QL=atgetfieldvalues(rerr,indAllQuad,'Length');          % quadrupole lengths
0157 QL(QL==0)=1;% thin lens magnets
0158 
0159 % convert response from KL to K as for dispersion response matrix
0160 % this is needed to merge the RM with the dispersion RM in the final
0161 % computation of correction.
0162 lengthsmat=repmat(QL',length(indBPM),1);
0163 respqx=respqx.*lengthsmat;
0164 respqz=respqz.*lengthsmat;
0165 
0166 [~,qkcor]=ismember(indQCor,indAllQuad);
0167 rdtQ=[...
0168     real(respqx(:,qkcor));...
0169     imag(respqx(:,qkcor));...
0170     real(respqz(:,qkcor));...
0171     imag(respqz(:,qkcor))];
0172 
0173 
0174 % skew RDT RM
0175 [respsx,respsz]=semrdtresp_mod(rerr,indBPM,indAllSkew);    % RDT response matrix assumes K=1
0176 SL=atgetfieldvalues(rerr,indAllSkew,'Length');          % quadrupole lengths
0177 SL(SL==0)=1;% thin lens magnets
0178 lengthsmat=repmat(SL',length(indBPM),1);
0179 respsx=respsx.*lengthsmat;
0180 respsz=respsz.*lengthsmat;
0181 
0182 [~,skcor]=ismember(indSCor,indAllSkew);
0183 rdtS=[...
0184     real(respsx(:,skcor));...
0185     imag(respsx(:,skcor));...
0186     real(respsz(:,skcor));...
0187     imag(respsz(:,skcor))];
0188 
0189 
0190 inCOD=[0 0 0 0 0 0]';
0191 [l,t,~]=atlinopt(r0,0,indBPM);
0192 refdispersion=zeros(2,length(indBPM));
0193 refdispersion(1,:)=arrayfun(@(a)a.Dispersion(1),l);
0194 refdispersion(2,:)=arrayfun(@(a)a.Dispersion(3),l);
0195 reftune=t;
0196 
0197 [KQnoer,KSnoer,~]=EquivalentGradientsFromAlignments6D(r0,inCOD);
0198 %KQnoer=atgetfieldvalues(r0,indAllQuad,'PolynomB',{1,2});
0199 %KSnoer=atgetfieldvalues(r0,indAllSkew,'PolynomA',{1,2});
0200 
0201 fx=respqx*KQnoer;
0202 fz=respqz*KQnoer;
0203 rdtvecq=[...
0204     real(fx);...
0205     imag(fx);...
0206     real(fz);...
0207     imag(fz)]';
0208 
0209 fx=respsx*KSnoer;
0210 fz=respsz*KSnoer;
0211 rdtvecs=[...
0212     real(fx);...
0213     imag(fx);...
0214     real(fz);...
0215     imag(fz)]';
0216 
0217 refrdt(1,:)=rdtvecq;
0218 refrdt(2,:)=rdtvecs;
0219 
0220 
0221 
0222 % get rdt vectors to correct
0223 [KQi,KSi,~]=EquivalentGradientsFromAlignments6D(rerr,inCOD);
0224 %KQ=atgetfieldvalues(rerr,indAllQuad,'PolynomB',{1,2});
0225 %KS=atgetfieldvalues(rerr,indAllSkew,'PolynomA',{1,2});
0226 
0227 fx=respqx*KQi;
0228 fz=respqz*KQi;
0229 rq0=[...
0230     real(fx);...
0231     imag(fx);...
0232     real(fz);...
0233     imag(fz)]';
0234 
0235 fx=respsx*KSi;
0236 fz=respsz*KSi;
0237 rs0=[...
0238     real(fx);...
0239     imag(fx);...
0240     real(fz);...
0241     imag(fz)]';
0242 
0243 
0244 alpha=mcf(rerr);
0245 indrfc=find(atgetcells(rerr,'Frequency'));
0246 
0247 % get initial dispersion
0248 
0249 d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0250 dx0=d(1,:);
0251 dy0=d(3,:);
0252 
0253 % get initial tune
0254 [~,t0,~]=atlinopt(rerr,0,1);
0255 
0256 %rerr0=rerr;
0257  qs0=atgetfieldvalues(rerr,indQCor,'PolynomB',{1,2});
0258  ss0=atgetfieldvalues(rerr,indSCor,'PolynomA',{1,2});
0259     
0260 % iterate correction
0261 Niter=size(neigSteerer,1);
0262 for iter=1:Niter
0263     
0264     if printouttext
0265         disp(['RDT Disp. Tune Steering iter ' num2str(iter,'%d, ') ...
0266             ' n-eig: ' num2str(neigSteerer(iter,:),'%d, ') ...
0267             ' alpha: ' num2str(wdisp,'%2.2f ')]);
0268     end
0269     
0270     % initial corrector strengths
0271     corq0=atgetfieldvalues(rerr,indQCor,'PolynomB',{1,2});
0272     cors0=atgetfieldvalues(rerr,indSCor,'PolynomA',{1,2});
0273     
0274     
0275     % get current rdt vectors to correct
0276     [KQe,KSe,~]=EquivalentGradientsFromAlignments6D(rerr,inCOD);
0277     %KQ=atgetfieldvalues(rerr,indAllQuad,'PolynomB',{1,2});
0278     %KS=atgetfieldvalues(rerr,indAllSkew,'PolynomA',{1,2});
0279     
0280     fx=respqx*KQe;
0281     fz=respqz*KQe;
0282     rq=[...
0283         real(fx);...
0284         imag(fx);...
0285         real(fz);...
0286         imag(fz)]';
0287     
0288     fx=respsx*KSe;
0289     fz=respsz*KSe;
0290     rs=[...
0291         real(fx);...
0292         imag(fx);...
0293         real(fz);...
0294         imag(fz)]';
0295     
0296     % get current dispersion
0297     d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0298     dx=d(1,:);
0299     dy=d(3,:);
0300     % get current tune
0301     [~,t,~]=atlinopt(rerr,0,1);
0302     
0303     
0304     % subtract reference orbit
0305     rq=rq-refrdt(1,:);
0306     rs=rs-refrdt(2,:);
0307     % subtract reference dispersion
0308     dx=dx-refdispersion(1,:);
0309     dy=dy-refdispersion(2,:);
0310     % subtract reference tune
0311     t=t-reftune;
0312     
0313     % weigths between RDT, tune and dispersion
0314     rq=rq*(1-wdisp(1)-wdisp(2));
0315     rs=rs*(1-wdisp(3));
0316     dx=dx*(wdisp(1));
0317     dy=dy*(wdisp(3));
0318     t=t*(wdisp(2));
0319     
0320     % build RMs
0321     if  correctflags(1) % mean0 no dpp
0322         RMQ=[rdtQ*(1-wdisp(1)-wdisp(2));drmQ{1}*(wdisp(1));tuneQ*(wdisp(2));ones(size(indQCor))];
0323         %RMQ=[rdtQ*(1-wdisp(1));drmQ{1}*(wdisp(1));ones(size(indQCor))];
0324         RMS=[rdtS*(1-wdisp(3));drmS{3}*(wdisp(3));ones(size(indSCor))];
0325     elseif ~correctflags(1) % no dpp no mean0
0326         RMQ=[rdtQ*(1-wdisp(1)-wdisp(2));drmQ{1}*(wdisp(1));tuneQ*(wdisp(2))];
0327         %RMQ=[rdtQ*(1-wdisp(1));drmQ{1}*(wdisp(1))];
0328         RMS=[rdtS*(1-wdisp(3));drmS{3}*(wdisp(3))];
0329     end
0330     
0331     % compute correction
0332     if correctflags(1) % mean = 0
0333         vecq=[rq dx t 0]';
0334         %vecq=[rq dx 0]';
0335         vecs=[rs dy 0]';
0336     else % no constraint on correctors mean
0337         vecq=[rq dx t]';
0338         %vecq=[rq dx]';
0339         vecs=[rs dy]';
0340     end
0341     
0342     dcq=qemsvd_mod(RMQ,vecq,neigSteerer(iter,1));
0343     dcs=qemsvd_mod(RMS,vecs,neigSteerer(iter,2));
0344     
0345     % get total correctors values and apply scaling
0346     
0347     qs=corq0-dcq*scalefactor;
0348     ss=cors0-dcs*scalefactor;
0349     
0350     % limit steerers strengths
0351     if ~isempty(steererlimit)
0352         qs(abs(qs)>steererlimit(1))=steererlimit(1);
0353         ss(abs(ss)>steererlimit(2))=steererlimit(2);
0354     end
0355     
0356     % apply correction in lattice
0357     rcor=rerr;
0358     rcor=atsetfieldvalues(rcor,indQCor,'PolynomB',{1,2},qs);
0359     rcor=atsetfieldvalues(rcor,indSCor,'PolynomA',{1,2},ss);
0360     
0361     % lattice start point for next iteration
0362     rerr=rcor;
0363 end
0364 
0365 
0366 % get current rdt vectors to correct
0367 [KQ,KS,~]=EquivalentGradientsFromAlignments6D(rcor,inCOD);
0368 %KQ=atgetfieldvalues(rcor,indQCor,'PolynomB',{1,2});
0369 %KS=atgetfieldvalues(rcor,indAllSkew,'PolynomA',{1,2});
0370 
0371 fx=respqx*KQ;
0372 fz=respqz*KQ;
0373 rqc=[...
0374     real(fx);...
0375     imag(fx);...
0376     real(fz);...
0377     imag(fz)]';
0378 
0379 fx=respsx*KS;
0380 fz=respsz*KS;
0381 rsc=[...
0382     real(fx);...
0383     imag(fx);...
0384     real(fz);...
0385     imag(fz)]';
0386 
0387 % get current dispersion
0388 d=finddispersion6Err(rcor,indBPM,indrfc,alpha,delta,inCOD);
0389 dxc=d(1,:);
0390 dyc=d(3,:);
0391 % get current tune
0392 [~,tc,~]=atlinopt(rcor,0,1);
0393 
0394 
0395 if printouttext
0396     % display results
0397     disp(['        before' '    ' '-->' '    ' 'after'])
0398     disp(['rq: ' num2str(std(rq0-refrdt(1,:))*1e3,'%3.3f') ' -> ' num2str(std(rqc-refrdt(1,:))*1e3,'%3.3f') '']);
0399     disp(['rs: ' num2str(std(rs0-refrdt(2,:))*1e3,'%3.3f') ' -> ' num2str(std(rsc-refrdt(2,:))*1e3,'%3.3f') '']);
0400     disp(['dX: ' num2str(std(dx0-refdispersion(1,:))*1e3,'%3.3f') ' -> ' num2str(std(dxc-refdispersion(1,:))*1e3,'%3.3f') 'mm'])
0401     disp(['dY: ' num2str(std(dy0-refdispersion(2,:))*1e3,'%3.3f') ' -> ' num2str(std(dyc-refdispersion(2,:))*1e3,'%3.3f') 'mm'])
0402     disp(['tX: ' num2str((t0(1)-reftune(1)),'%3.3f') ' -> ' num2str((tc(1)-reftune(1)),'%3.3f') ''])
0403     disp(['tY: ' num2str((t0(2)-reftune(2)),'%3.3f') ' -> ' num2str((tc(2)-reftune(2)),'%3.3f') ''])
0404     disp(['    ' 'min' '    ' 'mean' '    ' 'max'])
0405     disp(['qs:'  num2str([min(qs-qs0) mean(qs-qs0) max(qs-qs0)]*1e0,' %2.2f ') ' 1/m2'])
0406     disp(['ss:'  num2str([min(ss-ss0) mean(ss-ss0) max(ss-ss0)]*1e0,' %2.2f ') ' 1/m2'])
0407     disp(['dpp: ' num2str(inCOD(5))])
0408     
0409 %     figure;
0410 %     subplot(4,1,1);
0411 %     plot(rq0-refrdt(1,:),'r'); hold on;
0412 %     plot(rqc-refrdt(1,:),'b');
0413 %     ylabel('quad RDT')
0414 %     legend('before','after')
0415 %     subplot(4,1,2);
0416 %     plot(rs0-refrdt(2,:),'r'); hold on;
0417 %     plot(rsc-refrdt(2,:),'b');
0418 %     ylabel('skew RDT')
0419 %     subplot(4,1,3);
0420 %     plot(dx0-refdispersion(1,:),'r'); hold on;
0421 %     plot(dxc-refdispersion(1,:),'b');
0422 %     ylabel('Hor. disp.')
0423 %     subplot(4,1,4);
0424 %     plot(dy0-refdispersion(2,:),'r'); hold on;
0425 %     plot(dyc-refdispersion(2,:),'b');
0426 %     ylabel('Ver. disp.')
0427 %     saveas(gca,['RDTdispCor' strrep(num2str(wdisp,'%2.2f_'),'.','p') '.fig']);
0428 %     export_fig(['RDTdispCor' strrep(num2str(wdisp,'%2.2f_'),'.','p') '.jpg']);
0429 
0430     
0431 end
0432 
0433 
0434 end
atRDTdispersioncorrection

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
