Commit bc6672e9 authored by ncarmignani's avatar ncarmignani

Merge branch 'master' of gitlab.esrf.fr:BeamDynamics/matlaboperation

parents 0ef507fe 1ecfe835
This diff is collapsed.
......@@ -9,8 +9,8 @@ tic;
disp(['---- ID ' num2str(idnum,'%.2d') ' -----']);
setpropbump('h','id',rc,idnum,filedata);
setpropbump('v','id',rc,idnum,filedata);
%setpropbump('h','id',rc,idnum,filedata);
%setpropbump('v','id',rc,idnum,filedata);
disp(['---- BM ' num2str(idnum,'%.2d') ' -----']);
setpropbump('h','bm',rc,idnum,filedata);
......@@ -21,8 +21,8 @@ tic;
for bpmnum=1:10
disp(['---- BPM ' num2str(idnum,'%.2d') 'c' num2str(bpmnum,'%.2d') ' -----']);
setpropbump_bpmpos('h',bpmnum,rc,idnum,filedata);
setpropbump_bpmpos('v',bpmnum,rc,idnum,filedata);
% setpropbump_bpmpos('h',bpmnum,rc,idnum,filedata);
% setpropbump_bpmpos('v',bpmnum,rc,idnum,filedata);
end
......@@ -66,18 +66,27 @@ end
switch beamline
case 'id'
D = 5.3028;
% normalized quantities to maximum bump amplitude
shn = sh'./max(h);
hn = h./max(h);
case 'bm'
% normalized quantities to maximum bump amplitude
%shn = sh'./abs(h(bpmupstream-selb(1)+1))*2;
%hn = h./abs(h(bpmupstream-selb(1)+1))*2;
shn = sh'./max(h);
hn = h./max(h);
switch plane
case 'h'
D = 0.9120 / 1.224;
D = 0.9120 * 1.224;
case 'v'
D = 0.9120 / 0.8022;
D = 0.9120 * 0.8022;
end
end
% normalized quantities to maximum bump amplitude
shn = sh'./max(h);
hn = h./max(h);
fullnames= cellfun(@(a)(['srmag/' plane 'st-' a '']),hnames,'un',0);
......
......@@ -18,6 +18,10 @@ classdef DynamicAperture < RingControl %& GridScan
scraper_positions % positions to step scraper in [mm]
khdev % device name of hor. kicker to use for amplitue scan
kedev % device name of extraction. kicker to use to trigger injection bump
khdevCurr0 % initial Khdev Current, used only if kickmode = fullbump
kickmode % "fullbump" = full bump + reduce khdev, "singlekicker" = zero bump + increase khdev
tolhk % A
diithres % loss rate %
......@@ -62,7 +66,7 @@ classdef DynamicAperture < RingControl %& GridScan
obj.dfreq = 0; % on energy DA only by default
obj.scraper_positions = -[4:1:11]; % in [mm]
obj.scraper_positions = -[4.5:0.5:7]; % in [mm]
% measurement results
obj.x = NaN(length(obj.dfreq));
......@@ -70,7 +74,20 @@ classdef DynamicAperture < RingControl %& GridScan
obj.o = NaN(length(obj.dfreq));
obj.q = NaN(length(obj.dfreq),2);
obj.khdev = 'sr/ps-k1/1';
obj.khdev = 'sr/ps-k4/4';
obj.kedev = 'sy/ps-ke/1';
obj.kickmode= 'fullbump'; % or 'singlekicker'
switch obj.kickmode
case 'fullbump'
obj.khdevCurr0 = tango.Attribute([obj.khdev '/Current']).set;
case 'singlekicker'
obj.khdevCurr0 = 0;
otherwise
error('not appropriate kickmode, either fullbump or singlekicker');
end
obj.tolhk=1; % A
obj.diithres=0.01; % loss rate %
......@@ -89,7 +106,7 @@ classdef DynamicAperture < RingControl %& GridScan
otherwise
obj.kicker_amplitude_step =10;
obj.amplhk =[100:obj.kicker_amplitude_step:1500]; % initial kicker amplitudes to scan
obj.amplhk =[450:obj.kicker_amplitude_step:1500]; % initial kicker amplitudes to scan
end
obj.amplhk0 = obj.amplhk; % initial kicker amplitudes to scan
......@@ -132,7 +149,9 @@ classdef DynamicAperture < RingControl %& GridScan
disp(['delta frequency steps [Hz]: ' num2str(obj.dfreq)]);
disp(['Scraper positions [mm]: ' num2str(obj.scraper_positions)]);
disp(['measurement type: ' obj.measurement_type]);
disp(['Kicker used: ' obj.khdev]);
disp(['Kicker used: ' obj.khdev ' (' num2str(obj.khdevCurr0) ' A)']);
disp(['Ext. Kicker: ' obj.kedev]);
disp(['Kick mode: ' obj.kickmode]);
disp(['Kicker amplitudes: ' num2str(obj.amplhk)]);
disp(['---- ---- ---- ---- ---- ----']);
......@@ -172,32 +191,58 @@ classdef DynamicAperture < RingControl %& GridScan
k1dev = tango.Device(obj.khdev);
kEdev = tango.Device('sy/ps-ke/1');
k1dev.Current=val(1);
kEdev = tango.Device(obj.kedev);
switch obj.kickmode
case 'fullbump'
k1dev.Current= obj.khdevCurr0 -val(1);
case 'singlekicker'
k1dev.Current = val(1);
otherwise
error('not appropriate kickmode, either fullbump or singlekicker');
end
k1dev.Standby();
pause(3);
InitCounterMode = kEdev.CounterMode;
kEdev.CounterMode=int16(nave);
for i=1:numel(val)
fprintf('Kicking with %.3f A\n',val(i));
k1dev.Current=val(i);
switch obj.kickmode
case 'fullbump'
fprintf('Kicking with %.3f A\n',obj.khdevCurr0 -val(i));
k1dev.Current= obj.khdevCurr0 -val(i);
case 'singlekicker'
fprintf('Kicking with %.3f A\n',val(i));
k1dev.Current=val(i);
otherwise
error('not appropriate kickmode, either fullbump or singlekicker');
end
pause(3);
r1=obj.stored_current(); % mA
kEdev.On()
pause(8+nave);
kEdev.On(); % trgger kickers
pause(8+nave);
r2=obj.stored_current(); % mA
deltai=(r2-r1)/nave/r1;
deltai=(r2-r1)/nave/r1; % average relative lost current per shot
r=-deltai*sum(1./(1+(0:nave-1)*deltai))/nave;
cur(i)=r2;
dii(i)=r;
set(hh,'YData',dii);
% stop if current reduced more than threshold (avoid kill)
if (r > 1.2*threshold) break; end
% stop if more than 90% of the initial current is lost
if (r2 < i0/10) break; end
end
% kEdev.CounterMode=int16(0);
% kEdev.CounterMode=InitCounterMode;
kEdev.Standby();
......@@ -319,8 +364,12 @@ classdef DynamicAperture < RingControl %& GridScan
obj.scraper.set(scrappos,'wait'); % movable class takes care of waiting for setting to be reached
disp(['Scraper at : ' num2str(scrappos) ' mm']);
if ip==1
pause(60);
else
pause(10);
end
disp(['scan horizontal kicker']);
try
......@@ -369,12 +418,16 @@ classdef DynamicAperture < RingControl %& GridScan
disp(['Saved data in : ' [obj.dirname '/' 'scrappos' num2str(ip) '_df' num2str(ifreq) '.mat']])
f2=figure;
plot(curs,dii,'b');hold on;
plot([x1,x2],[y1,y2]*100,'k');
plot(curthres,obj.diithres,'ro');
xlabel('h kicker amplitude [A]')
plot(curs,dii*100,'.b--','LineWidth',2,'MarkerSize',12);hold on;
plot([x1,x2],[y1,y2]*100,'g--','LineWidth',2,'MarkerSize',12);
plot(curthres,obj.diithres*100,'ro','LineWidth',2,'MarkerSize',12);
xlabel('hor. kicker amplitude [A]')
ylabel('lost current on scraper [%]');
title(['Scraper at : ' num2str(scrappos) ' mm'])
grid on;
set(gcf,'color','white');
drawnow;
saveas(gca,fullfile(obj.dirname,['scrappos' num2str(ip) '.fig']));
......
......@@ -33,13 +33,13 @@ C2=polyfit(pi(sel2),hi(sel2),0);
figure;
%subplot(2,1,1);
plot(pi,hi,'b-'); hold on;
plot(dd.position,dd.HKval,'ro');
plot(pi,polyval(P1,pi),'m-');
plot(pi,polyval(C2,pi),'c-');
xlabel('Scraper position [mm]');
ylabel('Horizontal Kicker Amplitude [A]');
title('Horizontal Kicker Amplitude to lose 5% on the scraper');
plot(pi,hi,'b-','LineWidth',2,'MarkerSize',12); hold on;
plot(dd.position,dd.HKval,'ro','LineWidth',2,'MarkerSize',12);
plot(pi,polyval(P1,pi),'m--','LineWidth',2,'MarkerSize',12);
plot(pi,polyval(C2,pi),'c--','LineWidth',2,'MarkerSize',12);
xlabel({'','Scraper position [mm]'});
ylabel({'Horizontal Kicker Amplitude [A]',''});
title({'Horizontal Kicker Amplitude';'to loose 5% on the scraper'});
%subplot(2,1,2);
%plot(pi(1:end-1),diff(hi)./diff(pi)) ;%+(pi(2)-pi(1))/2
......@@ -51,7 +51,9 @@ m2=0; b2=C2;
xintersect = (b2-b1)/(m1-m2);
yintersect = m1*xintersect + b1;
plot(xintersect,yintersect,'ks');
plot(xintersect,yintersect,'ks','LineWidth',2,'MarkerSize',12);
grid on;
set(gcf,'color','white');
try
a=load('/operation/beamdyn/matlab/optics/sr/theory/alpha');
......
......@@ -92,9 +92,9 @@ obj.initial_coordinates = movable('');
obj.ke = [TANGO_HOST 'sy/ps-ke/1'];
obj.gun = [TANGO_HOST 'elin/beam/run'];
obj.rips = [TANGO_HOST 'sy/ps-rips/manager'];
obj.scraper = [TANGO_HOST 'srdiag/scraper/c04-int'];
obj.scraper = movable([TANGO_HOST 'srdiag/scraper/c04-int/Position'],...
obj.scraper = [TANGO_HOST 'srdiag/scraper/c11-int'];
obj.scraper = movable([TANGO_HOST 'srdiag/scraper/c11-int/Position'],...
'get_fun',getfunction,'set_fun',setfunction,...
'absolute',true,...
'calibration',1,'limits',[-25 -2.5]);
......
classdef atmodel < handle
classdef atmodel < matlab.mixin.Copyable
%Class providing access to optical values
%Used as a base class for tl2.model, sy.mode, sr.model
%
......@@ -51,6 +51,7 @@ classdef atmodel < handle
emitv_
espread_
blength_
dpp_
end
properties (SetAccess=private)
......@@ -69,6 +70,7 @@ classdef atmodel < handle
end
properties (Dependent)
tunes % Tunes
dpp % Off-momentum
end
properties
......@@ -91,6 +93,10 @@ classdef atmodel < handle
end
events
LatticeModified
end
methods (Static)
function pth=getdirectory(varargin)
pth='';
......@@ -260,17 +266,13 @@ classdef atmodel < handle
end
methods (Access=protected)
function recompute(this,atstruct)
function recompute(this,varargin)
%Compute the stored parameters
if nargin >= 2
this.ring=atstruct;
this.modified=true;
end
if this.modified
ndata=length(this.ring)+1;
this.ll_=findspos(this.ring,ndata);
this.alpha_=mcf(this.ring);
[lindata,avebeta,avemu,avedisp,~,~]=atavedata(this.ring,0.0,1:ndata);
[lindata,avebeta,avemu,avedisp,~,~]=atavedata(this.ring,this.dpp_,1:ndata);
this.inidata=lindata;
tune=lindata(end).mu/2/pi;
this.tunes_=tune;
......@@ -296,7 +298,6 @@ classdef atmodel < handle
idz=this.getid('oc');
this.vlave(idz,12)=setintegfield(this.ring(idz),'PolynomB',{4});
this.modified=false;
this.emith_=[]; % Triger computation of emittances
end
function v=setintegfield(at,varargin)
......@@ -308,7 +309,7 @@ classdef atmodel < handle
function recompute_emittance(this)
if isempty(this.emith_)
[~,params]=atx(this.ring,0.0,1);
[~,params]=atx(this.ring,this.dpp_,1);
this.emith_=params.modemittance(1);
% this.emitv_=max(params.modemittance(2),5.e-12);
this.emitv_=params.modemittance(2);
......@@ -329,8 +330,8 @@ classdef atmodel < handle
if ischar(code) % string code
lcode=lower(code);
if ~isfield(this.idx,lcode)
id=this.elselect(code);
try
id=this.elselect(lcode);
this.idx.(lcode)=id;
catch
warning('atmodel:getid','Cannot store %s in the cache',code);
......@@ -345,7 +346,7 @@ classdef atmodel < handle
function idx=elselect(this,code)
%Return indices of selected elements
switch code
switch lower(code)
case 'bpm'
idx=atgetcells(this.ring,'Class','Monitor');
case 'qp'
......@@ -371,12 +372,14 @@ classdef atmodel < handle
end
methods
function this=atmodel(atstruct)
function this=atmodel(atstruct,varargin)
%Object giving access to the optical parameters
%
%ATMODEL(ATring)
%
%ATring: AT structure
% [dpp,varargs]=getargs(varargin,0); %#ok<ASGLU>
dpp=0;
ringpar=atgetcells(atstruct,'Class','RingParam');
if any(ringpar)
ringparam=atstruct{find(ringpar,1)};
......@@ -386,9 +389,23 @@ classdef atmodel < handle
this.strname='?';
this.periods=1;
end
this.recompute(atstruct);
this.ring=atstruct;
this.dpp_=dpp;
end
function set.modified(this,val)
this.modified=val;
if val
this.emith_=[]; %#ok<MCSUP>
notify(this,'LatticeModified');
end
end
function set.ring(this,atstruct)
this.ring=atstruct;
this.idx=[]; %#ok<MCSUP> % Clear cache
this.modified=true; %#ok<MCSUP> % Force recomputation
end
% dependent properties computed by atlinopt
function ll=get.ll(this)
this.recompute();
......@@ -405,6 +422,13 @@ classdef atmodel < handle
function set.tunes(this,tunes)
this.settune(tunes/this.periods);
end
function dpp=get.dpp(this)
dpp=this.dpp_;
end
function set.dpp(this,dpp)
this.dpp_=dpp;
this.modified=true;
end
function nuv=get.nuv(this)
this.recompute();
nuv=this.periods*this.tunes_(2);
......@@ -486,6 +510,17 @@ classdef atmodel < handle
'UniformOutput',false);
end
function varargout=avedata(this,varargin)
%Return average optical parameters at selected locations
%
%[lindata1,lindata2,...]=LINDATA(CODE_1,CODE_2...)
%CODEn: One of the available locations, or AT index
%lindatan: average data at the entrance of selected elements
this.recompute();
varargout=cellfun(@(arg) this.vlave(this.getid(arg),:), varargin,...
'UniformOutput',false);
end
function varargout=plot(this,varargin)
%Plot the optical functions
%
......@@ -501,10 +536,12 @@ classdef atmodel < handle
end
end
function settune(this,tunes,varargin) %#ok<INUSD>
%Change tunes
function settune(this,varargin) %#ok<INUSD>
%SETTUNE Change thetunes of the model
%
%SETTUNE(TUNES)
%SETTUNE(DPP,TUNES)
% DPP: deltaP/P, default 0
% TUNES: 1x2 vector of desired tunes
error('atmodel:NotImplemented','No generic imlementation of retuning');
end
......@@ -530,8 +567,6 @@ classdef atmodel < handle
% the strength of all quadrupoles
idz=this.getid(code);
this.ring(idz)=atsetfieldvalues(this.ring(idz),varargin{:});
this.modified=true;
this.emith_=[];
end
function varargout=getfieldvalue(this,code,varargin)
......@@ -549,6 +584,25 @@ classdef atmodel < handle
varargout{1}=atgetfieldvalues(this.ring(idz),varargin{:});
end
function substitute(this,code,modfunc,varargin)
%Substitute selected element with new ones
%
%SUBSTITUTE(CODE,MOD_FUNC,VALUE1,VALUE2,...)
% CODE: One of the available locations, or AT index
% MOD_FUNC: Function handle. Selected elements are replaced by
% NEWELEM=MODIFIER(ELEM,VALUE1,VALUE2,...)
% VALUEn: Column cell array of argument n to the MODIFIER function.
% Size should be N_SELECTED x 1
%
%Example:
%>> SUBSTITUTE('qp',@set_strength,STRENGTHS) apply
% set_strength(strength) to each quadrupole:
%
idz=this.getid(code);
this.ring(idz)=cellfun(modfunc,this.ring(idz),varargin{:},...
'UniformOutput',false);
end
function data=bpm(this,varargin)
data=this.getave('bpm',varargin{:});
end
......
......@@ -48,7 +48,7 @@ defaultNIterations=4;
defaultNoiseLevelEvaluation=false;
defaultFixedChrom=false;
expectedSextupoleGroups={'sextcor','Sext3Fams','SF2cor','Debug_2Var'};
expectedSextupoleGroups={'sextcor','Sext3Fams','SF2cor','SF2cor_fast','Debug_2Var'};
expectedSextupoleGroupsSave={'sextcor'};
expectedObjectives={'TLT','BLM','TLT_nonorm','BLM_nonorm'};
......@@ -87,7 +87,7 @@ switch TLTorBLM
% error('update averaging time and noiseval for NewBLM reading')
%averagingtime=5;
objfun=@(xx)objBLM_nonorm(xx,devopt,averagingtime);
noiseval=4000; % noise evaluated on 50 measurements with this averaging time.
noiseval=400; % noise evaluated on 50 measurements with this averaging time.
case 'TLT_nonorm'
%averagingtime=5;
objfun=@(xx)objTouLT_nonorm(xx,devopt,averagingtime);
......@@ -183,7 +183,7 @@ else
dmat = eye(Nvar);
end
step = 0.01;
step = 0.1;
g_cnt = 0;
[x1,f1,nf]=powellmain(objfun,x0,step,dmat,0,Nit,'noplot',1000*Nit);
......
......@@ -86,15 +86,16 @@ while it < maxIt
% set(0,'CurrentFigure',gcf);
% end
% if ~isempty(g_data)
[~,~,~] = process_scandata(g_data,Nvar,vrange,'plot');
% end
%direction with largest decrease
if (fm - f1)*(1) > del,
del=(fm - f1)*(1);
k=ii;
fprintf('iteration %d, var %d: del = %f updated\n',it, ii, del);
% if ~isempty(g_data)
[~,~,~] = process_scandata(g_data,Nvar,vrange,'plot');
% end
end
fm=f1;
......
......@@ -32,7 +32,7 @@ end
%% evaluate the objective function:
% get the current values. if the change is too large, do it in steps
pact=cellfun(@(d)gettangoval(d),attributenames(devopt))';
maxchange=0.1;
maxchange=2;
if max(abs(p-pact))>maxchange
......@@ -67,7 +67,7 @@ else
end
% wait power supply
pause(2);
pause(.5);
obj=getTotalLosses(averagingtime);
......
......@@ -2,13 +2,13 @@ function [names,vrange,cellnamerange]=attributenames(typedevice)
% [names,range]=attributenames(typedevice)
%
% list attribute names and possible ranges for typedevice:
% sextcor, SF2cor, skewcor,
% sextcor, SF2cor, skewcor,
%
%see also:
switch typedevice
case 'injsextcor'
case 'injsextcor'
% +/- 2 A
sxt = tango.Attribute('srmag/m-s/all/Strengths');
......@@ -30,7 +30,7 @@ switch typedevice
cellnamerange=magnames;
case 'sextcor'
% +/- 2 A
......@@ -51,36 +51,39 @@ switch typedevice
% +/- 2 A
sxt = tango.Attribute('srmag/m-s/all/Strengths');
names = sxt.device.MagnetNames.read;
names = sxt.device.MagnetNames.read;
K = sxt.value;
magnames={};
indsxt = 1:length(names);
for is = sort(indsxt([2:6:end,5:6:end]))
magnames = [magnames,{[names{is} '/Strength'];[-0.2 +0.2]+K(is)}];
magnames = [magnames,{[names{is} '/Strength'];[-0.3 +0.3]+K(is)}];
end
cellnamerange=magnames;
case 'Debug_2Var'
case 'Debug_2Var'
% +/- 2 A
sxt = tango.Attribute('srmag/m-s/all/Strengths');
names = sxt.device.MagnetNames.read;
names = sxt.device.MagnetNames.read;
K = sxt.value;
magnames={};
indsxt = 1:length(names);
for is = sort(indsxt([2,191]))
magnames = [magnames,{[names{is} '/Strength'];[-0.2 +0.2]+K(is)}];
magnames = [magnames,{[names{is} '/Strength'];[-0.5 +0.5]+K(is)}];
end
cellnamerange=magnames;
case 'Sext3Fams'
case 'SF2cor_fast'
load('/mntdirect/_machfs/MDT/2020/2020_06_28/SextupolesOptimizations/SF2cor_names.mat')
case 'Sext3Fams'
% +/- 2 A
......@@ -88,7 +91,7 @@ switch typedevice
'srmag/m-sf2/all/DeltaCorrectionStrength',...
'srmag/m-sd1bd/all/DeltaCorrectionStrength'};
K=zeros(size(names));
magnames={};
indsxt = 1:length(names);
......@@ -103,18 +106,18 @@ switch typedevice
case 'quadcor'
% +/- 2 A
% +/- 2 A
% +/- 2 A
sxt = tango.Attribute('srmag/m-q/all/Strengths');
names = sxt.device.MagnetNames.read;
names = sxt.device.MagnetNames.read;
K = sxt.value;
magnames={};
indsxt = 1:length(names);
for is = indsxt
if strcmp(names{is}(1:9),'srmag/m-q')
magnames = [magnames,{[names{is} '/Strength'];[0.98 1.02]*K(is)}];
magnames = [magnames,{[names{is} '/Strength'];[0.98 1.02]*K(is)}];
end
end
......@@ -127,7 +130,7 @@ switch typedevice
% +/- 2 A
sxt = tango.Attribute('srmag/m-q/all/Strengths');
names = sxt.device.MagnetNames.read;
names = sxt.device.MagnetNames.read;