[4] | 1 | function [ENVELOPE, RMSDP, RMSBL] = ohmienvelope(RING,RADELEMINDEX,varargin) |
---|
| 2 | %OHMIENVELOPE calculates equilibrium beam envelope in a |
---|
| 3 | % circular accelerator using Ohmi's beam envelope formalism [1] |
---|
| 4 | % [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994) |
---|
| 5 | % |
---|
| 6 | % [ENVELOPE, RMSDP, RMSBL] = ONMIENVELOPE(RING,RADELEMINDEX,REFPTS) |
---|
| 7 | % |
---|
| 8 | % ENVELOPE is a structure with fields |
---|
| 9 | % Sigma - [SIGMA(1); SIGMA(2)] - RMS size [m] along |
---|
| 10 | % the principal axis of a tilted ellips |
---|
| 11 | % Assuming normal distribution exp(-(Z^2)/(2*SIGMA)) |
---|
| 12 | % Tilt - Tilt angle of the XY ellips [rad] |
---|
| 13 | % Positive Tilt corresponds to Corkscrew (right) |
---|
| 14 | % rotation of XY plane around s-axis |
---|
| 15 | % R - 6-by-6 equilibrium envelope matrix R |
---|
| 16 | % |
---|
| 17 | % RMSDP - RMS momentum spread |
---|
| 18 | % RMSBL - RMS bunch length[m] |
---|
| 19 | |
---|
| 20 | NumElements = length(RING); |
---|
| 21 | |
---|
| 22 | [MRING, MS, orbit] = findm66(RING,1:NumElements+1); |
---|
| 23 | |
---|
| 24 | B = cell(1,NumElements); % B{i} is the diffusion matrix of the i-th element |
---|
| 25 | BATEXIT = B; % BATEXIT{i} is the cumulative diffusion matrix from |
---|
| 26 | % element 0 to the end of the i-th element |
---|
| 27 | M = B; % 6-by-6 transfer matrix of the i-th element |
---|
| 28 | |
---|
| 29 | % calculate Radiation-Diffusion matrix B for elements with radiation |
---|
| 30 | for i = RADELEMINDEX |
---|
| 31 | B{i} = findmpoleraddiffmatrix(RING{i},orbit(:,i)); |
---|
| 32 | end |
---|
| 33 | |
---|
| 34 | % Calculate 6-by-6 linear transfer matrix in each element |
---|
| 35 | % near the equilibrium orbit |
---|
| 36 | for i = 1:NumElements |
---|
| 37 | ELEM = RING{i}; |
---|
| 38 | M{i} = findelemm66(ELEM,ELEM.PassMethod,orbit(:,i)); |
---|
| 39 | % Set Radiation-Diffusion matrix B to 0 in elements without radiation |
---|
| 40 | if isempty(B{i}) |
---|
| 41 | B{i} = zeros(6,6); |
---|
| 42 | end |
---|
| 43 | end |
---|
| 44 | % Calculate cumulative Radiation-Diffusion matrix for the ring |
---|
| 45 | BCUM = zeros(6,6); % Cumulative diffusion matrix of the entire ring |
---|
| 46 | |
---|
| 47 | for i = 1:NumElements |
---|
| 48 | BCUM = M{i}*BCUM*M{i}' + B{i}; |
---|
| 49 | BATEXIT{i} = BCUM; |
---|
| 50 | end |
---|
| 51 | % ------------------------------------------------------------------------ |
---|
| 52 | % Equation for the moment matrix R is |
---|
| 53 | % R = MRING*R*MRING' + BCUM; |
---|
| 54 | % We rewrite it in the form of Lyapunov equation to use MATLAB's LYAP function |
---|
| 55 | % AA*R + R*BB = -CC |
---|
| 56 | % where |
---|
| 57 | % AA = inv(MRING) |
---|
| 58 | % BB = -MRING' |
---|
| 59 | % CC = -inv(MRING)*BCUM |
---|
| 60 | % ----------------------------------------------------------------------- |
---|
| 61 | AA = inv(MRING); |
---|
| 62 | BB = -MRING'; |
---|
| 63 | CC = -inv(MRING)*BCUM; |
---|
| 64 | |
---|
| 65 | R = lyap(AA,BB,CC); % Envelope matrix at the ring entrance |
---|
| 66 | |
---|
| 67 | RMSDP = sqrt(R(5,5)); % R.M.S. energy spread |
---|
| 68 | RMSBL = sqrt(R(6,6)); % R.M.S. bunch lenght |
---|
| 69 | |
---|
| 70 | if nargin == 2 % no REFPTS |
---|
| 71 | ENVELOPE.R = R; |
---|
| 72 | elseif nargin == 3 |
---|
| 73 | REFPTS = varargin{1}; |
---|
| 74 | |
---|
| 75 | REFPTSX = REFPTS; |
---|
| 76 | % complete REFPTS with 1 and NumElements+1 if necessary |
---|
| 77 | if REFPTS(1)~=1 |
---|
| 78 | REFPTSX = [1 REFPTS]; |
---|
| 79 | end |
---|
| 80 | if REFPTS(end)~= NumElements+1 |
---|
| 81 | REFPTSX = [REFPTSX NumElements+1]; |
---|
| 82 | end |
---|
| 83 | % Now REFPTS has at least 2 ponts and the first one is the ring entrance |
---|
| 84 | |
---|
| 85 | NRX = length(REFPTSX); |
---|
| 86 | ENVELOPE = struct('Sigma',num2cell(zeros(2,NRX),1),... |
---|
| 87 | 'Tilt',0,'R',zeros(6)); |
---|
| 88 | |
---|
| 89 | ENVELOPE(1).R = R; |
---|
| 90 | |
---|
| 91 | for i=2:NRX |
---|
| 92 | ELEM = REFPTSX(i); |
---|
| 93 | ENVELOPE(i).R = MS(:,:,ELEM)*R*MS(:,:,ELEM)'+BATEXIT{ELEM-1}; |
---|
| 94 | end |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | if REFPTS(1)~=1 |
---|
| 98 | ENVELOPE = ENVELOPE(2:end); |
---|
| 99 | end |
---|
| 100 | if REFPTS(end)~= NumElements+1 |
---|
| 101 | ENVELOPE = ENVELOPE(1:end-1); |
---|
| 102 | end |
---|
| 103 | |
---|
| 104 | else |
---|
| 105 | error('Too many input arguments'); |
---|
| 106 | end |
---|
| 107 | |
---|
| 108 | for i=1:length(ENVELOPE) |
---|
| 109 | [U,DR] = eig(ENVELOPE(i).R([1 3],[1 3])); |
---|
| 110 | ENVELOPE(i).Tilt = asin((U(2,1)-U(1,2))/2); |
---|
| 111 | ENVELOPE(i).Sigma(1) = sqrt(DR(1,1)); |
---|
| 112 | ENVELOPE(i).Sigma(2) = sqrt(DR(2,2)); |
---|
| 113 | end |
---|