Commit af6f4509 authored by Simone Liuzzo's avatar Simone Liuzzo
Browse files

class for orbit correction (include absolute correction)

parent ee94b2ad
function ModelRM = MeasureTrajectoryRM(obj,varargin)
% measure trajectory RM (single sided)
% optional inputs:
% 'h',true
% 'v',true
% 'cv',false
% 'sp',false
% 'rf',false
simulator = strcmp(obj.machine,'ebs-simu');
switch_on_gun_rips = true;
if simulator
switch_on_gun_rips = false;
end
sfi = obj.savefile;
obj.savefile = false;
waitPStime = 3; % [s]
%
p = inputParser;
addRequired(p,'obj'); % default empty
addOptional(p,'h',true,@islogical); % default empty
addOptional(p,'v',true,@islogical); % default empty
addOptional(p,'cv',true,@islogical); % default empty
addOptional(p,'sp',true,@islogical); % default empty
addOptional(p,'rf',true,@islogical); % default empty
parse(p,obj,varargin);
h = p.Results.h;
v = p.Results.v;
cv = p.Results.cv;
sp = p.Results.sp;
measrf = p.Results.rf;
% h=false;
% v=false;
% cv=true;
% sp=true;
% deltarf = 50;
% measrf = false;
% start innjection once and for all
if switch_on_gun_rips
gun = tango.Device(obj.dev_gun);
rips = tango.Device(obj.dev_rips);
rips.StartRamping()
gun.On();
while ~strcmp(rips.State,'Running')
disp('rips moving!')
pause(1);
end
end
if ~simulator
deltak = 1e-4; % rad
else
deltak = 1e-6; % rad
end
t0 = obj.measuretrajectory;
enabledcorH = obj.getEnabledHsteerers;
enabledcorV = obj.getEnabledVsteerers;
% % for quick test
% enabledcorH(4:end)=false;
% enabledcorV(4:end)=false;
% H steerers RM
ModelRM.TrajHCor{1} = [];
ModelRM.TrajHCor{3} = [];
ModelRM.TrajVCor{1} = [];
ModelRM.TrajVCor{3} = [];
if h
hs0 = obj.getHSteerers;
for ik = 1:length(hs0) %
if enabledcorH(ik) % only if eneabled corrector
disp(['H steerer' num2str(ik)])
hs = hs0;
hs(ik) = hs(ik) + deltak;
obj.setHSteerers(hs);
pause(waitPStime);
t = obj.measuretrajectory;
obj.setHSteerers(hs0);
rmh(:,ik) = (t(1,:) - t0(1,:))./deltak; % m/rad
rmv(:,ik) = (t(2,:) - t0(2,:))./deltak;
else
disp(['H steerer' num2str(ik) ' is disabled'])
rmh(:,ik) = NaN(size(t0(1,:)));
rmv(:,ik) = NaN(size(t0(2,:)));
end
end
ModelRM.TrajHCor{1} = rmh;
ModelRM.TrajHCor{3} = rmv;
if measrf
% RF RM
disp('RF RM');
rf0 = obj.getRFFrequency;
rf = rf0 + deltarf;
obj.setRFFrequency(rf);
pause(waitPStime);
t = obj.measuretrajectory;
obj.setRFFrequency(rf0);
ModelRM.TrajHDPP = (t(1,:) - t0(1,:))./deltak;
ModelRM.TrajVDPP = (t(2,:) - t0(2,:))./deltak;
else
ModelRM.TrajHDPP = [];
ModelRM.TrajVDPP = [];
end
end
obj.ModelRM = ModelRM; %temporary store of ModelRM
% V steerers RM
if v
vs0 = obj.getVSteerers;
for ik = 1:length(vs0) %1:2%
if enabledcorV(ik) % only if eneabled corrector
disp(['V steerer' num2str(ik) ' is disabled'])
vs = vs0;
vs(ik) = vs(ik) + deltak;
obj.setVSteerers(vs);
pause(waitPStime);
t = obj.measuretrajectory;
obj.setVSteerers(vs0);
rmh(:,ik) = (t(1,:) - t0(1,:))./deltak;
rmv(:,ik) = (t(2,:) - t0(2,:))./deltak;
else
disp(['V steerer' num2str(ik) ' is disabled'])
rmh(:,ik) = NaN(size(t0(1,:)));
rmv(:,ik) = NaN(size(t0(2,:)));
end
end
ModelRM.TrajVCor{1} = rmh;
ModelRM.TrajVCor{3} = rmv;
end
obj.ModelRM = ModelRM; %temporary store of ModelRM
% septa
disp('septa RM')
if sp
hs0 = obj.getSepta;
for ik = 1:length(hs0)
hs = hs0;
hs(ik) = hs(ik) + deltak;
obj.setSepta(hs);
pause(waitPStime);
t = obj.measuretrajectory;
obj.setSepta(hs0);
ModelRM.RMSP(ik,:) = (t(1,:) - t0(1,:))./deltak;
end
end
obj.ModelRM = ModelRM; %temporary store of ModelRM
% CV TL2 steerers RM
disp('CV RM')
if cv
vs0 = obj.getTL2CV;
for ik = 1:length(vs0)
vs = vs0;
vs(ik) = vs(ik) + deltak;
obj.setTL2CV(vs);
pause(waitPStime);
t = obj.measuretrajectory;
obj.setTL2CV(vs0);
ModelRM.RMCV(ik,:) = (t(2,:) - t0(2,:))./deltak;
end
end
obj.ModelRM = ModelRM; %temporary store of ModelRM
if switch_on_gun_rips
pause(0.25); % rips still moving.
rips.StopRamping()
gun.Off(); % limit dose
%Ke.PulseNumber = pn; % set back initial number of pulses for extraction
end
obj.savefile = sfi;
obj.ModelRM = ModelRM;
end
classdef OrbitCorrection < RingControl
%OrbitCorrection
%
%see also: RingControl
properties
ModelRM % response matrix structure for trajectory correction
end
methods
% creator
function obj = OrbitCorrection(machine)
%OrbitCorrection Construct an instance of this class
%
% can be used for: ebs-simu, esrf-sr, esrf-sy
%
% Inherits all methods and properties from RingControl:
%
%see also: RingControl
obj@RingControl(machine);
obj.savefile = true;
obj.n_acquisitions = 1;
switch obj.machine
case 'ebs-simu'
b = load('/machfs/liuzzo/EBS/Simulator/commissioningtools/@OrbitCorrection/EBSOrbitRM.mat');
obj.ModelRM = b.ModelRM;
case 'esrf-sr'
case 'esrf-sy'
end
end
% methods defined in seprated functions
% measure response matrix
modelRM = MeasureOrbitRM(obj);
% compute response matrix
modelRM = SimulateOrbitRM(obj,varargin);
% compute RM
C = findrespmat(obj, PERTURB, PVALUE,FIELD,M,N, varargin)
% compute closed orbit with BPM errs
[t0] = findorbit(obj,ring,incod);
% correction function
[hs,vs,cv,sp,rf] = correction(obj,varargin);
end
methods (Static)
% test function
tryme
end
end
function ModelRM = SimulateOrbitRM(obj, varargin)
% SimulateOrbitRM computes a orbit response matrix for lattice
%
% INPUT (variable):
% 'indBPM' : AT bpm indexes (default: from obj.rmodel, Monitor)
% 'indHCor' : AT hor. correctors indexes (default: from obj.rmodel, KickAngle)
% 'indVCor' : AT ver. correctors indexes (default: from obj.rmodel, KickAngle)
% 'kick' : [rad] kick for RM computation
% 'delta' : [] energy variation for frequency change RM computation
% 'RF' : use RF (logic, default: true)
% 'rad' : use radiation (logic, default: true)
%
%
% OUTPUT:
% ModelRM structure with fields:
% struct with fields (640 BPM, 288 correctors in this case):
%
% OrbHCor: {[640×288 double] [] [640×288 double]}
% OrbVCor: {[640×288 double] [] [640×288 double]}
% OrbHDPP: [1×640 double]
% OrbVDPP: [1×640 double]
% kval: scalar, rad, kick for RM
% delta: scalar,[], energy deviation for RM
%
%
% ex:
% >> ft = OrbitCorrection('ebs-simu');
% >> RM = ft.SimulateOrbitRM('RF',false,'rad',true); % simulate trajectory RM with RF off.
% >> ft.ModelRM % display computed trajectory RMs
%
%
%
%see also: OrbitCorrection.findorbit OrbitCorrection.findrespmat
r0 = obj.rmodel; % ring lattice
%% parse inputs
p = inputParser();
addRequired(p,'obj');
addOptional(p,'indBPM',obj.indBPM,@isnumeric);
addOptional(p,'indHCor',obj.indHst,@isnumeric);
addOptional(p,'indVCor',obj.indVst,@isnumeric);
addOptional(p,'kick',1e-4,@isnumeric);
addOptional(p,'delta',1e-3,@isnumeric);
addOptional(p,'RF',true,@islogical);
addOptional(p,'rad',true,@islogical);
parse(p,obj,varargin{:});
%alpha = mcf(r0);
%f0 = obj.getRFFrequency;% atgetfieldvalues(r0,find(atgetcells(r0,'Frequency'),1,'first'),'Frequency');
inCOD = zeros(6,1);
indBPM = p.Results.indBPM;
indHorCor = p.Results.indHCor;
indVerCor = p.Results.indVCor;
kval = p.Results.kick;
delta = p.Results.delta;
rf_on = p.Results.RF;
radon = p.Results.rad;
indrfc = findcells(r0,'Frequency');
hn=r0{indrfc(1)}.HarmNumber;
Vt=sum(atgetfieldvalues(r0,indrfc,'Voltage'));
if ~rf_on & ~radon
warning('no rf and no radiation is not a realistic lattice condition')
r0 = atsetcavity(r0,0.0,0.0,hn);
elseif ~rf_on & radon
r0 = atsetcavity(r0,0.0,1.0,hn);
elseif rf_on & ~radon
warning('RF on but no radiation is not a realistic lattice condition')
r0 = atsetcavity(r0,Vt,0.0,hn);
elseif rf_on & radon
r0 = atsetcavity(r0,Vt,1.0,hn);
end
ModelRM.kval = kval;
ModelRM.delta = delta;
% % for quick test
% indHorCor(4:end)=[];
% indVerCor(4:end)=[];
K = [0,0];
CV = [0,0];
Septa = [0,0];
disp('Simulated orbit RM computation in progress ... ')
[t0] = obj.findorbit(r0);
%% horizontal
% trajectory
RMTH= obj.findrespmat(indHorCor,kval,'KickAngle',1,1);
RMTH{1}=RMTH{1}./kval;
RMTH{2}=RMTH{2}./kval;
RMTH{3}=RMTH{3}./kval;
RMTH{4}=RMTH{4}./kval;
disp(' --- computed horizontal orbit RM from model --- ');
% store data
ModelRM.OrbHCor=RMTH;
% temporary store of ModelRM
obj.ModelRM = ModelRM;
%% vertical
RMTV = obj.findrespmat(indVerCor,kval,'KickAngle',1,2);
RMTV{1}=RMTV{1}./kval;
RMTV{2}=RMTV{2}./kval;
RMTV{3}=RMTV{3}./kval;
RMTV{4}=RMTV{4}./kval;
disp(' --- computed vertical orbit RM from model --- ');
% store data
ModelRM.OrbVCor=RMTV;
% temporary store of ModelRM
obj.ModelRM = ModelRM;
%% dpp
% orbit RM dpp
inCODPdpp=inCOD;
inCODPdpp(5)=delta;% linepass is used.
% o=trajfun(r0,indBPM,inCODPdpp);
o=obj.findorbit(r0,inCODPdpp);
oxPdpp=o(1,:);
oyPdpp=o(3,:);
inCODMdpp=inCOD;
inCODMdpp(5)=-delta;
% o=trajfun(r0,indBPM,inCODMdpp);
o=obj.findorbit(r0,inCODMdpp);
oxMdpp=o(1,:);
oyMdpp=o(3,:);
dppH=(oxPdpp-oxMdpp)/2/delta;
dppV=(oyPdpp-oyMdpp)/2/delta;
disp(' --- computed orbit/dpp RM from model --- ');
% store data
ModelRM.OrbHDPP=dppH;
ModelRM.OrbVDPP=dppV;
% store simulated ModelRM
obj.ModelRM = ModelRM;
disp('A new response matrix has been set for correction');
end
function ModelRM = SimulateOrbitRM(obj, varargin)
% SimulateOrbitRM computes a trajectory response matrix for lattice
%
% INPUT (variable):
% 'indBPM' : AT bpm indexes (default: from obj.rmodel, Monitor)
% 'indHCor' : AT hor. correctors indexes (default: from obj.rmodel, KickAngle)
% 'indVCor' : AT ver. correctors indexes (default: from obj.rmodel, KickAngle)
% 'kick' : [rad] kick for RM computation
% 'delta' : [] energy variation for frequency change RM computation
% 'RF' : use RF (logic, default: true)
% 'rad' : use radiation (logic, default: true)
%
%
% OUTPUT:
% ModelRM structure with fields:
% struct with fields (640 BPM, 288 correctors in this case):
%
% TrajHCor: {[640×288 double] [] [640×288 double]}
% TrajVCor: {[640×288 double] [] [640×288 double]}
% TrajHDPP: [1×640 double]
% TrajVDPP: [1×640 double]
% RMSP: [2×640 double]
% RMCV: [2×640 double]
% kval: scalar, rad, kick for RM
% delta: scalar,[], energy deviation for RM
%
%
% ex:
% >> ft = FirstTurns('ebs-simu');
% >> ft.setTurns(4); % set turns to compute trajetory RM
% >> ft.injoff = -6e-3; % set the injected - stored beam offset.
% >> TrajRM = ft.SimulateTrajectoryRM('RF',false,'rad',true); % simulate trajectory RM with RF off.
% >> ft.ModelRM % display computed trajectory RMs
%
%
%
%see also: First2Turns.find2turnstrajectory6Err First2Turns.findrespmat2turns
r0 = obj.rmodel; % ring lattice
%% parse inputs
p = inputParser();
addRequired(p,'obj');
addOptional(p,'indBPM',obj.indBPM,@isnumeric);
addOptional(p,'indHCor',obj.indHst,@isnumeric);
addOptional(p,'indVCor',obj.indVst,@isnumeric);
addOptional(p,'kick',1e-4,@isnumeric);
addOptional(p,'delta',1e-3,@isnumeric);
addOptional(p,'RF',true,@islogical);
addOptional(p,'rad',true,@islogical);
parse(p,obj,varargin{:});
%alpha = mcf(r0);
%f0 = obj.getRFFrequency;% atgetfieldvalues(r0,find(atgetcells(r0,'Frequency'),1,'first'),'Frequency');
inCOD = zeros(6,1);
indBPM = p.Results.indBPM;
indHorCor = p.Results.indHCor;
indVerCor = p.Results.indVCor;
kval = p.Results.kick;
delta = p.Results.delta;
rf_on = p.Results.RF;
radon = p.Results.rad;
indrfc = findcells(r0,'Frequency');
hn=r0{indrfc(1)}.HarmNumber;
Vt=sum(atgetfieldvalues(r0,indrfc,'Voltage'));
if ~rf_on & ~radon
warning('no rf and no radiation is not a realistic lattice condition')
r0 = atsetcavity(r0,0.0,0.0,hn);
elseif ~rf_on & radon
r0 = atsetcavity(r0,0.0,1.0,hn);
elseif rf_on & ~radon
warning('RF on but no radiation is not a realistic lattice condition')
r0 = atsetcavity(r0,Vt,0.0,hn);
elseif rf_on & radon
r0 = atsetcavity(r0,Vt,1.0,hn);
end
ModelRM.kval = kval;
ModelRM.delta = delta;
% % for quick test
% indHorCor(4:end)=[];
% indVerCor(4:end)=[];
K = [0,0];
CV = [0,0];
Septa = [0,0];
disp('Simulated Trajectory RM computation in progress ... ')
[t0] = obj.findNturnstrajectory6Err(r0,K,CV,Septa);
%% horizontal
% trajectory
RMTH= obj.findrespmatNturns(indHorCor,kval,'KickAngle',1,1);
RMTH{1}=RMTH{1}./kval;
RMTH{2}=RMTH{2}./kval;
RMTH{3}=RMTH{3}./kval;
RMTH{4}=RMTH{4}./kval;
disp(' --- computed horizontal trajectory RM from model --- ');
% store data
ModelRM.TrajHCor=RMTH;
% temporary store of ModelRM
obj.ModelRM = ModelRM;
%% vertical
RMTV = obj.findrespmatNturns(indVerCor,kval,'KickAngle',1,2);
RMTV{1}=RMTV{1}./kval;
RMTV{2}=RMTV{2}./kval;
RMTV{3}=RMTV{3}./kval;
RMTV{4}=RMTV{4}./kval;
disp(' --- computed vertical trajectory RM from model --- ');
% store data
ModelRM.TrajVCor=RMTV;
% temporary store of ModelRM
obj.ModelRM = ModelRM;
%% dpp
% orbit RM dpp
inCODPdpp=inCOD;
inCODPdpp(5)=delta;% linepass is used.
% o=trajfun(r0,indBPM,inCODPdpp);
o=obj.findNturnstrajectory6Err(r0,K,CV,Septa,inCODPdpp);
oxPdpp=o(1,:);
oyPdpp=o(3,:);
inCODMdpp=inCOD;
inCODMdpp(5)=-delta;
% o=trajfun(r0,indBPM,inCODMdpp);
o=obj.findNturnstrajectory6Err(r0,K,CV,Septa,inCODMdpp);
oxMdpp=o(1,:);
oyMdpp=o(3,:);
dppH=(oxPdpp-oxMdpp)/2/delta;
dppV=(oyPdpp-oyMdpp)/2/delta;
disp(' --- computed trajectory/dpp RM from model --- ');
% store data
ModelRM.TrajHDPP=dppH;
ModelRM.TrajVDPP=dppV;
% temporary store of ModelRM
obj.ModelRM = ModelRM;
%% TL2 septa correctors RM
obj.ModelRM.RMCV=[];
obj.ModelRM.RMSP=[];
cv8rm = obj.findNturnstrajectory6Err(r0,K,CV+[kval 0],Septa)-t0;