Description of atQuadRDTdispersioncorrection

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
