OHMIENVELOPE calculates equilibrium beam envelope in a
circular accelerator using Ohmi's beam envelope formalism [1]
[1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994)
[ENVELOPE, RMSDP, RMSBL] = ONMIENVELOPE(RING,RADELEMINDEX,REFPTS)
ENVELOPE is a structure with fields
Sigma - [SIGMA(1); SIGMA(2)] - RMS size [m] along
the principal axis of a tilted ellips
Assuming normal distribution exp(-(Z^2)/(2*SIGMA))
Tilt - Tilt angle of the XY ellips [rad]
Positive Tilt corresponds to Corkscrew (right)
rotatiom of XY plane around s-axis
R - 6-by-6 equilibrium envelope matrix R
RMSDP - RMS momentum spread
RMSBL - RMS bunch length[m]
[ENVELOPE, RMSDP, RMSBL, M66, T, ORBIT] = OHMIENVELOPE(...)
Returns in addition the 6x6 transfer matrices and the closed orbit
from FINDM660001 function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts) 0002 %OHMIENVELOPE calculates equilibrium beam envelope in a 0003 % circular accelerator using Ohmi's beam envelope formalism [1] 0004 % [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994) 0005 % 0006 % [ENVELOPE, RMSDP, RMSBL] = ONMIENVELOPE(RING,RADELEMINDEX,REFPTS) 0007 % 0008 % ENVELOPE is a structure with fields 0009 % Sigma - [SIGMA(1); SIGMA(2)] - RMS size [m] along 0010 % the principal axis of a tilted ellips 0011 % Assuming normal distribution exp(-(Z^2)/(2*SIGMA)) 0012 % Tilt - Tilt angle of the XY ellips [rad] 0013 % Positive Tilt corresponds to Corkscrew (right) 0014 % rotatiom of XY plane around s-axis 0015 % R - 6-by-6 equilibrium envelope matrix R 0016 % 0017 % RMSDP - RMS momentum spread 0018 % RMSBL - RMS bunch length[m] 0019 % 0020 % [ENVELOPE, RMSDP, RMSBL, M66, T, ORBIT] = OHMIENVELOPE(...) 0021 % Returns in addition the 6x6 transfer matrices and the closed orbit 0022 % from FINDM66 0023 0024 NumElements = length(ring); 0025 if nargin<3, refpts=1; end 0026 0027 [mring, ms, orbit] = findm66(ring,1:NumElements+1); 0028 mt=squeeze(num2cell(ms,[1 2])); 0029 orb=num2cell(orbit,1)'; 0030 0031 zr={zeros(6,6)}; 0032 B=zr(ones(NumElements,1)); % B{i} is the diffusion matrix of the i-th element 0033 0034 % calculate Radiation-Diffusion matrix B for elements with radiation 0035 B(radindex)=cellfun(@findmpoleraddiffmatrix,... 0036 ring(radindex),orb(radindex),'UniformOutput',false); 0037 0038 % Calculate cumulative Radiation-Diffusion matrix for the ring 0039 BCUM = zeros(6,6); 0040 % Batbeg{i} is the cumulative diffusion matrix from 0041 % 0 to the beginning of the i-th element 0042 Batbeg=[zr;cellfun(@cumulb,ring,orb(1:end-1),B,'UniformOutput',false)]; 0043 0044 % ------------------------------------------------------------------------ 0045 % Equation for the moment matrix R is 0046 % R = MRING*R*MRING' + BCUM; 0047 % We rewrite it in the form of Lyapunov equation to use MATLAB's LYAP function 0048 % AA*R + R*BB = -CC 0049 % where 0050 % AA = inv(MRING) 0051 % BB = -MRING' 0052 % CC = -inv(MRING)*BCUM 0053 % ----------------------------------------------------------------------- 0054 AA = inv(mring); 0055 BB = -mring'; 0056 CC = -inv(mring)*BCUM; 0057 0058 R = lyap(AA,BB,CC); % Envelope matrix at the ring entrance 0059 0060 rmsdp = sqrt(R(5,5)); % R.M.S. energy spread 0061 rmsbl = sqrt(R(6,6)); % R.M.S. bunch lenght 0062 0063 [rr,tt,ss]=cellfun(@propag,mt(refpts),Batbeg(refpts),'UniformOutput',false); 0064 envelope=struct('R',rr,'Sigma',ss,'Tilt',tt); 0065 0066 nout=nargout-3; 0067 varargout=cell(1,nout); 0068 if nout>=3, varargout{3}=orbit(:,refpts); end 0069 if nout>=2, varargout{2}=ms(:,:,refpts); end 0070 if nout>=1, varargout{1}=mring; end 0071 0072 function btx=cumulb(elem,orbit,b) 0073 % Calculate 6-by-6 linear transfer matrix in each element 0074 % near the equilibrium orbit 0075 m=findelemm66(elem,elem.PassMethod,orbit); 0076 % Cumulative diffusion matrix of the entire ring 0077 BCUM = m*BCUM*m' + b; 0078 btx=BCUM; 0079 end 0080 0081 function [r,tilt,sigma]=propag(m,cumb) 0082 r=m*R*m'+cumb; 0083 [u,dr] = eig(r([1 3],[1 3])); 0084 tilt = asin((u(2,1)-u(1,2))/2); 0085 sigma=sqrt([dr(1,1);dr(2,2)]); 0086 end 0087 0088 end