Commit de79fa5d authored by Simone Liuzzo's avatar Simone Liuzzo

potentially different number of h and v steerers

parent 06243ba2
......@@ -22,9 +22,18 @@ khr=khrot(qemres.hlist);
khg=khgain(qemres.hlist);
kvr=kvrot(qemres.vlist);
kvg=kvgain(qemres.vlist);
bok=all(isfinite([resph resph2v respv2h respv]),2);
goodh = all(~isnan(resph)); % not disabled horizontal steerers
goodv = all(~isnan(respv)); % not disabled vertical steerers
bok=all(isfinite([...
resph(:,goodh)...
resph2v(:,goodh)...
respv2h(:,goodv)...
respv(:,goodv)]),2);
% Build model without errors
disp('computing model without errors')
[rh0,rh2v0]=qemresp(mach,qemres.ct,qemres.bpmidx,qemres.steerhidx(qemres.hlist),...
@(elem,ib,kick) qemsetsteerer(elem,kick,0),orbit0);
[rv2h0,rv0]=qemresp(mach,qemres.ct,qemres.bpmidx,qemres.steervidx(qemres.vlist),...
......@@ -39,12 +48,13 @@ bok=all(isfinite([resph resph2v respv2h respv]),2);
% qemres.bhscale*ones(size(qemres.bhgain)));
% [kvgain,bvgain]=qempanelgain(qemres.vlist,respv,rv,qemres.kvgain,...
% ones(size(qemres.bvgain)));
disp('computing steerers and BPM gains')
[khg,bhgain]=qempanelgain(resph,rh0,khg,qemres.bhscale*bhgain);
[kvg,bvgain]=qempanelgain(respv,rv0,kvg,bvgain);
for loop=2:5
[rhx,rh2vx,rv2hx,rvx]=qemcode(mach,qemres.ct,...
qemres.steerhidx(qemres.hlist),qemres.steerhidx(qemres.hlist),qemres.bpmidx(bok),...
qemres.steerhidx(qemres.hlist),qemres.steervidx(qemres.vlist),qemres.bpmidx(bok),...
khr(end,:),khg,kvr(end,:),kvg,[],[],[],orbit0(:,bok));
[resphx,resph2vx]=qembpmdecode(resph,resph2v,bhr(:,loop-1),bhgain,bvgain);
......
function [kgain,lgain]=qempanelgain(measresp,modelresp,kgini,bgini)
ok=isfinite(measresp);
bpmlist=sum(ok,2)>size(measresp,2)-1;
goodsteerer = all(~isnan(measresp));
bpmlist=sum(ok,2)>size(measresp(:,goodsteerer,1),2)-1; % sum only columns that are not NaN (disabled/frozen/fault correctors)
%nok=224-sum(bgini(~bpmlist));
sumgain=sum(bgini(bpmlist));
kgain=kgini;
......@@ -19,7 +20,7 @@ for count=2:5
[k2,l2]=meshgrid(kgain,lgain(bpmlist));
gain=(k2.*l2);
lincr=sum(gain.*xy,2)./sum(gain.*gain.*xx,2);
lincr=sum(gain(:,goodsteerer).*xy(:,goodsteerer),2)./sum(gain(:,goodsteerer).*gain(:,goodsteerer).*xx(:,goodsteerer),2);
% lincr=lincr*power(prod(lincr),-1/nok);
% lincr=lincr*(nok/sum(lincr)); % Keep <bgain>=0
lincr=lincr*(sumgain/sum(lgain(bpmlist).*lincr)); % Keep <bgain>=0
......
......@@ -4,68 +4,69 @@ function [rh,rh2v,rv2h,rv]=qemsteerdecode(mach,dct,qemres,resph,resph2v,respv2h,
%[RH,RH2V,RV2H,RV]=QEMSTEERDECODE(MACH,QEMRES,RHM,RH2VM,RV2HM,RVM,...
% KHROT,KHGAIN,KVROT,KVGAIN)
h_also_v_steerers = ismember(qemres.steerhidx,qemres.steervidx);
h_also_v = ismember(qemres.steerhidx,qemres.steervidx);
v_also_h = ismember(qemres.steervidx,qemres.steerhidx);
%remove disabled
h_also_v = h_also_v(qemres.hlist);
v_also_h = v_also_h(qemres.vlist);
honly = ~h_also_v;
vonly = ~v_also_h;
nbpm=size(resph,1);
chk=cos(khrot);
shk=sin(khrot);
cvk=cos(kvrot);
svk=sin(kvrot);
ck=cos(khrot(h_also_v_steerers)-kvrot);
hused = false(size(qemres.steerhidx));
hused(qemres.hlist)=true;
xh=hused(h_also_v_steerers);
qemres.hlist = find(hused(h_also_v_steerers)); %locally modify qemres.hlist
xv=false(size(qemres.steervidx));
xv(qemres.vlist)=true;
both=xh&xv;
honly=xh&(~xv);
vonly=xv&(~xh);
% both=false(1,96);
% honly=xh;
% vonly=xv;
vok=both(qemres.hlist);
% horizontal steerers gain and rotation
if sum(vok)>0
rh(:,vok)=(repmat(cvk(both).*khgain(both),nbpm,1).*resph(:,vok)-...
repmat(shk(both).*kvgain(both),nbpm,1).*respv2h(:,vok))./...
repmat(ck(both),nbpm,1);
rh2v(:,vok)=(repmat(cvk(both).*khgain(both),nbpm,1).*resph2v(:,vok)-...
repmat(shk(both).*kvgain(both),nbpm,1).*respv(:,vok))./...
repmat(ck(both),nbpm,1);
% horizontal and vertical steerers
if sum(h_also_v)>0
rh(:,h_also_v)=(repmat(cvk(v_also_h).*khgain(h_also_v),nbpm,1).*resph(:,h_also_v)-...
repmat(shk(h_also_v).*kvgain(v_also_h),nbpm,1).*respv2h(:,v_also_h))./...
repmat(chk(h_also_v),nbpm,1);
rh2v(:,h_also_v)=(repmat(cvk(v_also_h).*khgain(h_also_v),nbpm,1).*resph2v(:,h_also_v)-...
repmat(shk(h_also_v).*kvgain(v_also_h),nbpm,1).*respv(:,v_also_h))./...
repmat(chk(h_also_v),nbpm,1);
end
% only horizontal steerers
if sum(honly)>0
[rv2hmodel,rvmodel]=qemresp(mach,dct,qemres.bpmidx,qemres.steerhidx(honly),...
@(elem,ib,kick) qemsetsteerer(elem,0,kick),varargin{:});
rh(:,~vok)=(repmat(khgain(honly),nbpm,1).*resph(:,~vok)-...
rh(:,honly)=(repmat(khgain(honly),nbpm,1).*resph(:,honly)-...
repmat(shk(honly),nbpm,1).*rv2hmodel)./...
repmat(chk(honly),nbpm,1);
rh2v(:,~vok)=(repmat(khgain(honly),nbpm,1).*resph2v(:,~vok)-...
rh2v(:,honly)=(repmat(khgain(honly),nbpm,1).*resph2v(:,honly)-...
repmat(shk(honly),nbpm,1).*rvmodel)./...
repmat(chk(honly),nbpm,1);
end
hok=both(qemres.vlist);
% vertical steerers gain and rotation
if sum(hok)>0
rv2h(:,hok)=(repmat(chk(both).*kvgain(both),nbpm,1).*respv2h(:,hok)+...
repmat(svk(both).*khgain(both),nbpm,1).*resph(:,hok))./...
repmat(ck(both),nbpm,1);
rv(:,hok)=(repmat(chk(both).*kvgain(both),nbpm,1).*respv(:,hok)+...
repmat(svk(both).*khgain(both),nbpm,1).*resph2v(:,hok))./...
repmat(ck(both),nbpm,1);
% horizontal and vertical correctors
if sum(v_also_h)>0
rv2h(:,v_also_h)=(repmat(chk(h_also_v).*kvgain(v_also_h),nbpm,1).*respv2h(:,v_also_h)+...
repmat(svk(v_also_h).*khgain(h_also_v),nbpm,1).*resph(:,h_also_v))./...
repmat(cvk(v_also_h),nbpm,1);
rv(:,v_also_h)=(repmat(chk(h_also_v).*kvgain(v_also_h),nbpm,1).*respv(:,v_also_h)+...
repmat(svk(v_also_h).*khgain(h_also_v),nbpm,1).*resph2v(:,h_also_v))./...
repmat(cvk(v_also_h),nbpm,1);
end
% only vertical steerers
if sum(vonly)>0
[rhmodel,rh2vmodel]=qemresp(mach,dct,qemres.bpmidx,qemres.steervidx(vonly),...
@(elem,ib,kick) qemsetsteerer(elem,kick,0),varargin{:});
rv2h(:,~hok)=(repmat(kvgain(vonly),nbpm,1).*respv2h(:,~hok)+...
rv2h(:,vonly)=(repmat(kvgain(vonly),nbpm,1).*respv2h(:,vonly)+...
repmat(svk(vonly),nbpm,1).*rhmodel)./...
repmat(cvk(vonly),nbpm,1);
rv(:,~hok)=(repmat(kvgain(vonly),nbpm,1).*respv(:,~hok)+...
rv(:,vonly)=(repmat(kvgain(vonly),nbpm,1).*respv(:,vonly)+...
repmat(svk(vonly),nbpm,1).*rh2vmodel)./...
repmat(cvk(vonly),nbpm,1);
end
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment