Description of atSkewRDTdispersioncorrection

0001 function [rcor,inCOD,ss]=atSkewRDTdispersioncorrection(...
0002     rerr,...
0003     r0,...
0004     indBPM,...
0005     indSCor,...
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 %     indSCor,...        4) Nsx1 skew. 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: [Nh/2 Nv/2])
0026 %     correctflags,...   7) correct [ mean0](default: [ true])
0027 %     scalefactor,...    8) scale factor to correction (default: 1.0)
0028 %     [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     indSCor=finc(atgetcells(rerr,'Class','Sextupole'));
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(indSCor)/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(' ---wdispv=0.7'); end;
0097     wdisp=.7;
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             indSCor,...     %5 skew cor indexes
0126             [],...     %6 quad cor indexes
0127             [],...
0128             inCOD,...       %7 initial coordinates
0129             11 );       %8 specifiy rm to be computed
0130             
0131         
0132     
0133 end
0134 
0135 % load RM computed by getresponsematrices
0136 
0137 drmS=ModelRM.DispSCor;
0138 
0139 % quad RDT RM
0140 [~,~,ind]=EquivalentGradientsFromAlignments6D(rerr,inCOD);
0141 indAllSkew=ind;
0142 
0143 
0144 % skew RDT RM
0145 [respsx,respsz]=semrdtresp_mod(rerr,indBPM,indAllSkew);    % RDT response matrix assumes K=1
0146 SL=atgetfieldvalues(rerr,indAllSkew,'Length');          % quadrupole lengths
0147 SL(SL==0)=1;% thin lens magnets
0148 lengthsmat=repmat(SL',length(indBPM),1);
0149 respsx=respsx.*lengthsmat;
0150 respsz=respsz.*lengthsmat;
0151 
0152 [~,skcor]=ismember(indSCor,indAllSkew);
0153 rdtS=[...
0154     real(respsx(:,skcor));...
0155     imag(respsx(:,skcor));...
0156     real(respsz(:,skcor));...
0157     imag(respsz(:,skcor))];
0158 
0159 
0160 inCOD=[0 0 0 0 0 0]';
0161 [l,~,~]=atlinopt(r0,0,indBPM);
0162 refdispersion=arrayfun(@(a)a.Dispersion(3),l);
0163 
0164 
0165 [~,KSnoer,~]=EquivalentGradientsFromAlignments6D(r0,inCOD);
0166 
0167 fx=respsx*KSnoer;
0168 fz=respsz*KSnoer;
0169 rdtvecs=[...
0170     real(fx);...
0171     imag(fx);...
0172     real(fz);...
0173     imag(fz)]';
0174 
0175 refrdt(1,:)=rdtvecs;
0176 
0177 
0178 
0179 % get rdt vectors to correct
0180 [~,KSi,~]=EquivalentGradientsFromAlignments6D(rerr,inCOD);
0181 
0182 fx=respsx*KSi;
0183 fz=respsz*KSi;
0184 rs0=[...
0185     real(fx);...
0186     imag(fx);...
0187     real(fz);...
0188     imag(fz)]';
0189 
0190 
0191 alpha=mcf(rerr);
0192 indrfc=find(atgetcells(rerr,'Frequency'));
0193 
0194 % get initial dispersion
0195 
0196 d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0197 dy0=d(3,:);
0198 
0199 %rerr0=rerr;
0200  ss0=atgetfieldvalues(rerr,indSCor,'PolynomA',{1,2});
0201     
0202 % iterate correction
0203 Niter=size(neigSteerer,1);
0204 for iter=1:Niter
0205     
0206     if printouttext
0207         disp(['RDT Disp. Tune Steering iter ' num2str(iter,'%d, ') ...
0208             ' n-eig: ' num2str(neigSteerer(iter,:),'%d, ') ...
0209             ' alpha: ' num2str(wdisp,'%2.2f ')]);
0210     end
0211     
0212     % initial corrector strengths
0213     cors0=atgetfieldvalues(rerr,indSCor,'PolynomA',{1,2});
0214     
0215     
0216     % get current rdt vectors to correct
0217     [~,KSe,~]=EquivalentGradientsFromAlignments6D(rerr,inCOD);
0218     %KQ=atgetfieldvalues(rerr,indAllQuad,'PolynomB',{1,2});
0219     %KS=atgetfieldvalues(rerr,indAllSkew,'PolynomA',{1,2});
0220     
0221     fx=respsx*KSe;
0222     fz=respsz*KSe;
0223     rs=[...
0224         real(fx);...
0225         imag(fx);...
0226         real(fz);...
0227         imag(fz)]';
0228     
0229     % get current dispersion
0230     d=finddispersion6Err(rerr,indBPM,indrfc,alpha,delta,inCOD);
0231     dy=d(3,:);
0232  
0233     % subtract reference orbit
0234     rs=rs-refrdt(1,:);
0235     % subtract reference dispersion
0236     dy=dy-refdispersion(1,:);
0237     % subtract reference tune
0238 
0239     % weigths between RDT, tune and dispersion
0240     rs=rs*(1-wdisp(1));
0241     dy=dy*(wdisp(1));
0242     
0243     % build RMs
0244     if  correctflags(1) % mean0 no dpp
0245         RMS=[rdtS*(1-wdisp(1));drmS{3}*(wdisp(1));ones(size(indSCor))];
0246     elseif ~correctflags(1) % no dpp no mean0
0247         RMS=[rdtS*(1-wdisp(1));drmS{3}*(wdisp(1))];
0248     end
0249     
0250     % compute correction
0251     if correctflags(1) % mean = 0
0252         vecs=[rs dy 0]';
0253     else % no constraint on correctors mean
0254         vecs=[rs dy]';
0255     end
0256     
0257     dcs=qemsvd_mod(RMS,vecs,neigSteerer(iter,1));
0258     
0259     % get total correctors values and apply scaling
0260     
0261     ss=cors0-dcs*scalefactor;
0262     
0263     % limit steerers strengths
0264     if ~isempty(steererlimit)
0265         ss(abs(ss)>steererlimit(1))=steererlimit(2);
0266     end
0267     
0268     % apply correction in lattice
0269     rcor=rerr;
0270     rcor=atsetfieldvalues(rcor,indSCor,'PolynomA',{1,2},ss);
0271     
0272     % lattice start point for next iteration
0273     rerr=rcor;
0274 end
0275 
0276 
0277 % get current rdt vectors to correct
0278 [~,KS,~]=EquivalentGradientsFromAlignments6D(rcor,inCOD);
0279 
0280 fx=respsx*KS;
0281 fz=respsz*KS;
0282 rsc=[...
0283     real(fx);...
0284     imag(fx);...
0285     real(fz);...
0286     imag(fz)]';
0287 
0288 % get current dispersion
0289 d=finddispersion6Err(rcor,indBPM,indrfc,alpha,delta,inCOD);
0290 dyc=d(3,:);
0291 
0292 
0293 if printouttext
0294     % display results
0295     disp(['        before' '    ' '-->' '    ' 'after'])
0296     disp(['rs: ' num2str(std(rs0-refrdt(1,:))*1e3,'%3.3f') ' -> ' num2str(std(rsc-refrdt(1,:))*1e3,'%3.3f') '']);
0297     disp(['dY: ' num2str(std(dy0-refdispersion(1,:))*1e3,'%3.3f') ' -> ' num2str(std(dyc-refdispersion(1,:))*1e3,'%3.3f') 'mm'])
0298     disp(['    ' 'min' '    ' 'mean' '    ' 'max'])
0299     disp(['ss:'  num2str([min(ss-ss0) mean(ss-ss0) max(ss-ss0)]*1e0,' %2.2f ') ' 1/m2'])
0300     disp(['dpp: ' num2str(inCOD(5))])
0301     
0302 %     figure;
0303 %     subplot(2,1,1);
0304 %     plot(rs0-refrdt(1,:),'r'); hold on;
0305 %     plot(rsc-refrdt(1,:),'b');
0306 %     legend('before','after')
0307 %     subplot(2,1,2);
0308 %     plot(dy0-refdispersion(1,:),'r'); hold on;
0309 %     plot(dyc-refdispersion(1,:),'b');
0310     
0311 end
0312 
0313 
0314 end
atSkewRDTdispersioncorrection

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
