Skip to content

Matlab frequency control #560

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Apr 5, 2023
Merged

Matlab frequency control #560

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
0e88b32
Handle frequency change of dp optics in 6d
lfarv Feb 9, 2023
30202cd
trajectory response
lfarv Feb 14, 2023
a948db3
naming
lfarv Feb 17, 2023
44b962a
Added a test for off-momentum in 6D lattices
lfarv Mar 3, 2023
9fa3a87
Add tests off-momentum with dct and df
lfarv Apr 5, 2023
9e16dd1
Help improved
lfarv Apr 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 7 additions & 6 deletions atmat/atphysics/LinearOptics/atlinopt6.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,9 +80,9 @@
%
% See also atlinopt2 atlinopt4 tunechrom

[ringdata,elemdata]=wrapper6d(ring,@xlinopt6,varargin{:});
[ringdata,elemdata]=frequency_control(@xlinopt6,ring,varargin{:});

function [ringdata,elemdata] = xlinopt6(ring, is6d, varargin)
function [ringdata,elemdata] = xlinopt6(ring, varargin)
clight = PhysConstant.speed_of_light_in_vacuum.value; % m/s
[get_chrom, varargs]=getflag(varargin, 'get_chrom');
[get_w, varargs]=getflag(varargs, 'get_w');
Expand All @@ -91,6 +91,7 @@
[DPStep,~]=getoption(varargs,'DPStep');
[cavargs,varargs]=getoption(varargs,{'cavpts'});
[refpts,varargs]=getargs(varargs,1,'check',@(arg) isnumeric(arg) || islogical(arg));
is_6d=getoption(varargs,'is_6d',[]); % Always set by frequency_control, keep in varargs

if isempty(twiss_in) % Circular machine
[orbs,orbitin]=findorbit(ring,refpts,dpargs{:},varargs{:});
Expand Down Expand Up @@ -120,15 +121,15 @@
)';
ringdata=struct('tune',tunes,'damping_time',damping_times);

if is6d % 6D processing
if is_6d % 6D processing
[alpha,beta,disp]=cellfun(@output6,ri,'UniformOutput',false);
if get_w || get_chrom
frf=get_rf_frequency(ring);
DFStep=-DPStep*mcf(atradoff(ring))*frf;
rgup=atsetcavity(ring,'Frequency',frf+0.5*DFStep,cavargs{:});
rgdn=atsetcavity(ring,'Frequency',frf-0.5*DFStep,cavargs{:});
[~,o1P]=findorbit(rgup,'guess',orbitin,varargs{:});
[~,o1M]=findorbit(rgdn,'guess',orbitin,varargs{:});
[~,o1P]=findorbit6(rgup,'guess',orbitin,varargs{:});
[~,o1M]=findorbit6(rgdn,'guess',orbitin,varargs{:});
if get_w
[ringdata.chromaticity,w]=chrom_w(rgup,rgdn,o1P,o1M,refpts);
[elemdata.W]=deal(w{:});
Expand Down Expand Up @@ -180,7 +181,7 @@

function [vals,ms,phis,rmats,as]=build_1turn_map(ring,refpts,orbit,varargin)
% Build the initial distribution at entrance of the transfer line
if is6d
if is_6d
[mt,ms]=findm66(ring,refpts,'orbit',orbit,varargs{:});
else
[mt,ms]=findm44(ring,NaN,refpts,'orbit',orbit,varargs{:});
Expand Down
1 change: 1 addition & 0 deletions atmat/atphysics/LinearOptics/findm44.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@
varargs=getdparg(varargs);
[dp,varargs]=getoption(varargs,'dp',0.0);
[dpargs,varargs]=getoption(varargs,{'dct','df'});
[~,varargs]=getoption(varargs,'is_6d',[]); % Consume the is_6d option
[refpts,orbitin,varargs]=getargs(varargs,[],orbitin,'check',@(x) ~(ischar(x) || isstring(x))); %#ok<ASGLU>

if islogical(refpts)
Expand Down
98 changes: 56 additions & 42 deletions atmat/atphysics/LinearOptics/findm66.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [M66, varargout] = findm66(LATTICE, varargin)
function varargout = findm66(LATTICE, varargin)
%FINDM66 numerically finds the 6x6 transfer matrix of an accelerator lattice
% by differentiation of LINEPASS near the closed orbit
% FINDM66 uses FINDORBIT6 to search for the closed orbit in 6-d
Expand All @@ -7,6 +7,18 @@
% M66 = FINDM66(RING) finds the full one-turn 6-by-6
% matrix at the entrance of the first element
%
%[...]=FINDM66(RING,...,'dp',DP) Specify the momentum deviation when
% radiation is OFF (default: 0)
%
%[...]=FINDM66(RING,...,'dct',DCT) Specify the path lengthening when
% radiation is OFF (default: 0)
%
%[...]=FINDM66(RING,...,'df',DF) Specify the RF frequency deviation
% radiation is OFF (default: 0)
%
%[...]=FINDM66(RING,...,'orbit',ORBIT) Specify the orbit at the entrance
% of the ring, if known.
%
% [M66,T] = FINDM66(RING,REFPTS) in addition to M finds
% 6-by-6 transfer matrixes between entrances of
% the first element and each element indexed by REFPTS.
Expand Down Expand Up @@ -37,52 +49,54 @@
if ~iscell(LATTICE)
error('First argument must be a cell array');
end
NE = length(LATTICE);
[XYStep,varargs]=getoption(varargin,'XYStep'); % Step size for numerical differentiation %1.e-8
[DPStep,varargs]=getoption(varargs,'DPStep'); % Step size for numerical differentiation %1.e-6
[orbitin,varargs]=getoption(varargs,'orbit',[]);
[refpts,orbitin,varargs]=getargs(varargs,[],orbitin,'check',@(x) ~(ischar(x) || isstring(x))); %#ok<ASGLU>
[varargout{1:nargout}] = frequency_control(@xfindm66,LATTICE,varargin{:});

if islogical(refpts)
refpts(end+1:NE+1)=false;
elseif isnumeric(refpts)
refpts=setelems(false(1,NE+1),refpts);
else
error('REFPTS must be numeric or logical');
end
function [M66, varargout] = xfindm66(LATTICE, varargin)
NE = length(LATTICE);
[XYStep,varargs]=getoption(varargin,'XYStep'); % Step size for numerical differentiation %1.e-8
[DPStep,varargs]=getoption(varargs,'DPStep'); % Step size for numerical differentiation %1.e-6
[orbitin,varargs]=getoption(varargs,'orbit',[]);
[dpargs,varargs]=getoption(varargs,{'dp','dct','df'});
[is_6d,varargs]=getoption(varargs,'is_6d',[]); % Always set by frequency_control
[refpts,orbitin,varargs]=getargs(varargs,[],orbitin,'check',@(x) ~(ischar(x) || isstring(x))); %#ok<ASGLU>

if isempty(orbitin)
if check_radiation(LATTICE)
orbitin = findorbit6(LATTICE,'XYStep',XYStep,'DPStep',DPStep);
else
[~, orbitin] = findorbit4(LATTICE,0.0,'XYStep',XYStep);
end
end
if islogical(refpts)
refpts(end+1:NE+1)=false;
elseif isnumeric(refpts)
refpts=setelems(false(1,NE+1),refpts);
else
error('REFPTS must be numeric or logical');
end

refs=setelems(refpts,NE+1);
reqs=refpts(refs);
if isempty(orbitin)
[~, orbitin] = findorbit(LATTICE,'XYStep',XYStep,'DPStep',DPStep,'is_6d',is_6d,dpargs{:});
end

% Build a diagonal matrix of initial conditions
%scaling=2*XYStep*[1 0.1 1 0.1 1 1];
scaling=XYStep*[1 1 1 1 0 0] + DPStep*[0 0 0 0 1 1];
D6 = 0.5*diag(scaling);
% Add to the orbit_in. First 12 columns for derivative
% 13-th column is for closed orbit
RIN = orbitin + [D6 -D6 zeros(6,1)];
ROUT = linepass(LATTICE,RIN,refs);
TMAT3 = reshape(ROUT,6,13,[]);
M66 = (TMAT3(:,1:6,end)-TMAT3(:,7:12,end))./scaling;
refs=setelems(refpts,NE+1);
reqs=refpts(refs);

if nargout >= 2 % Calculate matrices at all REFPTS.
varargout{1} = (TMAT3(:,1:6,reqs)-TMAT3(:,7:12,reqs))./scaling;
% Return closed orbit if requested
if nargout >= 3
varargout{2}=squeeze(TMAT3(:,13,reqs));
end
end
% Build a diagonal matrix of initial conditions
%scaling=2*XYStep*[1 0.1 1 0.1 1 1];
scaling=XYStep*[1 1 1 1 0 0] + DPStep*[0 0 0 0 1 1];
D6 = 0.5*diag(scaling);
% Add to the orbit_in. First 12 columns for derivative
% 13-th column is for closed orbit
RIN = orbitin + [D6 -D6 zeros(6,1)];
ROUT = linepass(LATTICE,RIN,refs);
TMAT3 = reshape(ROUT,6,13,[]);
M66 = (TMAT3(:,1:6,end)-TMAT3(:,7:12,end))./scaling;

function mask=setelems(mask,idx)
mask(idx)=true;
end
if nargout >= 2 % Calculate matrices at all REFPTS.
varargout{1} = (TMAT3(:,1:6,reqs)-TMAT3(:,7:12,reqs))./scaling;
% Return closed orbit if requested
if nargout >= 3
varargout{2}=squeeze(TMAT3(:,13,reqs));
end
end

function mask=setelems(mask,idx)
mask(idx)=true;
end

end
end
68 changes: 35 additions & 33 deletions atmat/atphysics/Orbit/findorbit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,39 +30,41 @@
%
% See also FINDORBIT4, FINDSYNCORBIT, FINDORBIT6.

[orbitin,varargs]=getoption(varargin,'orbit',[]);
[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
[dp,varargs]=getoption(varargs,'dp',NaN);
[dct,varargs]=getoption(varargs,'dct',NaN);
[df,varargs]=getoption(varargs,'df',NaN);
if isempty(orbitin)
if check_6d(ring) % Radiation ON: 6D orbit
if isfinite(dp) || isfinite(dct) || isfinite(df)
warning('AT:linopt','In 6D, "dp", "dct" and "df" are ignored');
[orbs,orbitin]=frequency_control(@xfindorbit,ring,varargin{:});

function [orbs,orb0]=xfindorbit(ring,varargin)
[orb0,varargs]=getoption(varargin,'orbit',[]);
[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
[dp,varargs]=getoption(varargs,'dp',NaN);
[dct,varargs]=getoption(varargs,'dct',NaN);
[df,varargs]=getoption(varargs,'df',NaN);
[is_6d,varargs]=getoption(varargs,'is_6d',[]); % Always set by frequency_control
if isempty(orb0)
if is_6d % Radiation ON: 6D orbit
orb0=xorbit_6(ring,varargs{:});
elseif isfinite(df)
[cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
'cell_harmnumber','cell_revolution_frequency');
dct=-cell_l*df/(cell_frev*cell_h+df);
orb0=xorbit_ct(ring,dct,varargs{:});
elseif isfinite(dct)
orb0=xorbit_ct(ring,dct,varargs{:});
else
orb0=xorbit_dp(ring,dp,varargs{:});
end
args={'KeepLattice'};
else
args={};
end
orbitin=xorbit_6(ring,varargs{:});
elseif isfinite(df)
[cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
'cell_harmnumber','cell_revolution_frequency');
dct=-cell_l*df/(cell_frev*cell_h+df);
orbitin=xorbit_ct(ring,dct,varargs{:});
elseif isfinite(dct)
orbitin=xorbit_ct(ring,dct,varargs{:});
else
orbitin=xorbit_dp(ring,dp,varargs{:});
end
args={'KeepLattice'};
else
args={};
end

if islogical(refpts)
refpts=find(refpts);
end
if isempty(refpts)
% return only the fixed point at the entrance of RING{1}
orbs=orbitin;
else
orbs=linepass(ring,orbitin,refpts,args{:});
end
if islogical(refpts)
refpts=find(refpts);
end
if isempty(refpts)
% return only the fixed point at the entrance of RING{1}
orbs=orb0;
else
orbs=linepass(ring,orb0,refpts,args{:});
end
end
end
1 change: 1 addition & 0 deletions atmat/atphysics/Orbit/findorbit4.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@
[dct,varargs]=getoption(varargs,'dct',NaN);
[df,varargs]=getoption(varargs,'df',NaN);
[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
[~,varargs]=getoption(varargs,'is_6d',[]); % Consume the is_6d option
if isempty(orbitin)
if isfinite(df)
[cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
Expand Down
50 changes: 33 additions & 17 deletions atmat/atphysics/Orbit/findorbit6.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,17 @@
% from the range 1 to length(RING)+1.
% See further explanation of REFPTS in the 'help' for FINDSPOS
%
% FINDORBIT6(...,'dp',DP)
% Specify the off-momentum. The RF frequency will be adjusted to get the
% desired value
%
% FINDORBIT6(...,'dct',DCT)
% Specify the path lengthening. The RF frequency will be adjusted to get
% the desired value
%
% FINDORBIT6(...,'df',DF)
% Specify the RF frequency deviation
%
% FINDORBIT6(RING,REFPTS,GUESS)
% FINDORBIT6(...,'guess',GUESS) The search for the fixed point
% starts from initial condition GUESS. Otherwise the search starts from
Expand All @@ -56,22 +67,27 @@
if ~iscell(ring)
error('First argument must be a cell array');
end
[orbitin,varargs]=getoption(varargin,'orbit',[]);
[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
if isempty(orbitin)
orbitin=xorbit_6(ring,varargs{:});
args={'KeepLattice'};
else
args={};
end
[orb6,orbitin] = frequency_control(@xfindorbit6,ring,varargin{:});

if islogical(refpts)
refpts=find(refpts);
end
if isempty(refpts)
% return only the fixed point at the entrance of RING{1}
orb6=orbitin;
else
orb6 = linepass(ring,orbitin,refpts,args{:});
end
function[orb6,orbitin] = xfindorbit6(ring,varargin)
[orbitin,varargs]=getoption(varargin,'orbit',[]);
[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
[~,varargs]=getoption(varargs,'is_6d',[]); %% Consume the is_6d option
if isempty(orbitin)
orbitin=xorbit_6(ring,varargs{:});
args={'KeepLattice'};
else
args={};
end

if islogical(refpts)
refpts=find(refpts);
end
if isempty(refpts)
% return only the fixed point at the entrance of RING{1}
orb6=orbitin;
else
orb6 = linepass(ring,orbitin,refpts,args{:});
end
end
end
1 change: 1 addition & 0 deletions atmat/atphysics/Orbit/findsyncorbit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@
[dp,varargs]=getoption(varargs,'dp',NaN);
[df,varargs]=getoption(varargs,'df',NaN);
[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
[~,varargs]=getoption(varargs,'is_6d',[]); % Consume the is_6d option
if isempty(orbitin)
if isfinite(df)
[cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
Expand Down
13 changes: 7 additions & 6 deletions atmat/atphysics/TuneAndChromaticity/tunechrom.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,15 +40,16 @@
if cpl1 || ~cpl2
warning('AT:ObsoleteParameter','The "coupled" flag is ignored: coupling is always assumed');
end
[varargout{1:nargout}]=wrapper6d(ring,@xtunechrom,allargs{:});
[varargout{1:nargout}]=frequency_control(@xtunechrom,ring,allargs{:});

function [tune, chrom] = xtunechrom(ring,is6d,varargin)
function [tune, chrom] = xtunechrom(ring,varargin)
[oldchrom,varargs]=getflag(varargin,'chrom');
[chrom,varargs]=getflag(varargs,'get_chrom');
[dpargs,varargs]=getoption(varargs,{'orbit','dp','dct','df'});
[cavargs,varargs]=getoption(varargs,{'cavpts'});
[DPStep,~]=getoption(varargs,'DPStep');
if is6d
is_6d=getoption(varargs,'is_6d',[]); % Always set by frequency_control, keep in varargs
if is_6d
tunefunc=@tune6;
else
tunefunc=@tune4;
Expand All @@ -59,13 +60,13 @@
tune=tunefunc(ring,'orbit',orbitin,varargs{:});

if chrom || oldchrom || nargout == 2
if is6d
if is_6d
frf=get_rf_frequency(ring);
DFStep=-DPStep*mcf(atradoff(ring))*frf;
rgup=atsetcavity(ring,'Frequency',frf+0.5*DFStep,cavargs{:});
rgdn=atsetcavity(ring,'Frequency',frf-0.5*DFStep,cavargs{:});
[~,o1P]=findorbit(rgup,'guess',orbitin,varargs{:});
[~,o1M]=findorbit(rgdn,'guess',orbitin,varargs{:});
[~,o1P]=findorbit6(rgup,'guess',orbitin,varargs{:});
[~,o1M]=findorbit6(rgdn,'guess',orbitin,varargs{:});
deltap=o1P(5)-o1M(5);
else
dp=orbitin(5);
Expand Down
Loading