Commit c0c2e10b authored by Lee Carver's avatar Lee Carver

Now taking whole esrf_matlab folder to save time

parent bb3327c4
function [kicks,eff,first_bpm,nb_bpm]=bump_coefs(resp,sb,sk,periods,flag)
%[kicks,eff,first_bpm,nb_bpm]=BUMP_COEFS(resp,sb,sk,periods,averaging)
% computes local bumps
%
% resp: response matrix
% sb: BPM abscissa
% sk: steerer abscissa
% periods: number of periods in the machine
% averaging: (optional) uses average over periods if true
if nargin < 5, flag=false; end
[nbpm,nc]=size(resp); %#ok<ASGLU>
nkpercell=floor(nc/periods); % frf possibly added to real correctors
ncor=periods*nkpercell;
if nc > ncor
kf=nc*ones(ncor,1);
else
kf=[];
end
copyr=ones(1,periods);
% maxpu=7;
maxpu=50;
k=(1:ncor)';
refk=[circshift(k,1) k circshift(k,-1) kf];
nk=size(refk,2);
kicks=zeros(nk,ncor);
eff=NaN(maxpu,ncor);
first_bpm=zeros(ncor,1);
nb_bpm=zeros(ncor,1);
% compute non-normalized bumps (least square)
for i=1:ncor
[kicks(:,i),eff2,in_bpm] = local_bump(resp,sb,sk,refk(i,:));
neff=min([maxpu;length(eff2)]);
eff(1:neff,i)=eff2(1:neff);
nb_bpm(i)=neff;
if neff > 0, first_bpm(i)=in_bpm(1);end
end
% compute average values over periods
meaneff=mean2(reshape(eff,nkpercell*maxpu,periods),2);
avereff=reshape(meaneff*copyr,maxpu,ncor);
averkick=reshape(mean2(reshape(kicks,nkpercell*nk,periods),2)*copyr,nk,ncor);
undef=find(~isfinite(eff)); % replace missing with average
eff(undef)=avereff(undef);
[fs,js]=sort(-reshape(meaneff,maxpu,nkpercell));%#ok<ASGLU> %compute normalization values
k=reshape(js(1,:)'*copyr,1,ncor)+maxpu*(0:ncor-1);
if flag
best=avereff(k(:)); % most efficient AVERAGE
kicks=averkick./(ones(nk,1)*best'); % normalize kicks
eff=avereff'./(best*ones(1,maxpu)); % normalize efficiencies
else
best=eff(k(:)); % most efficient INDEPENDENT
kicks=kicks./(ones(nk,1)*best'); % normalize kicks
eff=eff'./(best*ones(1,maxpu)); % normalize efficiencies
end
function [varargout] = celldeal(cellv)
%CELLDEAL Splits a cell array in separate variables
%
%The last output variable contains the rest of the input cell array
lg=length(cellv);
if nargout < lg
varargout=[cellv(1:nargout-1) {cellv(nargout:end)}]; %#ok<VARARG>
else
varargout=[cellv repmat({''},1,nargout-lg)]; %#ok<VARARG>
end
end
function [celnum,id,kdx]=getcellnum(npercell,ncells,celloffset,arg1,arg2)
%GETCELLNUM gives cell number and index in cell
%
%[CELL,NUM,KDX]=GETCELLNUM(NPERCELL,NCELLS,CELLOFFSET,IDX) returns the cell
% number and the index in cell of the IDX'th elements of a vector of
% length NCELLS*NPERCELL, starting at the beginning of cell CELLOFFSET+1
%
% NPERCELL: Number of elements per cell
% NCELLS: Total number of cells
% CELLOFFSET: Number of cells to skip
% IDX: vector of indexes in the range (1..NCELLS*NPERCELL)
% starting in cell CELLOFFSET+1
%
%[CELL,NUM,KDX]=GETCELLNUM(NPERCELL,NCELLS,CELLOFFSET,MASK)
%
% MASK: logical mask length NCELLS*NPERCELL
%
%[CELL,NUM,KDX]=GETCELLNUM(NPERCELL,NCELLS,CELLOFFSET,[CELL NUM])
%[CELL,NUM,KDX]=GETCELLNUM(NPERCELL,NCELLS,CELLOFFSET,CELL,NUM])
%
%See also: GETMASK, GETINDEX
if islogical(arg1) % Logical mask
arg1=find(arg1(:));
end
npc=zeros(1,ncells);
npc(:)=reshape(npercell,[],1);
coff=[0 cumsum(npc)];
if nargin >= 5 % cell, id
celnum=zeros(max([size(arg1);size(arg2)]));
id=celnum;
celnum(:)=arg1;
id(:)=arg2;
kdx=reshape(coff(mod(celnum-celloffset-1,32)+1),size(celnum))+id;
elseif size(arg1,2) == 2 % [cell id]
id=arg1(:,2);
celnum=arg1(:,1);
kdx=reshape(coff(mod(celnum-celloffset-1,32)+1),size(celnum))+id;
else % idx
idx=arg1;
cid=arrayfun(@(id) find(coff(2:end)>=id,1),idx)-1;
id=idx-reshape(coff(cid+1),size(idx));
celnum=mod(cid+celloffset,ncells)+1;
kdx=mod(idx+coff(celloffset+1)-1,sum(npc))+1;
end
end
function [pth,args,mach] = getmodelpath(varargin)
%GETMODELPATH % Get the full path of an optics directory
%
%[PTH,ARGS,MACH]=GETMODELPATH(VARARGIN) scans the input arguments looking for:
%
%- 'tl2'|'sy'|'ebs'
% machine name, path defaults to $(APPHOME)/mach/optics/settings/theory
%- ['tl2'|'sy'|'ebs',]'opticsname'
% machine and optics name
%- ['tl2'|'sy'|'ebs',]'/machfs/appdata/sr/optics/settings/opticsname':
% full path of optics directory
%
%PTH: Full path of the optics directory
%ARGS: remaining arguments
%MACH: 'sr', 'sy' or 'ebs'
global APPHOME
narg=1;
mach='ebs';
pth='';
while narg <= nargin && isempty(pth)
arg=varargin{narg};
if ~ischar(arg) || any(strcmpi(arg,{'h','v','x','z','h2v','v2h'}))
break;
else
if any(strcmpi(arg,{'tl2','sy','ebs'}))
mach=lower(arg);
else
pth=arg;
end
end
narg=narg+1;
end
if isempty(pth)
pth=fullfile(APPHOME,'optics',mach,'theory');
elseif isempty(fileparts(pth))
pth=fullfile(APPHOME,'optics',mach,pth);
else
dirs=regexp(pth,filesep,'split');
if any(strcmp('sy',dirs))
mach='sy';
elseif any(strcmp('ebs',dirs))
mach=ebs;
end
end
args=varargin(narg:end);
end
function varargout=load_corcalib(varargin)
%LOAD_CORCALIB loads corrector calibration
%
%[H,V,S,Q]=LOAD_CORCALIB(OPTICS,...);
%
%The OPTICS argument may be one of the following:
%
%- 'sy', 'ebs'
%
%H,V : hor./vert. steerer calibration in radians/A
%
%[COEF1,...]=LOAD_CORCALIB(MACH,PLANE1,...) returns the selected calibs
%
%PLANE : 'h' for horizontal steerers [1x1]
% 'v' for horizontal steerers [1x1]
% 's' for skew quad correctors [32x1]
% 'q' for normal quad correctors [1x1]
%
%
mach = varargin{1};
if strcmp(mach,'ebs') % default values
disp('no calibrations');
elseif strcmp(mach,'sy')
disp('Using default calibration');
h=3.79587e-3;
v=3.93469e-3;
q=[];
s=[];
end
varargout={h,v,s,q};
function [ oscill,aveorbit ] = load_dd(plane,varargin)
%LOAD_DD Load turn-by-turn position
%
%OSCILL=LOAD_DD(PLANE)
% PLANE: h, H, x, X, 1 => horizontal
% v, V, z, Z, 2 => vertical
% s,S,0 => sum
% OSCILL: Turn-by-turn positions (nturnsx224), 1st bpm is C4-1
%
%OSCILL=LOAD_DD(PLANE,KICKER)
% Sends forst a "On" command on KICKER and waits for the data
%
%[OSCILL,COD]=LOAD_DD(PLANE)
%OSCILL: Beam oscillation with average subtracted
%COD: average beam position
%
bpm=tango.cache('srdiag/beam-position/all', @tango.SrBpmDevice);
oscill=bpm.read_dd(plane,varargin{:})';
nturns=size(oscill,1);
disp([int2str(nturns) ' turns']);
if nargout >= 2
aveorbit=mean2(oscill);
oscill=oscill-aveorbit(ones(nturns,1),:);
end
function [ intens ] = load_sumdd()
%LOAD_sumDD Load turn-by-turn sum signal
%
%INTENS=LOAD_DD()
%
% INTENS= Turn-by-turn summ signal (nturnsx224), 1st bpm is C4-1
%
bpm=tango.cache('srdiag/beam-position/all', @tango.SrBpmDevice);
intens=bpm.read_dd('s')';
function [kicks,eff,in,orbit]=local_bump(resp,sb,sk,columns,ampl)
%[kicks,eff,orbit]=LOCAL_BUMP(resp,sb,sk,kindex) computes local bumps
% resp: response matrix
% sb: BPM positions
% sk: kicker positions
% kindex: vector of 3 steerer indexes
if columns(1) > columns(3)
out=sb<sk(columns(1)) & sb>sk(columns(3));
in=[find(sb>sk(columns(1))) ; find(sb<sk(columns(3)))];
else
out=sb<sk(columns(1)) | sb>sk(columns(3));
in=find( sb>sk(columns(1)) & sb<sk(columns(3)) );
end
out=find( out(:) & all(isfinite(resp'))' );
kicks=[1;-resp(out,columns(2:end))\resp(out,columns(1))];
eff=resp(in,columns)*kicks;
if nargin >= 5
best=-sort(-eff);
scale=ampl/best(1);
kicks=kicks*scale;
eff=eff*scale;
end
if nargout >= 4
orbit=resp(:,columns)*kicks;
% error=std(orbit(out));
end
function [DGL,GL0,GLf,quadDevice]=movetune(Q1,Q0,varargin)
%function [DGL,GL0,GLf,quadDevice]=movetune(Q1,Q0,varargin)
%
% returns currents to move SR tunes to Q1 starting from Q0
%
% INPUTS:
%
% Required
% Q1: [Qx,Qy] wished tune (measured)
% Q0: [Qx,Qy] initial tune. default is the measured tunes
%
% Optional (any order, string, value, string, value, ...)
% 'mode' : 'QF1QD2' : set tune using main families
% 'Matching' : move tune and keep optics locked
% 'lattice' : measured RM folder to execute qemextract (not yet there) or
% AT lattice
% (default: '', uses model lattice)
% 'setcur' : true/false(default), sets currents in the SR
%
% OUTPUTS:
% DGL: delta quadruple integrated strengths
% GL0: initial quadruple integrated strengths for Q0
% GLf: final quadrupole integrated strengths for Q1
% quadDevice: tango devices names
%
% example 1): default usage
% [DGL,GL0,GLf]=movetune([.20,.40],[.21 .34]);
%
% example 2): change mode and use measured lattice
% (not yet done)
%
%see also: atfittune
hostname=getenv('TANGO_HOST');%['Tango://ebs-simu:10000/'];
% PROBLEM HERE. ADDPATH, functions needed not in +ebs
setenv('TANGO_HOST', hostname);
% parse inputs
p = inputParser;
expectedmode={'QF1QD2','Matching'};
defaultmode=expectedmode{1};
defaultlattice='';% no errors lattice
defaultsetcur=false;% no setting of currents
% PROBLEM HERE. defaults from tango.Device
% dd=tango.Device([hostname 'srdiag/beam-tune/1']);
% defaulttunetocorrect=[dd.Tune_h.read dd.Tune_v.read];
addRequired(p,'Q1',@isnumeric);
addRequired(p,'Q0',@isnumeric);
addOptional(p,'mode',defaultmode,@(x) any(validatestring(x,expectedmode)));
addOptional(p,'lattice',defaultlattice,@(x)(isstr(x) || iscell(x)));
addOptional(p,'setcur',defaultsetcur);
parse(p,Q1,Q0,varargin{:});
changetunemode = p.Results.mode;
lattice = p.Results.lattice;
setcur = p.Results.setcur;
disp(repmat('-',1,20))
disp(['Tunes from: ' num2str(Q0)])
disp([' to: ' num2str(Q1)])
disp(repmat('-',1,20))
% use lattice model if matching requested
if strcmp(changetunemode,expectedmode([2]))
disp('Matching mode. Using model lattice without errors.')
disp(['lattice: ' lattice ' will be ignored.'])
lattice='';
end
% load lattice
if isempty(lattice)
sr=ebs.model();
r=sr.ring;
elseif iscell(lattice)
r=lattice;
sr = ebs.model(r);
else
error('not implemented yet. not possible to use measured lattice model')
sr=ebs.model('/machfs/MDT/2020/Dec...');
r=sr.ring;
end
% [resp,~]=ebs.tuneresponse(sr,changetunemode);
% DGL = resp*(Q1-Q0)';
latmodel = sr;
tunemode = changetunemode;
% get quadurpole and DQ indexes
iq = latmodel.get(0,'qp'); % quadrupoles
iq8 = latmodel.get(0,'qf8d'); % tilted qf8d (varying dipole component, 50 slices)
idq = latmodel.get(0,'dq'); % dipole quadrupole (fixed dipole component, but may be splitted for markers)
allq_unsort = [iq;iq8;idq];
% get order of quadrupoles components in lattice
[allq,iqa] = sort(allq_unsort);
% get lattice with initial tunes
latinitunes=ebs.model(latmodel);%,'reduce',true);
latinitunes.settune(Q0,tunemode);
% get lattice with final tunes
latnewtunes=ebs.model(latmodel);%,'reduce',true);
latnewtunes.settune(Q1,tunemode);
% compute difference in gradients
DGL_=latnewtunes.get(1,'qp') - latinitunes.get(1,'qp');
D8GL=latnewtunes.get(1,'qf8d') - latinitunes.get(1,'qf8d');
DqGL=latnewtunes.get(1,'dq') - latinitunes.get(1,'dq');
DGL=[DGL_; D8GL; DqGL];
q_group = NaN(size(allq));
ii=1;
ig=1;
while ii<length(allq)
q_group(ii)=ig;
% if not subsequent or after subsequent (mid. marker), change group
if allq(ii+1)-allq(ii)~=1 && allq(ii+1)-allq(ii)~=2
ig=ig+1;
end
ii=ii+1;
end
Nq= q_group(end-1)+1;
q_group(end) = Nq;
function dgls=sumk(dgl,qg)
% sum integrated gradient change for magnets of the same group
% (slices)
% dgl are the delta gradient, qg the quadrupole number, both
% in lattice order.
nq=max(qg);
dgls=zeros(nq,1);
for inq=1:nq
dgls(inq)=sum(dgl(qg==inq));% sum gradients
end
end
DGLsum = sumk(DGL(iqa),q_group)';
% get present strengths
quadDevice = 'srmag/m-q/all';
q = tango.Device(quadDevice);
GL0 = q.CorrectionStrengths.read;
GLf = GL0 + DGLsum;
if setcur
save('InitialQuadCorrectionStrengths','GL0');
% q.CorrectionStrenghts = GLf;
% set currents slowly.
p=GLf';
pact=GL0';
maxchange=0.01; p=GLf';
if max(abs(p-pact))>maxchange
[pmax]=max(abs(p-pact));
NP=floor(pmax/maxchange)+1;
disp(['large excursion, sharing in ' num2str(NP) ' steps'])
ptemp=pact;
% disp('initial values')
% disp(pact')
% disp('final values')
% disp(p')
%
for stepset=1:NP
ptemp=ptemp+(p-pact)./(NP);
% disp(['current values, step ' num2str(stepset) '/' num2str(NP)])
% disp(ptemp')
%
q.CorrectionStrengths = ptemp';
pause(0.5);
end
q.CorrectionStrengths = p';
else
% set in a singe shot
q.CorrectionStrengths = p';
end
end
end
function varargout=nselect(step,func,i,varargin)
%NSELECT decimate calls to a function
%
persistent is
if i==1
is=step-1;
end
if is==0
[varargout{1:nargout}]=func(i,varargin{:});
is=step;
else
varargout=cell(nargout,1);
end
is=is-1;
end
function r=responsem(bpm,st,nu,periods,alphal)
%R=RESPONSEM(OBS,EXC,NU) builds response matrix from EXC to OBS
% OBS,EXC: [beta Phi/2PiNu]
%
%R=RESPONSEM(OBS,EXC,{NU,ALPHA}) computes the response for a fixed RF freq.
% OBS,EXC: [beta Phi/2PiNu dispersion]
if iscell(nu)
alphal=nu{2};
nu=nu{1};
end
nbpm=size(bpm,1);
nst=size(st,1);
a=sqrt(bpm(:,1)*st(:,1)');
if nargin < 4
b=0.5 - abs(bpm(:,2*ones(nst,1)) - st(:,2*ones(nbpm,1))');
else
a=repmat(a,periods,periods);
k=0:(periods-1);
fbpm=reshape((bpm(:,2*ones(1,periods))+k(ones(1,nbpm),:)),nbpm*periods,1);
fst=reshape((st(:,2*ones(1,periods))+k(ones(1,nst),:)),nst*periods,1);
b=0.5 - abs(fbpm(:,ones(nst*periods,1)) - fst(:,ones(nbpm*periods,1))');
nu=nu*periods;
end
r=a.*cos(b*(2*pi*nu))./(2*sin(pi*nu));
if size(bpm,2)>=3 && size(st,2)>= 3
dispcor=bpm(:,3)*st(:,3)'./alphal;
if nargin >4
dispcor=repmat(dispcor,periods,periods);
end
r=r+dispcor;
end
function savetext(fname,format,values,header)
%SAVETEXT(FILENAME,FORMAT,VALUES,HEADER) saves matrix in user defined format
%
% FORMAT should describe one row of the matrix
% HEADER (optional) is a string printed before the values
fid=fopen(fname,'w');
if (nargin >= 4),if isstr(header)
for s=header', fprintf(fid,'%s\n',s); end
end,end
fprintf(fid,format,values');
fclose(fid);
function varargout=selectplane(plane,varargin)
%idx=selectplane(plane)
% plane selector: h|H|x|X|1, v|V|z|Z|2
%idx: 1 for horizontal
% 2 for vertical
% error for other
%
%[result1,result2,...]=selectplane(plane,selection1,selection2,...)
% plane selector: h|H|x|X|1, v|V|z|Z|2, v2h|V2H|3, h2v|H2V|4
%idx: selection(1) for horizontal
% selection(2) for vertical
% selection(3) for v2h
% selection(4) for h2v
%
%[...]=selectplane(...,'Choices',{choice1,choice2},...)
% Gives the possible choices for the plane selector. Each choice is a
% string or a cell array of strings
% default: {{'X','H'},{'Z','V'},'V2H','H2V'}
%
selection={{1,2}};
[choices,arguments]=getoption(varargin,'Choices',{{'x','h'},{'z','v'},'v2h','h2v'});
selection(1:length(arguments))=arguments;
nmax=min([length(choices) cellfun(@length,selection)]);
idx=0;
if ischar(plane)
for i=1:nmax
if any(strcmpi(plane,choices{i}))
idx=i;
break;
end
end
elseif isscalar(plane) && isnumeric(plane)
idx=plane;
end
if idx <= 0 || idx > nmax
error('cs:plane','plane should be %s',allchoices(choices(1:nmax)));
end
varargout=cellfun(@(arg) arg{idx},selection,'UniformOutput',false);
function res=allchoices(choices)
k=cellfun(@expand,num2cell(1:length(choices)),choices,'UniformOutput',false);
res=cat(2,k{:});
res=res(3:end);
function res=expand(i,choice)
if ischar(choice)
l={cat(2,'|',lower(choice),'|',upper(choice))};
else
l=cellfun(@(o) cat(2,'|',lower(o),'|',upper(o)),choice,'UniformOutput',false);
end
res=cat(2,', ',int2str(i),l{:});
end
end
end
function store_bpm(fname,bpm,bpm_name)
fid=fopen(fname,'w');
for i=1:size(bpm,1)
fprintf(fid,'%s\t%.6f\t%.4f\t%.6f\t%.4f\t%.4f\t%.6f\n',bpm_name{i},bpm(i,1:6));
end
fclose(fid);
function store_bump(fname,bump,st2name)
nb=size(bump,1);
if nb <= 0, return, end
st1name={st2name{nb} st2name{1:nb-1}};
st3name={st2name{2:nb} st2name{1}};
fid=fopen(fname,'w');
for i=1:nb
listf=5+(1:bump(i,5));
fprintf(fid,'%s\t%.3f\t',st1name{i},bump(i,1));
fprintf(fid,'%s\t%.3f\t',st2name{i},bump(i,2));
fprintf(fid,'%s\t%.3f\t',st3name{i},bump(i,3));
fprintf(fid,'%.0f\t%.0f',bump(i,4:5));
fprintf(fid,'\t%.6f',bump(i,listf));
fprintf(fid,'\n');
end
fclose(fid);
function store_steerer(fname,st,st_name)
fid=fopen(fname,'w');
for i=1:size(st,1)
fprintf(fid,'%s\t%.6f\t%.4f\t%.6f\t%.5e\n',st_name{i},st(i,1:4));
end
fclose(fid);
function store_variable(fname,name,value)
%STORE_VARIABLE Adds a variable to a .mat file
%
%STORE_VARIABLE(FILENAME,VARNAME,VALUE)
%
%FILENAME: Name of a .mat file on the MATLAB path
%VARNAME: Variable name
%VALUE: Value of the variable
try
v=load(fname);
catch %#ok<CTCH>
v=struct;
end
v.(name)=value; %#ok<STRNU>
save(fname,'-struct','v');
function y=sum2(x,dim)
%SUM2 Sum of elements. Ignores NaNs
dimx=size(x);
if nargin < 2