qempanelset.m

function [qemres,semres,qemb]=qempanelset(dirname,full_partial)
%QEMPANELSET
%
%[QEMB,QEMRES,SEMRES]=QEMPANELSET((DIRNAME,full_partial)
%
% full_partial = 1 partial, 0 full
% 
%QEMB:     Table of structures with fields q,cor,bh,bv,eta
%   QEMB(1): with no errors
%   QEMB(2): with errors
%QEMRES:   Empty structure to save results
%

global MACHFS

initname=fullfile(MACHFS,'MDT','respmat','qeminit.mat');
try
    qemres=load(initname);
catch %#ok<*CTCH>
    qemres=struct();
end
qemres.bhscale=1;
if full_partial
    qemres.hlist=1:6:96;
    qemres.vlist=1:6:96;
else
    qemres.hlist=1:1:96;
    qemres.vlist=1:1:96;
end
%qemres.qpcode=6;        % 32 quad correctors
qemres.qpcode=106;      % 31 quad correctors (n4/c4 removed)
%qemres.skcode=7;		% 32 skew correctors
%qemres.skcode=10;		% 64 skew correctors
qemres.skcode=110;		% 65 skew correctors (s4/c4 added)
qemres.datadir=dirname;
qemres.opticsdir=fullfile(dirname,'optics','');

semres=struct();
semres.hlist=qemres.hlist;
semres.vlist=qemres.vlist;
semres.iaxemz=NaN(14,1);
semres.wdisp=0;
semres.bvscale=1;
semres.skewkeep=true(1,65);
semres.skewkeep(1)=false;       % Do not touch kicker bump compensation
semres.skewkeep(end)=false;     % Do not touch kicker bump compensation

mach=sr.model.getat(qemres.opticsdir,'reduce',true);

try
    scorhl=load_correction(fullfile(qemres.datadir,'sextcor.dat'),...
        3,qemres.opticsdir);
    mach=setsrmagsextus(mach,fullfile(qemres.datadir,'srmag.dat'),scorhl);
catch
    warning('qem:noinit','No sextupole corrector values in the data directory');
    mach=setsrmagsextus(mach,fullfile(qemres.datadir,'srmag.dat'));
end
try
    stgs=load(fullfile(qemres.datadir,'settings.mat'));
    qemres.tunes=stgs.tunes-fix(stgs.tunes);
catch
    qemres.tunes=[0.44 0.39];
end
qemres.at=reshape(atfittune(mach(:),qemres.tunes,'QD6','QF7'),size(mach));	% retune the model
%  qemres.at=mach;															% no retuning
try
    df=stgs.deltafrf;
catch
    df=0;
end
circ=findspos(qemres.at,length(qemres.at)+1);
[~,periods,~,harmnumber]=atenergy(qemres.at);
frf=harmnumber*2.997924e8/circ;
qemres.ct=-circ*df/(frf+df);
syncorb=findsyncorbit(qemres.at,qemres.ct);
qemres.dpp=syncorb(5,1);
[~,tunes]=atlinopt(qemres.at,qemres.dpp);
qemres.fractunes=periods*tunes-fix(periods*tunes);

dipfold=sr.fold(findcells(qemres.at(:),'Class','Bend'));
qemres.dipidx=reshape(dipfold([1 4],:),1,64);

qemres.qpidx=findcells(qemres.at(:),'Class','Quadrupole');
if numel(qemres.qpidx) == 258
    qemres.qdlowidx=qemres.qpidx([153 154]);
    qemres.qpidx([153 154])=[];
else
    qemres.qdlowidx=[];
end
qemres.qpfit=true(size(qemres.qpidx));  % keep all quads
qemres.ql=atgetfieldvalues(qemres.at(qemres.qpidx),'Length');

thinsext=atgetcells(qemres.at(:),'Class',@isthinsext,'ThinMultipole');
thicksext=atgetcells(qemres.at(:),'Class','Sextupole');
qemres.sextidx=reshape(find(thinsext|thicksext),1,[]);
sxl=atgetfieldvalues(qemres.at(qemres.sextidx),'Length');
sxl(~(sxl~=0))=1;   % Handles NaN correctly
qemres.sxl=sxl;

qemres.steerhidx=qemres.sextidx(selcor(8));
if ~isfield(qemres,'khgain'), qemres.khgain=ones(1,length(qemres.steerhidx)); end
if ~isfield(qemres,'khrot'), qemres.khrot=zeros(1,length(qemres.steerhidx)); end
qemres.hsteerl=qemres.sxl(selcor(8));

qemres.steervidx=qemres.sextidx(selcor(9));
if ~isfield(qemres,'kvgain'), qemres.kvgain=ones(1,length(qemres.steervidx)); end
if ~isfield(qemres,'kvrot'), qemres.kvrot=zeros(1,length(qemres.steervidx)); end
qemres.vsteerl=qemres.sxl(selcor(9));

elist=sort(selcor(qemres.qpcode));
qemres.qcoridx=qemres.sextidx(elist);
qemres.qcorl=qemres.sxl(elist);

qemres.bpmidx=findcells(qemres.at(:),'Class','Monitor');
nbpm=length(qemres.bpmidx);
if ~isfield(qemres,'bhgain'), qemres.bhgain=ones(nbpm,1); end
if ~isfield(qemres,'bvgain'), qemres.bvgain=ones(nbpm,1); end
if ~isfield(qemres,'bhrot'), qemres.bhrot=ones(nbpm,1); end
if ~isfield(qemres,'bvrot'), qemres.bvrot=ones(nbpm,1); end
qemres.brot=0.5*(qemres.bhrot+qemres.bvrot);

elist=sort(selcor(qemres.skcode));
qemres.skewidx=qemres.sextidx(elist);
qemres.skewl=qemres.sxl(elist);

b1idx=findcells(qemres.at(:),'FamName','B1H');
b2idx=findcells(qemres.at(:),'FamName','B2S');
qemres.id25idx=b1idx(22);
qemres.d9idx=b2idx(6);
qemres.d11lensidx=b2idx(8);
qemres.qmode=struct;

if nargout >= 3
    qemb1.kn=atgetfieldvalues(qemres.at(qemres.qpidx),'PolynomB',{2});
    qemb1.dipdelta=zeros(length(qemres.dipidx),1);
    qemb1.dxs=zeros(length(qemres.sextidx),1);
    qemb1.cor=zeros(length(qemres.qcoridx),1);
    qemb1.skewcor=zeros(length(qemres.skewidx),1);
    qemb1.dzs=zeros(length(qemres.sextidx),1);
    qemb1.diptilt=zeros(length(qemres.dipidx),1);
    qemb1.ks=zeros(length(qemres.qpidx),1);
    qemb1.at=qemat(qemres,qemb1,false);
    qemb(1)=qemoptics(qemres,qemb1);
    
    try
        [kn,qemb(2).dipdelta,qemb(2).dxs]=qemerrload(qemres,fullfile(qemres.opticsdir,'quadinierrors.mat'));
        qemb(2).kn=qemb(1).kn.*(1+srmodulation(kn));
    catch
        warning('qem:noinit','No initial errors in the optics directory');
        qemb(2).kn=qemb(1).kn;
        qemb(2).dipdelta=qemb(1).dipdelta;
        qemb(2).dxs=qemb(1).dxs;
    end
    qemb(2).cor=load_correction(fullfile(qemres.datadir,'magcor.dat'),...
        qemres.qpcode,qemres.opticsdir);
    qemb(2).skewcor=zeros(length(qemres.skewidx),1);
    qemb(2).dzs=zeros(length(qemres.sextidx),1);
    qemb(2).diptilt=zeros(length(qemres.dipidx),1);
    qemb(2).ks=zeros(length(qemres.qpidx),1);
    mach=atfittune(qemat(qemres,qemb(2),false),qemres.tunes,'QD6','QF7');      % retune the model
    qemb(2).kn=atgetfieldvalues(mach(qemres.qpidx),'PolynomB',{2});
end

    function ok=isthinsext(elem,fieldname,fieldvalue)
        if isfield(elem,fieldname) && strcmp(elem.(fieldname),fieldvalue)
            pb=find(elem.PolynomB ~= 0,1);
            ok=~isempty(pb) && (pb == 3);
        else
            ok=false;
        end
    end
end