Home > pubtools > LatticeTuningFunctions > correction > RDT > atRDTdispersionmeasuredcorrection.m

atRDTdispersionmeasuredcorrection

PURPOSE ^

function [...

SYNOPSIS ^

function [rcor,inCOD,qs,ss]=atRDTdispersionmeasuredcorrection(rerr,rfit,r0,indBPM,indQCor,indSCor,inCOD,neigSteerer,correctflags,scalefactor,wdisp,ModelRM,steererlimit,printouttext)

DESCRIPTION ^

 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, dispersion will be taken
                           from this lattice
     rfit,...           2) AT lattice with fitted errors
                           optics will be taken from this lattice
     r0, ...            3) 2xNbpm reference rdt to correct to
     indBPM,...         4) Nbx1 bpm indexes       (default: monitor)
     indQCor,...        5) Nqx1 quad. cor indexes (default: quadrupole)
     indSCor,...        6) Nsx1 skew. cor indexes (default: sextupole)
     inCOD,...          7) 6x1 initial COD guess  (default: 6x1 zero)
     neigSteerer,...    8) 2xNiter eigenvectors for correction H and V at
                          each iteration (default: [Nh/2 Nv/2])
     correctflags,...   9) correct [ mean0](default: [ true])
     scalefactor,...   10) scale factor to correction (default: 1.0)
     [wdisph wtunes wdispv],...
                       11) dispersion and tune weight:
                           dispersionH*wdisph and orbith*(1-wdisph-wtune)
                           dispersionV*wdispv and orbith*(1-wdispv)                          
                           (default: 0.7 0.1 0.7)
     ModelRM,...       12) ModelRM.Disp(N/S)Quad = 6x1 cell of dispersion 
                           response mat. if [] compute RM (default: [])
                           (default 0*2xNb, or from r0 if reftune is r0)
     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 rdt tune dispersion from r0 lattice
 - retrival of normal and skew quadrupole components also from alignment
   errors and rotations
 - use atsetfieldvalues, atgetcells


 http://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.14.034002

see also: qemsvd_mod findorbit6Err getresponsematrices

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Thu 24-Aug-2017 18:47:33 by m2html © 2005