Commit 25c8e56b authored by Simone Liuzzo's avatar Simone Liuzzo

initial sempanel commit

parent fbafc8eb
function [S,N]=dDxyDskew(r,indbpm,indquad,varargin)
% function [S,N]=dDxyDskew(r,indbpm,indquad,varargin)
%
% Computes the derivative of the dispersion with respect to dipole angles.
% units are [(m/[%])/rad]
%
% Inputs:
% r : AT lattice
% indbpm : BPM indexes
% indquad : quadrupole manget indexes (all the dipoles in r)
% magmodel : '1slice', 'thick' (default)
%
% Outputs:
% S: vertical dispersion derivative with respect to skew quadrupoles.
% N: horizontal dispersion derivative with respect to skew quadrupoles.
%
% Example:
% [S,N]=dDxyDskew(r,....
% find(atgetcells(r,'Class','Monitor'))',....
% find(atgetcells(r,'Class','Bend')),...
% 'magmodel','thick');
%
%see also: atlinopt
% parse inputs
expectedmodels={'thick','1slice'};
p=inputParser;
addRequired(p,'r');
addRequired(p,'indbpm');
addRequired(p,'indquad');
addParameter(p,'magmodel',expectedmodels{1},...
@(x) any(validatestring(x,expectedmodels)));
parse(p,r,indbpm,indquad,varargin{:});
magmodel = p.Results.magmodel;
switch magmodel
case '1slice'
% optics at center
machx2=cellfun(@(a)atdivelem(a,[0.5,0.5])',r,'un',0);
machx2=[machx2{:}, atmarker('end')]';
[lx2,~,~]=atlinopt(machx2,0,1:length(machx2)+1);
l=lx2(2:2:length(machx2)+1);
case 'thick'
% optics at entrance
[l,~,~]=atlinopt(r,0,1:length(r)+1);
end
nbpm=length(indbpm);
% dipoles angle and length
K_quad = cellfun(@(a)a.PolynomB(2),r(indquad),'un',1)';
len_quad = cellfun(@(a)a.Length,r(indquad),'un',1)';
% beta x at bpm and quadrupole
by_bpm = arrayfun(@(a)a.beta(2),l(indbpm));
by_quad = arrayfun(@(a)a.beta(2),l(indquad));
bx_bpm = arrayfun(@(a)a.beta(1),l(indbpm));
bx_quad = arrayfun(@(a)a.beta(1),l(indquad));
% alpha x at quadrupole
ay_quad = arrayfun(@(a)a.alpha(2),l(indquad));
ax_quad = arrayfun(@(a)a.alpha(1),l(indquad));
% dispersion at quadrupoles
dx_quad = arrayfun(@(a)a.Dispersion(1),l(indquad));
dxp_quad = arrayfun(@(a)a.Dispersion(2),l(indquad));
dy_quad = arrayfun(@(a)a.Dispersion(3),l(indquad));
dyp_quad = arrayfun(@(a)a.Dispersion(4),l(indquad));
% phase advance x at bpm and bends
phx_bpm = arrayfun(@(a)a.mu(1),l(indbpm));
phx_quad = arrayfun(@(a)a.mu(1),l(indquad));
phy_bpm = arrayfun(@(a)a.mu(2),l(indbpm));
phy_quad = arrayfun(@(a)a.mu(2),l(indquad));
% tunes including integer part
Q=l(end).mu/2/pi;
% define all quantities as matrices,
% quad related are identical rows,
% bpm related are identical columns
DXquad =repmat( dx_quad,nbpm,1);
DXPquad=repmat( dxp_quad,nbpm,1);
DYquad =repmat( dy_quad,nbpm,1);
DYPquad=repmat( dyp_quad,nbpm,1);
[PYquad,PYbpm] = meshgrid( phy_quad,phy_bpm);
[BYquad,BYbpm] = meshgrid(sqrt(by_quad),sqrt(by_bpm)); % square root here for speed
[PXquad,PXbpm] = meshgrid( phx_quad,phx_bpm);
[BXquad,BXbpm] = meshgrid(sqrt(bx_quad),sqrt(bx_bpm)); % square root here for speed
AYquad = repmat( ay_quad,nbpm,1);
AXquad = repmat( ax_quad,nbpm,1);
sqrtKq = sqrt(abs(K_quad));
sqrtKqL = sqrtKq.*len_quad;
Kq = repmat( (K_quad), nbpm,1);
aKq = repmat( abs(K_quad), nbpm,1);
signKq = Kq./aKq;
Kq32 = repmat( abs(K_quad).^(3/2), nbpm,1);
sqKq = repmat( sqrtKq, nbpm,1);
sKq = repmat( sin(sqrtKqL), nbpm,1);% sin/cos here for speed
cKq = repmat( cos(sqrtKqL), nbpm,1);
shKq = repmat( sinh(sqrtKqL), nbpm,1);% sinh/cosh here for speed
chKq = repmat( cosh(sqrtKqL), nbpm,1);
Lquad = repmat( len_quad, nbpm,1);
% define phase distance
dphY = (PYbpm-PYquad);
dphX = (PXbpm-PXquad);
tau_ymj = abs(dphY) - pi*Q(2);
tau_xmj = abs(dphX) - pi*Q(1);
function y=signtilde(x)
y=sign(x)-double(x==0);
end
% phase term
switch magmodel
case '1slice'
% hor disp/ d skew
JCmjx = BXquad.*DYquad.*cos(tau_xmj); % J_{C,mj}^{(Dy)}
% ver disp/ d skew
JCmjy = BYquad.*DXquad.*cos(tau_ymj); % J_{C,mj}^{(Dx)}
case 'thick'
% ver disp/ d skew
TSym = DXquad ./ (2.*aKq.*Lquad.*BYquad) .* ...
(sKq.*shKq + signKq.*cKq.*chKq - signKq) + ...
DXPquad ./ (2.*Kq32.*Lquad.*BYquad) .* ...
(sKq.*chKq - cKq.*shKq) ;
TCym = DXquad .*BYquad ./ (2.*sqKq.*Lquad) .* ...
(sKq.*chKq + cKq.*shKq) + ...
-AYquad .*DXPquad ./ (2.*Kq32.*Lquad.*BYquad) .* ...
(sKq.*chKq - cKq.*shKq) + ...
-AYquad .*DXquad ./ (2.*aKq.*Lquad.*BYquad) .* ...
(sKq.*shKq + signKq.* cKq.*chKq - signKq) + ...
+BYquad .*DXPquad ./ (2.*aKq.*Lquad) .* ...
(sKq.*shKq - signKq.* cKq.*chKq + signKq) ;
JCmjy = TCym.*cos(tau_ymj) + TSym.*signtilde(dphY).*sin(tau_ymj);
% hor disp/ d skew
% ver disp/ d skew
TSxm = DYquad ./ (2.*aKq.*Lquad.*BXquad) .* ...
(sKq.*shKq - signKq.*cKq.*chKq + signKq) + ...
DYPquad ./ (2.*Kq32.*Lquad.*BXquad) .* ...
(sKq.*chKq - cKq.*shKq) ;
TCxm = DYquad .*BXquad ./ (2.*sqKq.*Lquad) .* ...
(sKq.*chKq + cKq.*shKq) + ...
-AXquad .*DYPquad ./ (2.*Kq32.*Lquad.*BXquad) .* ...
(sKq.*chKq - cKq.*shKq) + ...
-AXquad .*DYquad ./ (2.*aKq.*Lquad.*BXquad) .* ...
(sKq.*shKq - signKq.* cKq.*chKq + signKq) + ...
+BXquad .*DYPquad ./ (2.*aKq.*Lquad) .* ...
(sKq.*shKq + signKq.* cKq.*chKq - signKq) ;
JCmjx = TCxm.*cos(tau_xmj) + TSxm.*signtilde(dphX).*sin(tau_xmj);
end
% dispersion derivative respect to bending angles
S = BYbpm./(2*sin(pi*Q(2))) .* ( JCmjy );
N = BXbpm./(2*sin(pi*Q(1))) .* ( JCmjx );
end
function [rh,rv]=dispresp(mach,rotidx,bpmidx)
dCT=1.0e-8;
[rh1,rv1,dpp1]=rotresp(mach,dCT,rotidx,bpmidx);
[rh0,rv0,dpp0]=rotresp(mach,0,rotidx,bpmidx);
norm=1/(dpp1-dpp0);
rh=(rh1-rh0)*norm;
rv=(rv1-rv0)*norm;
function [rh,rv,dpp]=rotresp(mach,dCT,rotidx,bpmidx)
rot=1.e-3;
norm=1/rot;
orbit=findsyncorbit(mach,dCT,bpmidx);
dpp=orbit(5);
refh=orbit(1,:)*norm;
refv=orbit(3,:)*norm;
rh=zeros(length(bpmidx),length(rotidx));
rv=zeros(length(bpmidx),length(rotidx));
for i=1:length(rotidx)
orbit=findsyncorbit(atsettilt(mach,rotidx(i),rot),dCT,bpmidx);
rh(:,i)=orbit(1,:)*norm - refh;
rv(:,i)=orbit(3,:)*norm - refv;
end
function newmach = semaddid(mach,cell,lmid)
%SEMADDID Add skew components at the location of ids
%NEWMACH=SEMADDID(MACH,CELL,LMID)
% creates a new AT structure with thin skew quads,
% labelled IDH1, IDH2, IDH3... at a distance lmid
% from the center of the straight section
%
% CELL: number of the ID straight section (defaut: all sections)
% LMID: Distance of IDs from the center of the straight section
% (default [0 1.8]
%
%See also: SEMADDIDMARKER
if nargin < 3, lmid=[0 1.8]; end
if nargin < 2, cell=0; end
if cell > 0
period=mod(cell-4,32)+1;
mach2=semaddidmarker(mach,skewq('ID%c%d'),lmid);
[o1,o2]=idbpm(mach,period);
[n1,n2]=idbpm(mach2,period);
if n2 > n1
newmach=[mach(1:o1) mach2(n1+1:n2-1) mach(o2:end)]';
else
newmach=[mach2(1:n2-1) mach(o2:o1) mach2(n1+1:end)]';
end
else
newmach=semaddidmarker(mach,skewq('ID%c%d'),lmid);
end
end
function sk=skewq(famnam)
sk.FamName=famnam;
sk.Length=0;
sk.PassMethod='ThinMPolePass';
sk.PolynomA=[0 0];
sk.PolynomB=[0 0];
sk.MaxOrder=1;
end
function [n1,n2]=idbpm(mach,period)
bpmlist=findcells(mach(:),'BetaCode','PU');
bpmlist=reshape(bpmlist([224 1:223]),7,32);
n1=bpmlist(1,period);
n2=bpmlist(2,period);
end
function newmach = semaddidmarker(mach,marker,loc)
%SEMADDIDMARKER Add skew components at the location of ids
%NEWMACH=SEMADDIDMARKER(MACH,MARKER,loc)
% creates a new AT structure with new elements inserted in straights,
%
% MARKER: new element to be inserted
% LOC: Distance of markers from middle of straight section
%
if isempty(marker)
marker=struct('FamName','MK%c%d','Length',0,'PassMethod','IdentityPass');
end
if loc(1) == 0
dist=diff(loc(:))';
else
dist=diff([0;loc(:)])';
end
namepat=marker.FamName;
lsub=loc(end);
high=findcells(mach(:),'FamName','SDH[I1]'); % "1" for 6m sections
low=findcells(mach(:),'FamName','SDL[O1]'); % "1" for 6m sections
if isempty(high) || isempty(low), return; end
sdh=atelem(mach{high(1)},'FamName','SDHK');
sdl=atelem(mach{low(2)},'FamName','SDLK');
seqh={};
seql={};
idn=0;
for lg=dist(end:-1:1)
idn=idn+1;
seqh=[seqh;{atelem(marker,'FamName',sprintf(namepat,'H',idn));atelem(sdh,'Length',lg)}]; %#ok<AGROW>
seql=[seql;{atelem(marker,'FamName',sprintf(namepat,'L',idn));atelem(sdl,'Length',lg)}]; %#ok<AGROW>
end
if loc(1) == 0
idn=idn+1;
seqh=[seqh;{atelem(marker,'FamName',sprintf(namepat,'H',idn))}];
seql=[seql;{atelem(marker,'FamName',sprintf(namepat,'L',idn))}];
end
for lg=dist
idn=idn+1;
seqh=[seqh;{atelem(sdh,'Length',lg);atelem(marker,'FamName',sprintf(namepat,'H',idn))}]; %#ok<AGROW>
seql=[seql;{atelem(sdl,'Length',lg);atelem(marker,'FamName',sprintf(namepat,'L',idn))}]; %#ok<AGROW>
end
t1=seqh(2*length(dist)+1:end,:);
t2=seql;
t3=seqh(1:2*length(dist),:);
ncells=size(mach,2);
lh=getcellstruct(mach,'Length',high);
ll=getcellstruct(mach,'Length',low);
newmach=setcellstruct(setcellstruct( mach,'Length',high,lh-lsub),'FamName',high,'SDHX');
newmach=setcellstruct(setcellstruct(newmach,'Length',low ,ll-lsub),'FamName',low ,'SDLX');
newmach=[t1(:,ones(1,ncells));newmach(high(1):low(1),:);...
t2(:,ones(1,ncells));newmach(low(2):high(2),:);...
t3(:,ones(1,ncells));newmach(high(2)+1:end,:)];
end
function mach=sembuildat(mach,qemres,kn,klncor,dipdelta,ks,klscor,diptilt)
%SEMBUILDAT Updates the AT structure with strengths, tilts
%
%MACH=SEMBUILDAT(MACH0,QEMRES,KN,NCOR,KS,SCOR,DIPTILT,DIPDELTA)
%
%MACH : initial machine AT structure
%QEMRES : fields bpmidx,qpidx,dipidx,qcoridx,skewidx,qcorl,skewl
%KN : Normal quadrupole strengths
%NCOR : Normal corrector strengths
%KS : Skew quadrupole strengths
%SCOR : Skew corrector strengths
%DIPTILT: Dipole tilt
%DIPTILT: Dipole field error
%
%MACH : Resulting AT structure
%
if nargin >= 8
mach=atsettilt(mach,qemres.dipidx,diptilt);
end
if nargin >= 7
mach=setcellstruct(mach,'PolynomA',qemres.skewidx,klscor./qemres.skewl,2);
end
if nargin >= 6
[k,tilts]=semtilt(kn,ks);
mach=atsettilt(mach,qemres.qpidx,tilts);
elseif nargin >= 3
k=kn;
end
if nargin >= 5
bend=getcellstruct(mach,'BendingAngle',qemres.dipidx);
mach=setcellstruct(mach,'BendingAngle',qemres.dipidx,(1+dipdelta).*bend);
end
if nargin >= 4
mach=setcellstruct(mach,'PolynomB',qemres.qcoridx,klncor./qemres.qcorl,2);
end
if nargin >= 3
mach=setcellstruct(mach,'K',qemres.qpidx,k);
mach=setcellstruct(mach,'PolynomB',qemres.qpidx,k,2);
end
function [fun,options]=semcorargs(value,mach,bpmidx,skewidx,options)
switch value
case 3 % rms of vertical size in all cells
bpmlist=reshape(1:224,7,32);
allidx=bpmidx(bpmlist(3,:));
[radmach,radindex,cavindex]=atradon(mach); %#ok<NASGU>
fun=@(x) fitrmsiax(x,radmach,skewidx,allidx,radindex);
options=optimset(options,'TolFun', 1.e-8);
case 4 % average mode emittance at BPM locations
[radmach,radindex,cavindex]=atradon(mach); %#ok<NASGU>
fun=@(x) fitmode(x,mach,radmach,skewidx,bpmidx,radindex);
options=optimset(options,'TolFun', 1.e-13);
case 5 % max of vertical size in all cells
bpmlist=reshape(1:224,7,32);
allidx=bpmidx(bpmlist(3,:));
[radmach,radindex,cavindex]=atradon(mach); %#ok<NASGU>
fun=@(x) fitmaxiax(x,radmach,skewidx,allidx,radindex);
options=optimset(options,'TolFun', 1.e-8);
case 6 % rms gamma
options=optimset(options,'TolFun', 1.e-10);
fun=@(x) fitstdgamma(x,mach,skewidx,bpmidx);
case 7 % rms vertical dispersion
fun=@(x) fitdisp(x,mach,skewidx,bpmidx);
options=optimset(options,'TolFun', 1.e-7);
otherwise
fun=[];
options=[];
end
function f=fitstdgamma(skewcor,mach,skewidx,bpmidx)
mch=setcellstruct(mach,'PolynomA',skewidx,skewcor,2);
lindata=linopt(mch,0,bpmidx);
gamma=cat(1,lindata.gamma)-1;
%f=max(abs(gamma));
f=sqrt(sum(gamma.*gamma)/length(bpmidx));
end
function f=fitdisp(skewcor,mach,skewidx,bpmidx)
mch=setcellstruct(mach,'PolynomA',skewidx,skewcor,2);
[lindata,tunes,xsi]=atlinopt(mch,0,bpmidx); %#ok<ASGLU,NASGU>
eta=cat(2,lindata.Dispersion)';
f=std(eta(:,3),1);
end
function f=fitrmsiax(skewcor,mach,skewidx,iaxidx,radidx)
mch=setcellstruct(mach,'PolynomA',skewidx,skewcor,2);
[envelope,espread,blength]=ohmienvelope(mch,radidx,iaxidx); %#ok<ASGLU,NASGU>
bm66=reshape(cat(2,envelope.R),6,6,[]);
vsize=squeeze(bm66(3,3,:));
f=sqrt(sum(vsize)/length(vsize));
end
function f=fitmaxiax(skewcor,mach,skewidx,iaxidx,radidx)
mch=setcellstruct(mach,'PolynomA',skewidx,skewcor,2);
[envelope,espread,blength]=ohmienvelope(mch,radidx,iaxidx); %#ok<ASGLU,NASGU>
beam66=reshape(cat(2,envelope.R),6,6,[]);
vsize=squeeze(beam66(3,3,:));
f=max(sqrt(vsize));
end
function f=fitmode(skewcor,mach,radmach,skewidx,bpmidx,radidx)
mch=setcellstruct(mach,'PolynomA',skewidx,skewcor,2);
lindata=linopt(mch,0,bpmidx);
mch=setcellstruct(radmach,'PolynomA',skewidx,skewcor,2);
[envelope,espread,blength]=ohmienvelope(mch,radidx,bpmidx); %#ok<ASGLU,NASGU>
modemit=NaN(length(lindata),2);
for i=1:length(lindata)
[beamA,beamB]=beam44(lindata(i)); % betatron 4x4 coupled matrices
bm66=envelope(i).R;
siginv=inv(bm66);
bm44=inv(siginv(1:4,1:4));
modemit(i,:)=([beamA(:) beamB(:)]\bm44(:))';
end
emittance=mean(modemit);
f=emittance(2);
end
end
function [ks,is]=semcorfit(semres1,semres2,opticsdir,neig)
% analyzes the difference of 2 response matrices in terms of corrector
% change
% ORBIT.*W = RESPONSE * (KS./V)
if nargin < 4, neig=64; end
resp=[semres2.respv(:)-semres1.respv(:);...
semres2.resph(:)-semres1.resph(:);...
semres2.frespz-semres1.frespz];
orbitrange=size(semres1.skewresponse,1)-224;
v=sr.fold(std(semres1.skewresponse(1:orbitrange,:),1,1)); % normalize response matrix
v=sr.unfold(repmat(1./sqrt(sum(v.*v,2)/32),1,32));
w=[ones(orbitrange,1);semres1.wdisp*ones(224,1)];
sk=semsolve4(semres1.skewresponse.*(w*v'),resp.*w,neig);
ks=mean(sk,2).*v;
if nargin >= 3
coef=load_corcalib(opticsdir,10);
is=ks./coef(:);
end
function stop=semcorout(x,iter,state,handles)
global interrupt
persistent tstart
stop=false;
switch state
case 'init'
interrupt=false;
tstart=tic;
set(handles.text3,'String',sprintf('Running (%d-%g)',...
iter.iteration,iter.fval));
case 'iter'
drawnow;
stop=interrupt;
if toc(tstart) > 5
tstart=tic;
set(handles.text3,'String',sprintf('Running (%d-%g)',...
iter.iteration,iter.fval));
end
end
function [resp,tunes]=semderiv(mach,dpp,fitlist,blist,hslist,vslist,dispfunc)
%SEMDERIV Compute derivatives of response matrix
%
%[resp,tunes]=semderiv(mach,fitidx,bpmidx,hsteeridx,vsteeridx)
%
%MACH: AT machine structure
%FITIDX: Index of varying skew quadrupoles
%BPMIDX: Index of BPMs
%HSTEERIDX: Index of horizontal steerers
%VSTEERIDX: Index of vertical steerers
%
%RESP: response matrix for skew quads: [nbpms*(nhsteer+nvsteer) length(fitidx)]
%TUNES: initial tunes
if nargin < 7, dispfunc=@(i,itot) 0; end
nq=length(fitlist);
nbpm=length(blist);
nhst=length(hslist);
nvst=length(vslist);
[v,j,kdx]=unique([fitlist blist hslist vslist]); %#ok<ASGLU>
qidx=kdx(1:nq);
bidx=kdx(nq+(1:nbpm));
hsidx=kdx(nq+nbpm+(1:nhst));
vsidx=kdx(nq+nbpm+nhst+(1:nvst));
alphal=findspos(mach,length(mach)+1)*mcf(mach);
[lindata,beta,mu,eta]=atavedata(mach,dpp,[v length(mach)+1]);
mutot=lindata(end).mu;
tunes=mutot/2/pi;
phase=mu./mutot(ones(length(lindata),1),:);
hkq=responsem([beta(qidx,1) phase(qidx,1) eta(qidx,1)],[beta(hsidx,1) phase(hsidx,1) eta(hsidx,1)],{tunes(1),alphal});
vqb=responsem([beta(bidx,2) phase(bidx,2)],[beta(qidx,2) phase(qidx,2)],tunes(2));
vkq=responsem([beta(qidx,2) phase(qidx,2)],[beta(vsidx,2) phase(vsidx,2)],tunes(2));
hqb=responsem([beta(bidx,1) phase(bidx,1) eta(bidx,1)],[beta(qidx,1) phase(qidx,1) eta(qidx,1)],{tunes(1),alphal});
resp=zeros(nbpm*(nhst+nvst),nq);
for ib=1:nq
rv=vqb(:,ib)*hkq(ib,:); % ATTENTION
rh=hqb(:,ib)*vkq(ib,:);
resp(:,ib)=[rv(:);rh(:)];
dispfunc(ib,nq);
end
function [resp,tunes]=semderivAnalytic(mach,dpp,fitlist,blist,hslist,vslist,dispfunc)
%SEMDERIV Compute derivatives of response matrix
%
%[resp,tunes]=semderiv(mach,fitidx,bpmidx,hsteeridx,vsteeridx)
%
%MACH: AT machine structure
%FITIDX: Index of varying skew quadrupoles
%BPMIDX: Index of BPMs
%HSTEERIDX: Index of horizontal steerers
%VSTEERIDX: Index of vertical steerers
%
%RESP: response matrix for skew quads: [nbpms*(nhsteer+nvsteer) length(fitidx)]
%TUNES: initial tunes
if nargin < 7, dispfunc=@(i,itot) 0; end
nq=length(fitlist);
nbpm=length(blist);
nhst=length(hslist);
nvst=length(vslist);
% [v,j,kdx]=unique([fitlist blist hslist vslist]); %#ok<ASGLU>
% qidx=kdx(1:nq);
% bidx=kdx(nq+(1:nbpm));
% hsidx=kdx(nq+nbpm+(1:nhst));
% vsidx=kdx(nq+nbpm+nhst+(1:nvst));
% alphal=findspos(mach,length(mach)+1)*mcf(mach);
%
% [lindata,beta,mu,eta]=atavedata(mach,dpp,[v length(mach)+1]);
% mutot=lindata(end).mu;
% tunes=mutot/2/pi;
% phase=mu./mutot(ones(length(lindata),1),:);
%
% hkq=responsem([beta(qidx,1) phase(qidx,1) eta(qidx,1)],[beta(hsidx,1) phase(hsidx,1) eta(hsidx,1)],{tunes(1),alphal});
% vqb=responsem([beta(bidx,2) phase(bidx,2)],[beta(qidx,2) phase(qidx,2)],tunes(2));
%
% vkq=responsem([beta(qidx,2) phase(qidx,2)],[beta(vsidx,2) phase(vsidx,2)],tunes(2));
% hqb=responsem([beta(bidx,1) phase(bidx,1) eta(bidx,1)],[beta(qidx,1) phase(qidx,1) eta(qidx,1)],{tunes(1),alphal});
disp('ANALYTIC QEMPANEL Response V2')
bndidx=find(atgetcells(mach,'BendingAngle'));
%bndidx=sort(bndidx([1:4:end,4:4:end]));% only long part
%% get anlytical RM response
[...
dX_dq, dY_dq, ...
dDx_dq, dDx_db, Dx, ...
dXY_ds, dYX_ds, ...
dDx_ds, dDy_ds,dDy_da...
]=CalcRespXXRespMat_thick_V2(...
mach',dpp,...
blist',... % bpm
hslist,... % correctors
fitlist,... % quadrupoles
bndidx,... % any index, unused.
fitlist);
resp=zeros(nbpm*(nhst+nvst),nq);
respa=resp;
LQ=atgetfieldvalues(mach,fitlist,'Length');
qind=find(atgetcells(mach,'Class','Quadrupole'))'; % response matrix computed always at all quadrupoles
quadforresponse=find(ismember(qind,fitlist)); % quadrupoles to use for fit amoung all
for iq=1:nq
% rv=vqb(:,ib)*hkq(ib,:); % ATTENTION
% rh=hqb(:,ib)*vkq(ib,:);
% resp(:,ib)=[rv(:);rh(:)];
ib = quadforresponse(iq); % use only selected quadrupoles
rva=dYX_ds(:,:,ib); % ANALYTIC
rha=dXY_ds(:,:,ib);
respa(:,iq)=[rva(:);rha(:)]./LQ(iq); %FACTOR LQ length of quadrupoles to make the 2 response identical.
dispfunc(iq,nq);
end
% % to test several times, run in command line >> global semres; semres=rmfield(semres,'quadresponse');
% figure;
% plot(resp(:));
% hold on;
% plot(respa(:));
% legend('analytic sempanel','analytic x L_{Quad}');
% storefigure('skewderivCompTest')
% use analytic for fit.
resp=respa;
\ No newline at end of file
function mess=semdisp(qemb,qemres,semres)
nbpm=length(qemres.bpmidx);
alllist=1:nbpm; % all BPMS
blist=circshift(reshape(alllist,[],32),[0 3]);
iaxcells=[5;10;11;14;18;21;25;26;29;31;3];
inairlist=blist(3,iaxcells); % IAX
d9=nbpm+1; % D9 pinhole
d11lens=nbpm+2; % D911 lens
id25=nbpm+3; % ID25 pinhole
iaxname=[arrayfun(@(i) sprintf('C%i',i),iaxcells,'UniformOutput',false);...
{'D9';'ID11-LENS';'ID25'}];
emittances=qemb.emittances;
sbpm=cat(1,qemb.lindata.SPos);
sskew=findspos(qemb.at(:),qemres.skewidx);
spinhole=[qemb.d9data.SPos;qemb.d11lensdata.SPos;qemb.id25data.SPos];
mess={...
sprintf('em. H [nm]: %7.3f %7.3f %7.3f',1.e9*mean(emittances(alllist,1:2:end)));...
sprintf('em. V [pm]: %7.3f %7.3f %7.3f',1.e12*mean(emittances(alllist,2:2:end)));...
sprintf('V. dispersion [m]: %g',qemb.pm.alpha*std(qemb.frespz,1));...
'in-air V [pm]:';...
sprintf(' %3d:%8.4f %3d:%8.4f\n',[iaxcells 1.e12*emittances(inairlist,6)]');...
sprintf('aver.:%8.4f\n',1.e12*mean(emittances(inairlist,6)));...
'pinhole V [pm]:';...
sprintf(' D9:%8.4f',1.e12*emittances(d9,6));...
sprintf(' D11:%8.4f',1.e12*emittances(d11lens,6));...
sprintf('ID25:%8.4f',1.e12*emittances(id25,6))};
figure(1); % emittances
plot(sbpm,1.e12*emittances(alllist,2:2:end));
hold on
plot(sbpm(inairlist),1.e12*emittances(inairlist,end),'ro','MarkerFaceColor','r');
plot(spinhole,1.e12*emittances([d9 d11lens id25],end),'ko','MarkerFaceColor','k');
hold off
ax=axis;
axis([0 844.39 0 max([ax(4);1.e-12])]);
plotskew(sskew);
xlabel('s [m]');
ylabel('\epsilon_z [pm]');
legend('Eigen emittance','projected emittance','measured emittance');
grid on
figure(2) % Comparison model / iax
bar(1.e12*[semres.iaxemz qemb.emittances([inairlist d9 d11lens id25],6)]);
set(gca,'XTick',1:length(iaxname),'XTickLabel',iaxname);
legend('measurement','model');
grid on
function plotskew(sskew)
ns=size(sskew,2);
ylim=get(gca,'Ylim');
hold('on');
plot([sskew;sskew],repmat([0 0;-.1 .1]*ylim',1,ns),'k');
hold('off');
function semdispresp(resp,hlist,vlist)
nh=length(hlist);
nv=length(vlist);
for col=1:size(resp,2)
orbvh(:,col)=std2(reshape(resp(:,col),224,nh+nv),1)';
end
subplot(2,1,1);
bar(orbvh(1:nh,:));
xlabel('H steerer');
ylabel('rms V response');
legend('Measured','Model');
subplot(2,1,2);
bar(orbvh(nh+(1:nv),:));
xlabel('V steerer');
ylabel('rms H response');
legend('Measured','Model');
function deltaset=semerrfit(nhst,nvst,resp,dresp,mode,okfit,okbpm)
%QEMPANELFIT Vary lattice parameters to fit the measured response matrix
%DELTAV=QEMPANELFIT(NHST,NVST,RESP,DRESP,MODE,OKVAR,OKBPM)
%
%NHST: Number of horizontal kicks
%NHST: Number of vertical kicks
%RESP: Deviation of measured reponse matrices from the model
%DRESP: Derivatives of response matrices vs. variable elements
%MODE: Structure allowing control of SVD inversion. (default:)
%OKVAR: Select a subset of variable parameters (default: all)
%OKBPM: Select a subset of valid BPMS (default:all)
%
%DELTAV: Variation of selected parameters
[no,nq]=size(dresp);
nbpm=no/(nhst+nvst+1);
orbitrange=nbpm*(nhst+nvst);
if nargin<7, okbpm=true(nbpm,1); end
if nargin<6, okfit=true(nq,1); end
if nargin<5 || isempty(mode), mode=struct; end
if ~isfield(mode,'nsets'), mode.nsets=4; end