atcorrectorbit

 Closed orbit correction.

 function [...
    rcor,...           1) corrected lattice
    inCOD,...          2) initial COD (dpp is stored here)
    hs,vs...           3,4) total steerers strengths after correction
    ]=atcorrectorbit(...
     rerr,...          1) AT lattice to correct
     indBPM,...        2) Nbx1 bpm indexes
     indHCor,...       3) Nhx1 hor. cor indexes
     indVCor,...       4) Nvx1 ver. cor indexes
     inCOD,...         5) 6x1 initial COD guess
     neigSteerer,...   6) 2xNiter eigenvectors for correction H and V at
                          each iteration (default: [Nh/2 Nv/2])
     correctflags,...  7) correct [dpp mean0](default: [true true])
     scalefactor,...   8) scale factor to correction (default: 1.0)
     ModelRM,...       9) ModelRM.Orb(H/V)Cor = 4x1 cell of orbit response matrix
                          ModelRM.Orb(H/V)DPP = 6x1 array of orbit
                          response to dpp
                          if [] compute RM (default: [])
     reforbit,...      10) 2xNbpm reference orbit to correct to (default 0*2xNb)
     steererlimit      11) 2x1 limit of steerers abs(steerer)<steererlimit
                           (default: [], no limits)
     printouttext      12) 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 orbit refx refy
 use atsetfieldvalues, atgetcells


see also: qemsvd_mod findorbit6Err getresponsematrices

0001 function [rcor,inCOD,hs,vs]=atcorrectorbit(...
0002     rerr,...
0003     indBPM,...
0004     indHCor,...
0005     indVCor,...
0006     inCOD,...
0007     neigSteerer,...
0008     correctflags,...
0009     scalefactor,...
0010     ModelRM,...
0011     reforbit,...
0012     steererlimit,...
0013     printouttext)
0014 %
0015 % Closed orbit correction.
0016 %
0017 % function [...
0018 %    rcor,...           1) corrected lattice
0019 %    inCOD,...          2) initial COD (dpp is stored here)
0020 %    hs,vs...           3,4) total steerers strengths after correction
0021 %    ]=atcorrectorbit(...
0022 %     rerr,...          1) AT lattice to correct
0023 %     indBPM,...        2) Nbx1 bpm indexes
0024 %     indHCor,...       3) Nhx1 hor. cor indexes
0025 %     indVCor,...       4) Nvx1 ver. cor indexes
0026 %     inCOD,...         5) 6x1 initial COD guess
0027 %     neigSteerer,...   6) 2xNiter eigenvectors for correction H and V at
0028 %                          each iteration (default: [Nh/2 Nv/2])
0029 %     correctflags,...  7) correct [dpp mean0](default: [true true])
0030 %     scalefactor,...   8) scale factor to correction (default: 1.0)
0031 %     ModelRM,...       9) ModelRM.Orb(H/V)Cor = 4x1 cell of orbit response matrix
0032 %                          ModelRM.Orb(H/V)DPP = 6x1 array of orbit
0033 %                          response to dpp
0034 %                          if [] compute RM (default: [])
0035 %     reforbit,...      10) 2xNbpm reference orbit to correct to (default 0*2xNb)
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 orbit refx refy
0048 % use atsetfieldvalues, atgetcells
0049 %
0050 %
0051 %see also: qemsvd_mod findorbit6Err getresponsematrices
0052 
0053 
0054 
0055 % response matrix kicks
0056 kval=1e-5;
0057 delta=1e-3;
0058 
0059 alpha=mcf(rerr);
0060 indrfc=find(atgetcells(rerr,'Frequency'));
0061 f0=rerr{indrfc(1)}.Frequency;
0062 
0063 % default arguments
0064 if nargin<12
0065     printouttext=true;
0066 end
0067 if nargin<11
0068     steererlimit=[];
0069 end
0070 
0071 if nargin<4
0072     if printouttext
0073         disp('get BPM and Correctors indexes'); end;
0074     indBPM=finc(atgetcells(rerr,'Class','Monitor'));
0075     indHCor=finc(atgetcells(rerr,'iscorH','H'));
0076     indVCor=finc(atgetcells(rerr,'iscorV','V'));
0077 end
0078 
0079 if nargin<5
0080     inCOD=[0 0 0 0 0 0]';
0081 end
0082 
0083 if nargin<6
0084     neigSteerer=[length(indHCor) length(indVCor)]/2;
0085 end
0086 
0087 if nargin<7
0088     correctflags=[true true];
0089 end
0090 
0091 if nargin<8
0092     if printouttext
0093         disp(' --- scale set to 1.0'); end;
0094     scalefactor=1.0;
0095 end
0096 
0097 if nargin<9
0098     if printouttext, disp(' --- computing orbit Response matrix'); end;
0099     ModelRM=[];
0100 end
0101 
0102 if nargin<10
0103     if printouttext, disp(' --- reference orbit = 0'); end;
0104     reforbit=zeros(2,length(indBPM));
0105 end
0106 
0107 if scalefactor<0 || scalefactor>1
0108     if printouttext
0109         disp(' --- scale factor out of range. Set to 1.0'); end;
0110     scalefactor=1.0;
0111 end
0112 
0113 if correctflags(1) % dpp correction
0114     rmsel=[1 2 3];
0115 else
0116     rmsel=[1 2];
0117 end
0118 
0119 % load or compute response matrix
0120 if isempty(ModelRM)
0121     % get orbit RM
0122     if printouttext
0123         disp('get orbit RM'); end;
0124     
0125     ModelRM=getresponsematrices(...
0126         rerr,...          %1 AT lattice
0127         indBPM,...      %2 bpm indexes in at lattice
0128         indHCor,...     %3 h cor indexes
0129         indVCor,...     %4 v cor indexes
0130         [],...     %5 skew cor indexes
0131         [],...     %6 quad cor indexes
0132         [],...     %7 sext cor indexes
0133         inCOD,...       %8 initial coordinates
0134         rmsel...      %9 specifiy rm to be computed
0135         );
0136     
0137     
0138     if ~correctflags(1) % dpp correction
0139         ModelRM.OrbHDPP=[];
0140         ModelRM.OrbVDPP=[];
0141     end
0142     
0143 end
0144 
0145 ormH=ModelRM.OrbHCor;
0146 ormV=ModelRM.OrbVCor;
0147 % kval=ModelRM.kval;
0148 dppH=ModelRM.OrbHDPP;
0149 dppV=ModelRM.OrbVDPP;
0150 % delta=ModelRM.delta;
0151 
0152 % get initial orbit
0153 o=findorbit6Err(rerr,indBPM,inCOD);
0154 ox0=o(1,:);
0155 oy0=o(3,:);
0156 
0157 %rerr0=rerr;
0158 
0159 % iterate correction
0160 Niter=size(neigSteerer,1);
0161 for iter=1:Niter
0162     
0163     if printouttext
0164         disp(['Orbit correction iter ' num2str(iter,'%d, ') 'n-eig: ' num2str(neigSteerer(iter,:),'%d, ')]);
0165     end
0166     
0167     % initial corrector strengths
0168     corh0=atgetfieldvalues(rerr,indHCor,'PolynomB',{1,1});
0169     corv0=atgetfieldvalues(rerr,indVCor,'PolynomA',{1,1});
0170     
0171     % get current orbit
0172     o=findorbit6Err(rerr,indBPM,inCOD);
0173     ox=o(1,:);
0174     oy=o(3,:);
0175     
0176     % subtract reference orbit
0177     ox=ox-reforbit(1,:);
0178     oy=oy-reforbit(2,:);
0179     
0180     % build RMs
0181     if correctflags(1) && correctflags(2) % dpp and mean0
0182         RMH=[ [ormH{1};ones(size(indHCor))] [dppH';0] ];
0183         RMV=[ [ormV{3};ones(size(indVCor))] [dppV';0] ];
0184     elseif correctflags(1) && ~correctflags(2)% dpp no mean 0
0185         RMH=[ ormH{1} dppH' ];
0186         RMV=[ ormV{3} dppV' ];
0187     elseif ~correctflags(1) && correctflags(2) % mean0 no dpp
0188         RMH=[ormH{1};ones(size(indHCor))];
0189         RMV=[ormV{3};ones(size(indVCor))];
0190     elseif ~correctflags(1) && ~correctflags(2) % no dpp no mean0
0191         RMH=ormH{1};
0192         RMV=ormV{3};
0193     end
0194     
0195     % compute correction
0196     if correctflags(2) % mean 0
0197         dch=qemsvd_mod(RMH,[ox';0],neigSteerer(iter,1));
0198         dcv=qemsvd_mod(RMV,[oy';0],neigSteerer(iter,2));
0199     else % no constraint on correctors mean
0200         dch=qemsvd_mod(RMH,ox',neigSteerer(iter,1));
0201         dcv=qemsvd_mod(RMV,oy',neigSteerer(iter,2));
0202     end
0203     
0204     
0205     % get total correctors values and apply scaling
0206     if correctflags(1)
0207         hs=corh0-dch(1:end-1)*scalefactor;
0208         vs=corv0-dcv(1:end-1)*scalefactor;
0209         % energy deviation
0210         dd=-dch(end);
0211     else
0212         hs=corh0-dch*scalefactor;
0213         vs=corv0-dcv*scalefactor;
0214     end
0215     
0216     % limit steerers strengths
0217     if ~isempty(steererlimit)
0218         hs(abs(hs)>steererlimit(1))=steererlimit(1);
0219         vs(abs(vs)>steererlimit(2))=steererlimit(2);
0220     end
0221     
0222     rtest=atsetfieldvalues(rerr,indHCor,'PolynomB',{1,1},hs);
0223     rtest=atsetfieldvalues(rtest,indVCor,'PolynomA',{1,1},vs);
0224     
0225    
0226     if correctflags(1)
0227         
0228         rtest=atsetfieldvalues(rtest,indrfc,'Frequency',f0-alpha*(dd)*f0);
0229         
0230         if printouttext
0231             disp(['Delta RF : ' num2str(-alpha*(dd)*f0) ' Hz']);
0232         end
0233     end
0234     
0235     %[~,t,~]=atlinopt(rtest,0,1);
0236     t=tunechrom(rtest,0);
0237     if not(isnan(t(1)) || isnan(t(2)))
0238         % apply correction in lattice
0239         rcor = rtest;
0240     else 
0241         rcor = rerr;
0242         if printouttext
0243             disp('Orbit correction fails, stop at previous step.')
0244         end
0245         break
0246     end
0247     
0248     % lattice start point for next iteration
0249     rerr=rcor;
0250 end
0251 
0252 % get current orbit
0253 o=findorbit6Err(rcor,indBPM,inCOD);
0254 oxc=o(1,:);
0255 oyc=o(3,:);
0256 Lh=atgetfieldvalues(rcor,indHCor,'Length');
0257 Lv=atgetfieldvalues(rcor,indVCor,'Length');
0258 hsL=hs.*Lh;
0259 vsL=vs.*Lv;
0260 
0261 if printouttext
0262     % display results
0263     disp(['      before' '    ' '-->' '    ' 'after'])
0264     disp(['oX: ' num2str(std(ox0-reforbit(1,:))*1e6,'%3.3f') ' -> ' num2str(std(oxc-reforbit(1,:))*1e6,'%3.3f') 'um']);
0265     disp(['oY: ' num2str(std(oy0-reforbit(2,:))*1e6,'%3.3f') ' -> ' num2str(std(oyc-reforbit(2,:))*1e6,'%3.3f') 'um']);
0266     disp(['    ' 'min' '    ' 'mean' '    ' 'max'])
0267     disp(['hs:'  num2str([min(hsL) mean(hsL) max(hsL)]*1e3,' %2.2f ') ' mrad'])
0268     disp(['vs:'  num2str([min(vsL) mean(vsL) max(vsL)]*1e3,' %2.2f ') ' mrad'])
0269     disp(['dpp: ' num2str(inCOD(5))])
0270 end