1 | function [M44, varargout] = findm44(LATTICE,DP,varargin) |
---|
2 | %FINDM44 numerically finds the 4x4 transfer matrix of an accelerator lattice |
---|
3 | % for a particle with relative momentum deviation DP |
---|
4 | % |
---|
5 | % IMPORTANT!!! FINDM44 assumes constant momentum deviation. |
---|
6 | % PassMethod used for any element in the LATTICE SHOULD NOT |
---|
7 | % 1.change the longitudinal momentum dP |
---|
8 | % (cavities , magnets with radiation, ...) |
---|
9 | % 2.have any time dependence (localized impedance, fast kickers, ...) |
---|
10 | % |
---|
11 | % M44 = FINDM44(LATTICE,DP) finds a full one-turn |
---|
12 | % matrix at the entrance of the first element |
---|
13 | % !!! With this syntax FINDM44 assumes that the LATTICE |
---|
14 | % is a ring and first finds the closed orbit |
---|
15 | % |
---|
16 | % [M44,T] = FINDM44(LATTICE,DP,REFPTS) also returns |
---|
17 | % 4-by-4 transfer matrixes between entrance of |
---|
18 | % the first element and each element indexed by REFPTS. |
---|
19 | % T is 4-by-4-by-length(REFPTS) 3 dimensional array |
---|
20 | % so that the set of indexes (:,:,i) selects the 4-by-4 |
---|
21 | % matrix at the i-th reference point. |
---|
22 | % |
---|
23 | % Note: REFPTS is an array of increasing indexes that |
---|
24 | % select elements from range 1 to length(LATTICE)+1. |
---|
25 | % See further explanation of REFPTS in the 'help' for FINDSPOS |
---|
26 | % When REFPTS= [ 1 2 .. ] the first point is the entrance of the |
---|
27 | % first element and T(:,:,1) - identity matrix |
---|
28 | % |
---|
29 | % Note: REFPTS is allowed to go 1 point beyond the |
---|
30 | % number of elements. In this case the last point is |
---|
31 | % the EXIT of the last element. If LATTICE is a RING |
---|
32 | % it is also the entrance of the first element |
---|
33 | % after 1 turn: T(:,:,end) = M |
---|
34 | % |
---|
35 | % [M44, T] = FINDM44(LATTICE,DP,REFPTS,ORBITIN) - Does not search for |
---|
36 | % closed orbit. Instead the ORBITIN,a 1-by-6 vector of initial |
---|
37 | % conditions is used: [x0, px0, y0, py0, DP, 0]' where |
---|
38 | % the same DP as argument 2. The sixth component is ignored. |
---|
39 | % This syntax is useful to specify the entrance orbit |
---|
40 | % if LATTICE is not a ring or to avoid recomputing the |
---|
41 | % closed orbit if is already known. |
---|
42 | % |
---|
43 | % [M44, T] = FINDM44(LATTICE,DP,REFPTS,ORBITIN,'full') - same as above except |
---|
44 | % matrixes returned in T are full 1-turn matrixes at the entrance of each |
---|
45 | % element indexed by REFPTS. |
---|
46 | % |
---|
47 | % [M44, T, orbit] = FINDM44(...) in addition returns |
---|
48 | % at REFPTS the closed orbit calculated along the |
---|
49 | % way with findorbit4 |
---|
50 | % |
---|
51 | % See also LINEPASS, LATTICEPASS, FINDORBIT4, FINDSPOS |
---|
52 | |
---|
53 | % ************************************************************************* |
---|
54 | % The numerical differentiation in FINDM44 uses symmetric form |
---|
55 | % |
---|
56 | % F(x+delta) - F(x-delta) |
---|
57 | % -------------------------------------- |
---|
58 | % 2*delta |
---|
59 | % |
---|
60 | % with optimal differentiation step delta given by !!!! DO LATER |
---|
61 | % The relative error in the derivative computed this way |
---|
62 | % is !!!!!!!!!!!!!!!!! DO LATER |
---|
63 | % Reference: Numerical Recipes. |
---|
64 | |
---|
65 | |
---|
66 | if ~iscell(LATTICE) |
---|
67 | error('First argument must be a cell array'); |
---|
68 | end |
---|
69 | |
---|
70 | NE = length(LATTICE); |
---|
71 | |
---|
72 | |
---|
73 | switch nargin |
---|
74 | case 5 % FINDM44(LATTICE,DP,REFPTS,ORBITIN,'full') |
---|
75 | if(lower(varargin{3})=='full') |
---|
76 | FULLFLAG = 1; |
---|
77 | REFPTS = varargin{1}; |
---|
78 | R0 = varargin{2}; |
---|
79 | R0(5) = DP; |
---|
80 | R0(6)= 0; |
---|
81 | else |
---|
82 | error('Fifth argument - unknown option') |
---|
83 | end |
---|
84 | case 4 % FINDM44(LATTICE,DP,REFPTS,ORBITIN) |
---|
85 | FULLFLAG = 0; |
---|
86 | REFPTS = varargin{1}; |
---|
87 | R0 = varargin{2}; |
---|
88 | R0(5) = DP; |
---|
89 | R0(6)= 0; |
---|
90 | case 3 % FINDM44(LATTICE,DP,REFPTS) |
---|
91 | FULLFLAG = 0; |
---|
92 | REFPTS = varargin{1}; |
---|
93 | R0 = [findorbit4(LATTICE,DP);DP;0]; |
---|
94 | case 2 % FINDM44(LATTICE,DP) |
---|
95 | REFPTS = NE+1; |
---|
96 | FULLFLAG = 0; |
---|
97 | R0 = [findorbit4(LATTICE,DP);DP;0]; |
---|
98 | otherwise |
---|
99 | error('Incorrect number of input arguments'); |
---|
100 | end |
---|
101 | |
---|
102 | NR = length(REFPTS); |
---|
103 | |
---|
104 | |
---|
105 | % Determine step size to use for numerical differentiation |
---|
106 | global NUMDIFPARAMS |
---|
107 | |
---|
108 | if isfield(NUMDIFPARAMS,'XYStep') |
---|
109 | d = NUMDIFPARAMS.XYStep'; |
---|
110 | else |
---|
111 | % optimal differentiation step - Numerical Recipes |
---|
112 | %d = 6.055454452393343e-006; |
---|
113 | |
---|
114 | % the same value as in dt/dl in findm66.m |
---|
115 | % to test why the chrom. from findm44 & findm66 are different |
---|
116 | % need to convert back to the numerical recipes value |
---|
117 | % when the test is finished... |
---|
118 | d = 1e-8; |
---|
119 | end |
---|
120 | |
---|
121 | |
---|
122 | % Put together matrix of initial conditions |
---|
123 | |
---|
124 | D = d*eye(4); |
---|
125 | % First 8 columns for derivative |
---|
126 | % 9-th column is for closed orbit |
---|
127 | % R0 is the closed orbit |
---|
128 | RM = [[R0 R0 R0 R0 R0 R0 R0 R0] + [D -D; zeros(2,8)],R0]; |
---|
129 | |
---|
130 | if nargout < 2 |
---|
131 | % Calculate M44 at the first element only. Use linepass |
---|
132 | TMAT = linepass(LATTICE,RM); |
---|
133 | M44 = (TMAT(1:4,1:4)-TMAT(1:4,5:8))/(2*d); |
---|
134 | return |
---|
135 | else |
---|
136 | % Calculate matrices at all REFPTS. Use linepass |
---|
137 | % Need to include the exit of the LATTICE to REFPTS array |
---|
138 | if(REFPTS(NR)~=NE+1) |
---|
139 | REFPTS = [REFPTS NE+1]; |
---|
140 | NR1 = NR+1; |
---|
141 | else |
---|
142 | NR1 = NR; |
---|
143 | end |
---|
144 | |
---|
145 | TMAT = linepass(LATTICE,RM,REFPTS); |
---|
146 | TMAT3 = reshape(TMAT(1:4,:),4,9,NR1); |
---|
147 | M44 = (TMAT3(1:4,1:4,NR1)-TMAT3(1:4,5:8,NR1))/(2*d); |
---|
148 | |
---|
149 | MSTACK = (TMAT3(:,1:4,1:NR)-TMAT3(:,5:8,1:NR))/(2*d); |
---|
150 | |
---|
151 | if FULLFLAG |
---|
152 | S2 = [0 1;-1 0]; |
---|
153 | S4 = [S2, zeros(2);zeros(2),S2]; % symplectic identity matrix |
---|
154 | for k =1:NR |
---|
155 | T = MSTACK(:,:,k); |
---|
156 | varargout{1}(:,:,k) = T*M44*S4'*T'*S4; |
---|
157 | end |
---|
158 | else |
---|
159 | varargout{1}=MSTACK; |
---|
160 | end |
---|
161 | % return the closed orbit if requested |
---|
162 | if nargout == 3 |
---|
163 | varargout{2}=squeeze(TMAT3(:,9,1:NR)); |
---|
164 | end |
---|
165 | end |
---|