Description of atRDTdispersionmeasuredcorrection

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
