Commit ae944e19 authored by Simone Liuzzo's avatar Simone Liuzzo

missing files for bumps in SY

parent dc45ab68
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 [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 [y,xok]=mean2(x,varargin)
%MEAN2 Average or mean value. Ignores NaN and outliers.
% For vectors, MEAN2(X) is the mean value of the elements in X. For
% matrices, MEAN2(X) is a row vector containing the mean value of
% each column. For N-D arrays, MEAN2(X) is the mean value of the
% elements along the first non-singleton dimension of X.
%
% MEAN2(X,DIM) takes the mean along the dimension DIM of X.
%
% Example: If X = [0 1 2
% 3 4 5]
%
% then mean2(X,1) is [1.5 2.5 3.5] and mean2(X,2) is [1
% 4]
%
% MEAN2(X,DIM,LIMIT) Eliminates points with deviation larger than LIMIT*STD(X)
%
% [Y,OK]=MEAN2(X,...) returns a logical array OK indicating the accepted values
%
% See also MEDIAN, STD2, MIN2, MAX2, VAR, COV, MODE.
[y,xs,xok]=statok2(x,1,varargin{:});
function [xm,xs,xok]=statok2(x,varargin)
dimx=size(x);
if nargin < 3, dim=min(find(dimx ~= 1)); else dim=varargin{2}; end
if nargin < 2, weight=1; else weight=varargin{1}; end
if isempty(dim), dim = 1; end
if dim == 1
if nargin < 4, thresh=NaN; else thresh=varargin{3}; end
nc=prod(dimx(2:end));
xp=reshape(x,dimx(1),nc);
xokp=isfinite(xp);
[xpm,xps]=deal(NaN*ones(1,nc));
for j=1:nc
xq=xp(:,j);
ok=xokp(:,j);
nok=sum(ok);
go=nok>0;
if nok
while go
xqm=sum(xq(ok),1)/nok;
xqs=sqrt(var(xq(ok),weight,1));
xq0=xq;
xq0(ok)=xqm;
out=abs(xq-xq0)>(thresh*xqs);
ok=ok & (~out);
nok=sum(ok);
go=(sum(out)>0) && (nok>0);
end
xpm(j)=xqm;
xps(j)=xqs;
xokp(:,j)=ok;
else
[xpm(j),xps(j)]=deal(NaN);
xokp(:,j)=false;
end
end
xm=reshape(xpm,[1 dimx(2:end)]);
xs=reshape(xps,[1 dimx(2:end)]);
xok=reshape(xokp,dimx);
else
perm=1:ndims(x); perm(dim)=[]; perm=[dim perm];
[xm1,xs1,xok1]= statok2(permute(x,perm),weight,1,varargin{3:end});
xm=ipermute(xm1,perm);
xs=ipermute(xs1,perm);
xok=ipermute(xok1,perm);
end
function [y,ok] = std2(varargin)
%STD2 Standard deviation, ignoring NaN values.
% For vectors, Y = STD2(X) returns the standard deviation. For matrices,
% Y is a row vector containing the standard deviation of each column. For
% N-D arrays, STD operates along the first non-singleton dimension of X.
%
% STD2 normalizes Y by N, where N is the sample size. It produces the square
% root of the second moment of the sample about its mean.
% STD2(X,1) is the same as STD(X).
%
% Y = STD2(X,0) normalizes by N-1. This is the sqrt of an unbiased estimator
% of the variance of the population from which X is drawn, as long as X
% consists of independent, identically distributed samples.
%
% Y = STD2(X,FLAG,DIM) takes the standard deviation along the dimension
% DIM of X. Pass in FLAG==1 to use the default normalization by N, or
% 0 to use N-1.
%
% Example: If X = [4 -2 1
% 9 5 7]
% then std2(X,0,1) is [3.5355 4.9497 4.2426] and std2(X,0,2) is [3.0
% 2.0]
% See also MEAN2.
if nargin > 3
x=varargin{1};
w=varargin{2};
dim=varargin{3};
ok=isfinite(x);
tile = ones(1,max(ndims(x),dim)); tile(dim) = size(x,dim);
threshold=varargin{4};
while true
y=sqrt(var2(x,w,dim));
overs=abs(x-repmat(mean2(x,dim),tile)) > repmat(threshold*y,tile);
if ~any(overs(:)), break; end
ok(overs)=false;
x(overs)=NaN;
end
else
y=sqrt(var2(varargin{:}));
end
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);
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