ATLINOPT performs linear analysis of the COUPLED lattices
LinData = ATLINOPT(RING,DP,REFPTS) is a MATLAB structure array with fields
ElemIndex - ordinal position in the RING
SPos - longitudinal position [m]
ClosedOrbit - closed orbit column vector with
components x, px, y, py (momentums, NOT angles)
Dispersion - dispersion orbit position vector with
components eta_x, eta_prime_x, eta_y, eta_prime_y
calculated with respect to the closed orbit with
momentum deviation DP. Only if chromaticity is required.
M44 - 4x4 transfer matrix M from the beginning of RING
to the entrance of the element for specified DP [2]
A - 2x2 matrix A in [3]
B - 2x2 matrix B in [3]
C - 2x2 matrix C in [3]
gamma - gamma parameter of the transformation to eigenmodes
mu - [ mux, muy] horizontal and vertical betatron phase
beta - [betax, betay] vector
alpha - [alphax, alphay] vector
All values are specified at the entrance of each element specified in REFPTS.
REFPTS is an array of increasing indexes that select elements
from the range 1 to length(LINE)+1. Defaults to 1 (initial point)
See further explanation of REFPTS in the 'help' for FINDSPOS
[LinData,NU] = ATLINOPT() returns a vector of linear tunes
[nu_u , nu_v] for two normal modes of linear motion [1]
[LinData,NU, KSI] = ATLINOPT() returns a vector of chromaticities ksi = d(nu)/(dP/P)
[ksi_u , ksi_v] - derivatives of [nu_u , nu_v]
LinData = ATLINOPT(RING,DP,REFPTS,ORBITIN) does not search for closed orbit.
instead ORBITIN is used
Difference with linopt: Fractional tunes 0<=tune<1
Dispersion output (if chromaticity is required)
Alpha output
Phase advance output
Option to skip closed orbit search
See also ATX ATMODUL FINDSPOS TWISSRING TUNECHROM
[1] D.Edwars,L.Teng IEEE Trans.Nucl.Sci. NS-20, No.3, p.885-888, 1973
[2] E.Courant, H.Snyder
[3] D.Sagan, D.Rubin Phys.Rev.Spec.Top.-Accelerators and beams, vol.2 (1999)0001 function [lindata, varargout] = atlinopt(RING,DP,varargin) 0002 %ATLINOPT performs linear analysis of the COUPLED lattices 0003 % 0004 % LinData = ATLINOPT(RING,DP,REFPTS) is a MATLAB structure array with fields 0005 % 0006 % ElemIndex - ordinal position in the RING 0007 % SPos - longitudinal position [m] 0008 % ClosedOrbit - closed orbit column vector with 0009 % components x, px, y, py (momentums, NOT angles) 0010 % Dispersion - dispersion orbit position vector with 0011 % components eta_x, eta_prime_x, eta_y, eta_prime_y 0012 % calculated with respect to the closed orbit with 0013 % momentum deviation DP. Only if chromaticity is required. 0014 % M44 - 4x4 transfer matrix M from the beginning of RING 0015 % to the entrance of the element for specified DP [2] 0016 % A - 2x2 matrix A in [3] 0017 % B - 2x2 matrix B in [3] 0018 % C - 2x2 matrix C in [3] 0019 % gamma - gamma parameter of the transformation to eigenmodes 0020 % mu - [ mux, muy] horizontal and vertical betatron phase 0021 % beta - [betax, betay] vector 0022 % alpha - [alphax, alphay] vector 0023 % 0024 % All values are specified at the entrance of each element specified in REFPTS. 0025 % REFPTS is an array of increasing indexes that select elements 0026 % from the range 1 to length(LINE)+1. Defaults to 1 (initial point) 0027 % See further explanation of REFPTS in the 'help' for FINDSPOS 0028 % 0029 % [LinData,NU] = ATLINOPT() returns a vector of linear tunes 0030 % [nu_u , nu_v] for two normal modes of linear motion [1] 0031 % 0032 % [LinData,NU, KSI] = ATLINOPT() returns a vector of chromaticities ksi = d(nu)/(dP/P) 0033 % [ksi_u , ksi_v] - derivatives of [nu_u , nu_v] 0034 % 0035 % LinData = ATLINOPT(RING,DP,REFPTS,ORBITIN) does not search for closed orbit. 0036 % instead ORBITIN is used 0037 % 0038 % Difference with linopt: Fractional tunes 0<=tune<1 0039 % Dispersion output (if chromaticity is required) 0040 % Alpha output 0041 % Phase advance output 0042 % Option to skip closed orbit search 0043 % 0044 % See also ATX ATMODUL FINDSPOS TWISSRING TUNECHROM 0045 % 0046 % [1] D.Edwars,L.Teng IEEE Trans.Nucl.Sci. NS-20, No.3, p.885-888, 1973 0047 % [2] E.Courant, H.Snyder 0048 % [3] D.Sagan, D.Rubin Phys.Rev.Spec.Top.-Accelerators and beams, vol.2 (1999) 0049 0050 0051 global NUMDIFPARAMS 0052 0053 NE = length(RING); 0054 if nargin >= 3 0055 if islogical(varargin{1}) 0056 REFPTS=varargin{1}(:); 0057 REFPTS(end+1:NE+1)=false; 0058 isel=1:sum(REFPTS); 0059 elseif isnumeric(varargin{1}) 0060 isel=varargin{1}; 0061 [r2,ir2]=sort(isel(:)); 0062 REFPTS=setelems(false(NE+1,1),r2); 0063 isel(ir2)=1:numel(isel); 0064 else 0065 error('REFPTS must be numeric or logical'); 0066 end 0067 else 0068 REFPTS=setelems(false(NE+1,1),1); 0069 isel=1; 0070 end 0071 0072 spos = findspos(RING,REFPTS); 0073 [M44, MS, orb] = findm44(RING,DP,REFPTS,varargin{2:end}); 0074 0075 % Calculate A,B,C, gamma at the first element 0076 M =M44(1:2,1:2); 0077 N =M44(3:4,3:4); 0078 m =M44(1:2,3:4); 0079 n =M44(3:4,1:2); 0080 0081 % 2-by-2 symplectic matrix 0082 S = [0 1; -1 0]; 0083 H = m + S*n'*S'; 0084 t = trace(M-N); 0085 0086 g = sqrt(1 + sqrt(t*t/(t*t+4*det(H))))/sqrt(2); 0087 G = diag([g g]); 0088 C = -H*sign(t)/(g*sqrt(t*t+4*det(H))); 0089 A = G*G*M - G*(m*S*C'*S' + C*n) + C*N*S*C'*S'; 0090 B = G*G*N + G*(S*C'*S'*m + n*C) + S*C'*S'*M*C; 0091 [MSA,MSB,gamma,CL,AL,BL]=cellfun(@analyze,num2cell(MS,[1 2]),'UniformOutput',false); 0092 0093 0094 [BX,AX,MX]=lop(reshape(cat(3,MSA{:}),2,2,[]),A); 0095 [BY,AY,MY]=lop(reshape(cat(3,MSB{:}),2,2,[]),B); 0096 0097 %tunes = [MX(end),MY(end)]/2/pi; 0098 0099 if nargout >= 2 0100 cos_mu_x = trace(A)/2; 0101 cos_mu_y = trace(B)/2; 0102 tns = acos([cos_mu_x cos_mu_y])/2/pi; 0103 if A(1,2) < 0, tns(1)=1-tns(1); end 0104 if B(1,2) < 0, tns(2)=1-tns(2); end 0105 varargout{1}=tns; 0106 end 0107 0108 dispargs={}; 0109 if nargout >= 3 0110 if isfield(NUMDIFPARAMS,'DPStep') 0111 dDP = NUMDIFPARAMS.DPStep'; 0112 else 0113 dDP = 1e-6; 0114 end 0115 refs=false(1,length(RING)+1); 0116 % Calculate tunes for DP+dDP 0117 [orbP,o1P]=findorbit4(RING,DP+0.5*dDP,REFPTS); 0118 [orbM,o1M]=findorbit4(RING,DP-0.5*dDP,REFPTS); 0119 dispersion = (orbP-orbM)/dDP; 0120 [LD, tunesP] = atlinopt(RING,DP+0.5*dDP,refs,o1P); %#ok<ASGLU> 0121 [LD, tunesM] = atlinopt(RING,DP-0.5*dDP,refs,o1M); %#ok<ASGLU> 0122 varargout{2} = (tunesP - tunesM)/dDP; 0123 dispargs={'Dispersion',num2cell(dispersion,1)'}; 0124 end 0125 ld = struct('ElemIndex',num2cell(find(REFPTS)),... 0126 'SPos',num2cell(spos)',... 0127 'ClosedOrbit',num2cell(orb,1)',... 0128 dispargs{:},... 0129 'M44',squeeze(num2cell(MS,[1 2])),... 0130 'gamma',gamma(:),... 0131 'C',squeeze(CL),... 0132 'A',squeeze(AL),... 0133 'B',squeeze(BL),... 0134 'beta',num2cell([BX,BY],2),... 0135 'alpha',num2cell([AX,AY],2),... 0136 'mu',num2cell([MX,MY],2)); 0137 0138 lindata=reshape(ld(isel),size(isel)); 0139 0140 function [MSA,MSB,gamma,CL,AL,BL]=analyze(MS) 0141 M12 =MS(1:2,1:2); 0142 N12 =MS(3:4,3:4); 0143 m12 =MS(1:2,3:4); 0144 n12 =MS(3:4,1:2); 0145 0146 gamma = sqrt(det(n12*C+G*N12)); 0147 MSA = (G*M12-m12*S*C'*S')/gamma; 0148 MSB = (n12*C+G*N12)/gamma; 0149 0150 CL=(M12*C+G*m12)*S*MSB'*S'; 0151 AL=MSA*A*S*MSA'*S'; 0152 BL=MSB*B*S*MSB'*S'; 0153 end 0154 0155 function mask=setelems(mask,idx) 0156 mask(idx)=true; 0157 end 0158 0159 function UP = BetatronPhaseUnwrap(P) 0160 % unwrap negative jumps in betatron 0161 %JUMPS = [0; diff(P)] < -1.e-5; 0162 JUMPS = [0; diff(P)] < -1.e-3; 0163 UP = P+cumsum(JUMPS)*2*pi; 0164 end 0165 0166 function [beta,alpha,phase]=lop(MS,A0) 0167 sinmu = sign(A0(1,2))*sqrt(-A0(1,2)*A0(2,1)-(A0(1,1)-A0(2,2))^2/4); 0168 ax = (A0(1,1)-A0(2,2))/2/sinmu; 0169 bx = A0(1,2)/sinmu; 0170 bbb=squeeze(MS(1,2,:)); 0171 aaa=squeeze(MS(1,1,:))*bx-bbb*ax; 0172 0173 beta = (aaa.^2 + bbb.^2)/bx; 0174 alpha = -(aaa.*squeeze(MS(2,1,:)*bx-MS(2,2,:)*ax) + bbb.*squeeze(MS(2,2,:)))/bx; 0175 try 0176 phase = atan2(bbb,aaa); 0177 catch %#ok<CTCH> 0178 phase=NaN(size(beta)); 0179 end 0180 phase = BetatronPhaseUnwrap(phase); 0181 end 0182 0183 end