FINDM44 numerically finds the 4x4 transfer matrix of an accelerator lattice
for a particle with relative momentum deviation DP
IMPORTANT!!! FINDM44 assumes constant momentum deviation.
PassMethod used for any element in the LATTICE SHOULD NOT
1.change the longitudinal momentum dP
(cavities , magnets with radiation, ...)
2.have any time dependence (localized impedance, fast kickers, ...)
M44 = FINDM44(LATTICE,DP) finds a full one-turn
matrix at the entrance of the first element
!!! With this syntax FINDM44 assumes that the LATTICE
is a ring and first finds the closed orbit
[M44,T] = FINDM44(LATTICE,DP,REFPTS) also returns
4-by-4 transfer matrixes between entrance of
the first element and each element indexed by REFPTS.
T is 4-by-4-by-length(REFPTS) 3 dimensional array
so that the set of indexes (:,:,i) selects the 4-by-4
matrix at the i-th reference point.
Note: REFPTS is an array of increasing indexes that
select elements from range 1 to length(LATTICE)+1.
See further explanation of REFPTS in the 'help' for FINDSPOS
When REFPTS= [ 1 2 .. ] the fist point is the entrance of the
first element and T(:,:,1) - identity matrix
Note: REFPTS is allowed to go 1 point beyond the
number of elements. In this case the last point is
the EXIT of the last element. If LATTICE is a RING
it is also the entrance of the first element
after 1 turn: T(:,:,end) = M
[M44,T] = FINDM44(LATTICE,DP,REFPTS,ORBITIN) - Does not search for
closed orbit. Instead the ORBITIN,a 1-by-6 vector of initial
conditions is used: [x0, px0, y0, py0, DP, 0]' where
the same DP as argument 2. The sixth component is ignored.
This syntax is useful to specify the entrance orbit
if LATTICE is not a ring or to avoid recomputting the
closed orbit if is already known.
[M44,T] = FINDM44(LATTICE,DP,REFPTS,ORBITIN,'full') - same as above except
matrixes returned in T are full 1-turn matrixes at the entrance of each
element indexed by REFPTS.
[M44,T,orbit] = FINDM44(...) in addition returns
at REFPTS the closed orbit calculated along the
way with findorbit4
See also FINDM66, FINDORBIT40001 function [M44, varargout] = findm44(LATTICE,DP,varargin) 0002 %FINDM44 numerically finds the 4x4 transfer matrix of an accelerator lattice 0003 % for a particle with relative momentum deviation DP 0004 % 0005 % IMPORTANT!!! FINDM44 assumes constant momentum deviation. 0006 % PassMethod used for any element in the LATTICE SHOULD NOT 0007 % 1.change the longitudinal momentum dP 0008 % (cavities , magnets with radiation, ...) 0009 % 2.have any time dependence (localized impedance, fast kickers, ...) 0010 % 0011 % M44 = FINDM44(LATTICE,DP) finds a full one-turn 0012 % matrix at the entrance of the first element 0013 % !!! With this syntax FINDM44 assumes that the LATTICE 0014 % is a ring and first finds the closed orbit 0015 % 0016 % [M44,T] = FINDM44(LATTICE,DP,REFPTS) also returns 0017 % 4-by-4 transfer matrixes between entrance of 0018 % the first element and each element indexed by REFPTS. 0019 % T is 4-by-4-by-length(REFPTS) 3 dimensional array 0020 % so that the set of indexes (:,:,i) selects the 4-by-4 0021 % matrix at the i-th reference point. 0022 % 0023 % Note: REFPTS is an array of increasing indexes that 0024 % select elements from range 1 to length(LATTICE)+1. 0025 % See further explanation of REFPTS in the 'help' for FINDSPOS 0026 % When REFPTS= [ 1 2 .. ] the fist point is the entrance of the 0027 % first element and T(:,:,1) - identity matrix 0028 % 0029 % Note: REFPTS is allowed to go 1 point beyond the 0030 % number of elements. In this case the last point is 0031 % the EXIT of the last element. If LATTICE is a RING 0032 % it is also the entrance of the first element 0033 % after 1 turn: T(:,:,end) = M 0034 % 0035 % [M44,T] = FINDM44(LATTICE,DP,REFPTS,ORBITIN) - Does not search for 0036 % closed orbit. Instead the ORBITIN,a 1-by-6 vector of initial 0037 % conditions is used: [x0, px0, y0, py0, DP, 0]' where 0038 % the same DP as argument 2. The sixth component is ignored. 0039 % This syntax is useful to specify the entrance orbit 0040 % if LATTICE is not a ring or to avoid recomputting the 0041 % closed orbit if is already known. 0042 % 0043 % [M44,T] = FINDM44(LATTICE,DP,REFPTS,ORBITIN,'full') - same as above except 0044 % matrixes returned in T are full 1-turn matrixes at the entrance of each 0045 % element indexed by REFPTS. 0046 % 0047 % [M44,T,orbit] = FINDM44(...) in addition returns 0048 % at REFPTS the closed orbit calculated along the 0049 % way with findorbit4 0050 % 0051 % See also FINDM66, FINDORBIT4 0052 0053 % ************************************************************************* 0054 % The numerical differentiation in FINDM44 uses symmetric form 0055 % 0056 % F(x+delta) - F(x-delta) 0057 % -------------------------------------- 0058 % 2*delta 0059 % 0060 % with optimal differentiation step delta given by !!!! DO LATER 0061 % The relative error in the derivative computed this way 0062 % is !!!!!!!!!!!!!!!!! DO LATER 0063 % Reference: Numerical Recipes. 0064 0065 0066 if ~iscell(LATTICE) 0067 error('First argument must be a cell array'); 0068 end 0069 0070 NE = length(LATTICE); 0071 0072 if nargin >= 5 % FINDM44(LATTICE,DP,REFPTS,ORBITIN,'full') 0073 if strcmpi(varargin{3},'full') 0074 FULLFLAG=1; 0075 else 0076 error('Fifth argument - unknown option'); 0077 end 0078 else 0079 FULLFLAG=0; 0080 end 0081 if nargin >= 4 && ~isempty(varargin{2}) % FINDM44(LATTICE,DP,REFPTS,ORBITIN) 0082 R0 = varargin{2}; 0083 R0(5) = DP; 0084 R0(6) = 0; 0085 else 0086 R0 = [findorbit4(LATTICE,DP);DP;0]; 0087 end 0088 if nargin >= 3 % FINDM44(LATTICE,DP,REFPTS) 0089 if islogical(varargin{1}) 0090 REFPTS=varargin{1}; 0091 REFPTS(end+1:NE+1)=false; 0092 elseif isnumeric(varargin{1}) 0093 REFPTS=setelems(false(1,NE+1),varargin{1}); 0094 else 0095 error('REFPTS must be numeric or logical'); 0096 end 0097 else 0098 REFPTS=false(1,NE+1); 0099 end 0100 refs=setelems(REFPTS,NE+1); 0101 reqs=REFPTS(refs); 0102 0103 0104 % Determine step size to use for numerical differentiation 0105 global NUMDIFPARAMS 0106 % Transverse 0107 if isfield(NUMDIFPARAMS,'XYStep') 0108 dt = NUMDIFPARAMS.XYStep'; 0109 else 0110 % optimal differentiation step - Numerical Recipes 0111 dt = 6.055454452393343e-006; 0112 end 0113 0114 % Build a diagonal matrix of initial conditions 0115 D4 = [0.5*dt*eye(4);zeros(2,4)]; 0116 % Add to the orbit_in. First 8 columns for derivative 0117 % 9-th column is for closed orbit 0118 RIN = R0(:,ones(1,9)) + [D4 -D4 zeros(6,1)]; 0119 ROUT = linepass(LATTICE,RIN,refs); 0120 TMAT3 = reshape(ROUT(1:4,:),4,9,[]); 0121 M44 = (TMAT3(:,1:4,end)-TMAT3(:,5:8,end))./dt; 0122 0123 if nargout >= 2 % Calculate matrixes at all REFPTS. 0124 MSTACK = (TMAT3(:,1:4,reqs)-TMAT3(:,5:8,reqs))./dt; 0125 0126 if FULLFLAG 0127 S2 = [0 1;-1 0]; 0128 S4 = [S2, zeros(2);zeros(2),S2]; % symplectic identity matrix 0129 v=cellfun(@rotate,num2cell(MSTACK,[1 2]),'UniformOutput',false); 0130 varargout{1}=cat(3,v{:}); 0131 else 0132 varargout{1}=MSTACK; 0133 end 0134 % return the closed orbit if requested 0135 if nargout == 3 0136 varargout{2}=squeeze(TMAT3(:,9,reqs)); 0137 end 0138 0139 end 0140 0141 function mask=setelems(mask,idx) 0142 mask(idx)=true; 0143 end 0144 0145 function w=rotate(v) 0146 w=v*M44*S4'*v'*S4; 0147 end 0148 end