TWISSRING calculates linear optics functions for an UNCOUPLED ring
[TwissData, tune] = TWISSRING(LATTICE,DP) calculates twiss parameters
and closed orbit coordinates at the RING entrance assuming
constant energy deviation DP.
[TwissData, tune] = TWISSRING(LATTICE,DP,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, tune, chrom] = TWISSRING(...,'chrom', DDP) also calculates
linear dispersion and chromaticity. 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 and chromaticity. If not supplied
the default value of 1e-8 is used.
Note: To resolve the integer part of the tune
and the uncertainty of acos(trace(M)/2) it is necessary to
supply sufficient number of REFPTS properly spaced in betatron phase.
TwisData is a 1-by-REFPTS (1-by-1) structure array with fields
(Some are the same as in the output of LINOPT)
ElemIndex - integer (element number) 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 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 RING
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 MATLAB function 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 TWISSLINE, LINOPT, TUNECHROM.0001 function [TD, varargout] = twissring(RING,DP,varargin) 0002 %TWISSRING calculates linear optics functions for an UNCOUPLED ring 0003 % 0004 % [TwissData, tune] = TWISSRING(LATTICE,DP) calculates twiss parameters 0005 % and closed orbit coordinates at the RING entrance assuming 0006 % constant energy deviation DP. 0007 % 0008 % [TwissData, tune] = TWISSRING(LATTICE,DP,REFPTS) calculates Twiss parameters 0009 % and closed orbit coordinates at specified reference points REFPTS. 0010 % 0011 % Note: REFPTS is an array of increasing indexes that 0012 % select elements from range 1 to length(LATTICE)+1. 0013 % See further explanation of REFPTS in the 'help' for FINDSPOS 0014 % 0015 % [TwissData, tune, chrom] = TWISSRING(...,'chrom', DDP) also calculates 0016 % linear dispersion and chromaticity. Dispersion is returned as one 0017 % of the fields in TwissData. 0018 % !!! Last argument DDP is a momentum deviation on top 0019 % of DP (the second argument) used to calculate and normalize 0020 % dispersion and chromaticity. If not supplied 0021 % the default value of 1e-8 is used. 0022 % 0023 % Note: To resolve the integer part of the tune 0024 % and the uncertainty of acos(trace(M)/2) it is necessary to 0025 % supply sufficient number of REFPTS properly spaced in betatron phase. 0026 % 0027 % TwisData is a 1-by-REFPTS (1-by-1) structure array with fields 0028 % (Some are the same as in the output of LINOPT) 0029 % ElemIndex - integer (element number) in the RING 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 RING 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 MATLAB function 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 TWISSLINE, LINOPT, TUNECHROM. 0052 0053 NE=length(RING); 0054 DDP_default = 1e-8; 0055 % Process input arguments 0056 switch nargin 0057 case 2 0058 REFPTS=NE+1; 0059 CHROMFLAG=0; 0060 case 3 0061 if isnumeric(varargin{1}) 0062 REFPTS = varargin{1}; 0063 CHROMFLAG = 0; 0064 elseif ischar(varargin{1}) & strncmp(lower(varargin{1}),'chrom',5) 0065 CHROMFLAG = 1; 0066 REFPTS = NE+1; 0067 DDP = DDP_default; 0068 else 0069 error('Third argument must be a numeric array or string'); 0070 end 0071 case 4 0072 if isnumeric(varargin{1}) 0073 REFPTS = varargin{1}; 0074 if ischar(varargin{2}) & strncmp(lower(varargin{2}),'chrom',5) 0075 CHROMFLAG = 1; 0076 DDP = DDP_default; 0077 else 0078 error('Fourth argument - wrong type'); 0079 end 0080 elseif ischar(varargin{1}) & strncmp(lower(varargin{1}),'chrom',5) 0081 CHROMFLAG = 1; 0082 REFPTS = NE+1; 0083 if isnumeric(varargin{2}) 0084 DDP = varargin{2}; 0085 else 0086 error('Fourth argument - wrong type'); 0087 end 0088 end 0089 case 5 0090 if isnumeric(varargin{1}) 0091 REFPTS = varargin{1}; 0092 else 0093 error('Third argument - wrong type'); 0094 end 0095 if ischar(varargin{2}) & strncmp(lower(varargin{2}),'chrom',5) 0096 CHROMFLAG = 1; 0097 else 0098 error('Fourth argument - wrong type'); 0099 end 0100 if isnumeric(varargin{3}) 0101 DDP = varargin{3}; 0102 else 0103 error('Fifth argument - wrong type'); 0104 end 0105 otherwise 0106 error('Wrong number of arguments'); 0107 end 0108 0109 % Include the endpoint if it is not already in REFPTS 0110 if REFPTS(end)==NE+1 0111 [M44, MS, orb] = findm44(RING,DP,REFPTS); 0112 else 0113 [M44, MS, orb] = findm44(RING,DP,[REFPTS,NE+1]); 0114 end 0115 0116 0117 0118 0119 cos_mu_x = (M44(1,1)+M44(2,2))/2; 0120 cos_mu_y = (M44(3,3)+M44(4,4))/2; 0121 0122 sin_mu_x = sign(M44(1,2))*sqrt(-M44(1,2)*M44(2,1)-(M44(1,1)-M44(2,2))^2/4); 0123 sin_mu_y = sign(M44(3,4))*sqrt(-M44(3,4)*M44(4,3)-(M44(3,3)-M44(4,4))^2/4); 0124 0125 0126 ax = (M44(1,1)-M44(2,2))/2/sin_mu_x; 0127 ay = (M44(3,3)-M44(4,4))/2/sin_mu_y; 0128 0129 bx = M44(1,2)/sin_mu_x; 0130 by = M44(3,4)/sin_mu_y; 0131 0132 BX = squeeze((MS(1,1,:)*bx-MS(1,2,:)*ax).^2 + MS(1,2,:).^2)/bx; 0133 BY = squeeze((MS(3,3,:)*by-MS(3,4,:)*ay).^2 + MS(3,4,:).^2)/by; 0134 0135 0136 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; 0137 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; 0138 0139 MX = atan(squeeze( MS(1,2,:)./(MS(1,1,:)*bx-MS(1,2,:)*ax))); 0140 MY = atan(squeeze(MS(3,4,:)./(MS(3,3,:)*by-MS(3,4,:)*ay))); 0141 0142 MX = BetatronPhaseUnwrap(MX); 0143 MY = BetatronPhaseUnwrap(MY); 0144 0145 tune = [MX(end),MY(end)]/2/pi; 0146 0147 NR = length(REFPTS); 0148 % Build TD only for points originally referenced in REFPTS 0149 TD = struct('ElemIndex',num2cell(REFPTS),... 0150 'SPos',num2cell(findspos(RING,REFPTS)),... 0151 'ClosedOrbit',num2cell(orb(:,1:NR),1),... 0152 'M44', squeeze(num2cell(MS(:,:,1:NR),[1 2]))',... 0153 'beta', num2cell([BX(1:NR),BY(1:NR)],2)',... 0154 'alpha', num2cell([AX(1:NR),AY(1:NR)],2)',... 0155 'mu', num2cell([MX(1:NR),MY(1:NR)],2)'); 0156 0157 0158 if CHROMFLAG 0159 [TD_DDP tune_DDP] = twissring(RING,DP+DDP,REFPTS); 0160 DORBIT = reshape(cat(1,TD_DDP.ClosedOrbit),4,[]); 0161 DISPERSION = num2cell((DORBIT-orb(:,1:NR))/DDP,1); 0162 [TD.Dispersion] = deal( DISPERSION{:}); 0163 end 0164 0165 if nargout>1 0166 varargout{1}=tune; 0167 end 0168 if nargout==3 & CHROMFLAG 0169 varargout{2} = (tune_DDP-tune)/DDP; 0170 end 0171 0172 function UP = BetatronPhaseUnwrap(P) 0173 % unwrap negative jumps in betatron 0174 DP = diff(P); 0175 JUMPS = [0; diff(P)] < 0; 0176 UP = P+cumsum(JUMPS)*pi; 0177 0178