TWISSLINE calculates linear optics functions for an UNCOUPLED transport line
TwissData = TWISSLINE(LATTICE,DP,TWISSDATAIN) propagates twiss
parameters and closed orbit coordinates from the LINE entrance
given by TWISSDATAIN assuming constant energy deviation DP.
TWISSDATAIN is a 1-by-1 structure with the same field names
as the return argument. (See below)
!!! IMPORTANT: Since TWISSLINE does not search for closed orbit
its value at the entrance must be supplied in the
ClosedOrbit field of TWISSDATAIN structure.
TwissData = TWISSLINE(LATTICE,DP,TWISSDATAIN,REFPTS) calculates Twiss parameters
and closed orbit coordinates at specified reference points REFPTS
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
TwissData = TWISSLINE(...,'chrom', DDP) also calculates
linear dispersion. Dispersion is returned as one
of the fields in TwissData.
!!! Last argument DDP is a momentum deviation on top
of DP (the second argument) used to calculate and normalize
dispersion. If not supplied
the default value of 1e-8 is used.
TwisData is a 1-by-REFPTS (1-by-1 if no REFPTS specified) structure array with fields:
ElemIndex - integer (element number) in the LINE
SPos - longitudinal position [m]
ClosedOrbit - closed orbit column vector with
components x, px, y, py (momentums, NOT angles)
Dispersion - dispersion orbit position 4-by-1 vector with
components [eta_x, eta_prime_x, eta_y, eta_prime_y]'
calculated with respect to the closed orbit with
momentum deviation DP
M44 - 4x4 transfer matrix M from the beginning of LINE
to the entrance of the element for specified DP [2]
beta - [betax, betay] horizontal and vertical Twiss parameter beta
alpha - [alphax, alphay] horizontal and vertical Twiss parameter alpha
mu - [mux, muy] horizontal and vertical betatron phase
!!! NOT 2*PI normalized
Use CAT to get the data from fields of TwissData into MATLAB arrays.
Example:
>> TD = twissring(THERING,0,1:length(THERING));
>> BETA = cat(1,TD.beta);
>> S = cat(1,TD.SPos);
>> plot(S,BETA(:,1))
See also TWISSRING, LINOPT, TUNECHROM.0001 function [TD, varargout] = twissline(LINE,DP,TWISSDATAIN,varargin) 0002 %TWISSLINE calculates linear optics functions for an UNCOUPLED transport line 0003 % 0004 % TwissData = TWISSLINE(LATTICE,DP,TWISSDATAIN) propagates twiss 0005 % parameters and closed orbit coordinates from the LINE entrance 0006 % given by TWISSDATAIN assuming constant energy deviation DP. 0007 % TWISSDATAIN is a 1-by-1 structure with the same field names 0008 % as the return argument. (See below) 0009 % !!! IMPORTANT: Since TWISSLINE does not search for closed orbit 0010 % its value at the entrance must be supplied in the 0011 % ClosedOrbit field of TWISSDATAIN structure. 0012 % 0013 % TwissData = TWISSLINE(LATTICE,DP,TWISSDATAIN,REFPTS) calculates Twiss parameters 0014 % and closed orbit coordinates at specified reference points REFPTS 0015 % 0016 % Note: REFPTS is an array of increasing indexes that 0017 % select elements from range 1 to length(LATTICE)+1. 0018 % See further explanation of REFPTS in the 'help' for FINDSPOS 0019 % 0020 % TwissData = TWISSLINE(...,'chrom', DDP) also calculates 0021 % linear dispersion. Dispersion is returned as one 0022 % of the fields in TwissData. 0023 % !!! Last argument DDP is a momentum deviation on top 0024 % of DP (the second argument) used to calculate and normalize 0025 % dispersion. If not supplied 0026 % the default value of 1e-8 is used. 0027 % 0028 % TwisData is a 1-by-REFPTS (1-by-1 if no REFPTS specified) structure array with fields: 0029 % ElemIndex - integer (element number) in the LINE 0030 % SPos - longitudinal position [m] 0031 % ClosedOrbit - closed orbit column vector with 0032 % components x, px, y, py (momentums, NOT angles) 0033 % Dispersion - dispersion orbit position 4-by-1 vector with 0034 % components [eta_x, eta_prime_x, eta_y, eta_prime_y]' 0035 % calculated with respect to the closed orbit with 0036 % momentum deviation DP 0037 % M44 - 4x4 transfer matrix M from the beginning of LINE 0038 % to the entrance of the element for specified DP [2] 0039 % beta - [betax, betay] horizontal and vertical Twiss parameter beta 0040 % alpha - [alphax, alphay] horizontal and vertical Twiss parameter alpha 0041 % mu - [mux, muy] horizontal and vertical betatron phase 0042 % !!! NOT 2*PI normalized 0043 % 0044 % Use CAT to get the data from fields of TwissData into MATLAB arrays. 0045 % Example: 0046 % >> TD = twissring(THERING,0,1:length(THERING)); 0047 % >> BETA = cat(1,TD.beta); 0048 % >> S = cat(1,TD.SPos); 0049 % >> plot(S,BETA(:,1)) 0050 % 0051 % See also TWISSRING, LINOPT, TUNECHROM. 0052 0053 0054 DDP_default = 1e-8; 0055 NE=length(LINE); 0056 0057 % Process input arguments 0058 switch nargin 0059 case 3 0060 REFPTS=NE+1; 0061 CHROMFLAG=0; 0062 case 4 0063 if isnumeric(varargin{1}) 0064 REFPTS = varargin{1}; 0065 CHROMFLAG = 0; 0066 elseif ischar(varargin{1}) & strncmp(lower(varargin{1}),'chrom',5) 0067 CHROMFLAG = 1; 0068 REFPTS = NE+1; 0069 DDP = DDP_default; 0070 else 0071 error('Third argument must be a numeric array or string'); 0072 end 0073 case 5 0074 if isnumeric(varargin{1}) 0075 REFPTS = varargin{1}; 0076 if ischar(varargin{2}) & strncmp(lower(varargin{2}),'chrom',5) 0077 CHROMFLAG = 1; 0078 DDP = DDP_default; 0079 else 0080 error('Fourth argument - wrong type'); 0081 end 0082 elseif ischar(varargin{1}) & strncmp(lower(varargin{1}),'chrom',5) 0083 CHROMFLAG = 1; 0084 REFPTS = NE+1; 0085 if isnumeric(varargin{2}) 0086 DDP = varargin{2}; 0087 else 0088 error('Fourth argument - wrong type'); 0089 end 0090 end 0091 case 6 0092 if isnumeric(varargin{1}) 0093 REFPTS = varargin{1}; 0094 else 0095 error('Fourth argument - wrong type'); 0096 end 0097 if ischar(varargin{2}) & strncmp(lower(varargin{2}),'chrom',5) 0098 CHROMFLAG = 1; 0099 else 0100 error('Fifth argument - wrong type'); 0101 end 0102 if isnumeric(varargin{3}) 0103 DDP = varargin{3}; 0104 else 0105 error('Sixth argument - wrong type'); 0106 end 0107 otherwise 0108 error('Wrong number of arguments'); 0109 end 0110 0111 0112 0113 0114 if isfield(TWISSDATAIN,'alpha') 0115 ax = TWISSDATAIN(end).alpha(1); 0116 ay = TWISSDATAIN(end).alpha(2); 0117 else 0118 error('TWISSDATAIN structure does not have field ''alpha'''); 0119 end 0120 0121 if isfield(TWISSDATAIN,'beta') 0122 bx = TWISSDATAIN(end).beta(1); 0123 by = TWISSDATAIN(end).beta(2); 0124 else 0125 error('TWISSDATAIN structure does not have field ''beta'''); 0126 end 0127 0128 if isfield(TWISSDATAIN,'mu') 0129 mux = TWISSDATAIN(end).mu(1); 0130 muy = TWISSDATAIN(end).mu(2); 0131 else 0132 error('TWISSDATAIN structure does not have field ''mu'''); 0133 end 0134 0135 R0 = [TWISSDATAIN(end).ClosedOrbit;DP;0]; 0136 0137 [M44, MS, orb] = findm44(LINE,DP,REFPTS,R0); 0138 0139 BX = squeeze((MS(1,1,:)*bx-MS(1,2,:)*ax).^2 + MS(1,2,:).^2)/bx; 0140 BY = squeeze((MS(3,3,:)*by-MS(3,4,:)*ay).^2 + MS(3,4,:).^2)/by; 0141 0142 AX = -squeeze((MS(1,1,:)*bx-MS(1,2,:)*ax).*(MS(2,1,:)*bx-MS(2,2,:)*ax) + MS(1,2,:).*MS(2,2,:))/bx; 0143 AY = -squeeze((MS(3,3,:)*by-MS(3,4,:)*ay).*(MS(4,3,:)*by-MS(4,4,:)*ay) + MS(3,4,:).*MS(4,4,:))/by; 0144 0145 MX = atan(squeeze( MS(1,2,:)./(MS(1,1,:)*bx-MS(1,2,:)*ax))); 0146 MY = atan(squeeze(MS(3,4,:)./(MS(3,3,:)*by-MS(3,4,:)*ay))); 0147 0148 MX = BetatronPhaseUnwrap(MX); 0149 MY = BetatronPhaseUnwrap(MY); 0150 0151 TD = struct('ElemIndex',num2cell(REFPTS),... 0152 'SPos',num2cell(findspos(LINE,REFPTS)),... 0153 'ClosedOrbit',num2cell(orb,1),... 0154 'M44', squeeze(num2cell(MS,[1 2]))',... 0155 'beta', num2cell([BX,BY],2)',... 0156 'alpha', num2cell([AX,AY],2)',... 0157 'mu', num2cell([MX,MY],2)'); 0158 0159 0160 if CHROMFLAG 0161 TWISSDATAIN_DDP = TWISSDATAIN(end); 0162 TWISSDATAIN_DDP.ClosedOrbit = TWISSDATAIN_DDP.ClosedOrbit+TWISSDATAIN_DDP.Dispersion(:)*DDP; 0163 TD_DDP = twissline(LINE,DP+DDP,TWISSDATAIN_DDP,REFPTS); 0164 DORBIT = reshape(cat(1,TD_DDP.ClosedOrbit),4,length(cat(1,TD_DDP.ClosedOrbit))/4); 0165 DISPERSION = num2cell((DORBIT-orb)/DDP,1); 0166 [TD.Dispersion] = deal( DISPERSION{:}); 0167 0168 end 0169 0170 0171 0172 function UP = BetatronPhaseUnwrap(P) 0173 % unwrap negative jumps in betatron phase 0174 DP = diff(P); 0175 JUMPS = [0; diff(P)] < -1e-3; % modified! was 0! 0176 UP = P+cumsum(JUMPS)*pi; 0177 0178