source: MML/trunk/at/doc_html/at/atphysics/ohmienvelope.html @ 4

Last change on this file since 4 was 4, checked in by zhangj, 10 years ago

Initial import--MML version from SOLEIL@2013
File size: 9.9 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2                "http://www.w3.org/TR/REC-html40/loose.dtd">
3<html>
4<head>
5  <title>Description of ohmienvelope</title>
6  <meta name="keywords" content="ohmienvelope">
7  <meta name="description" content="OHMIENVELOPE calculates equilibrium beam envelope in a">
8  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
9  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
10  <meta name="robots" content="index, follow">
11  <link type="text/css" rel="stylesheet" href="../../m2html.css">
12</head>
13<body>
14<a name="_top"></a>
15<div><a href="../../index.html">Home</a> &gt;  <a href="../index.html">at</a> &gt; <a href="index.html">atphysics</a> &gt; ohmienvelope.m</div>
16
17<!--<table width="100%"><tr><td align="left"><a href="../../index.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
18<td align="right"><a href="index.html">Index for at/atphysics&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->
19
20<h1>ohmienvelope
21</h1>
22
23<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
24<div class="box"><strong>OHMIENVELOPE calculates equilibrium beam envelope in a</strong></div>
25
26<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
27<div class="box"><strong>function [ENVELOPE, RMSDP, RMSBL] = ohmienvelope(RING,RADELEMINDEX,varargin) </strong></div>
28
29<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
30<div class="fragment"><pre class="comment">OHMIENVELOPE calculates equilibrium beam envelope in a
31 circular accelerator using Ohmi's beam envelope formalism [1]
32 [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994)
33 
34 [ENVELOPE, RMSDP, RMSBL] = ONMIENVELOPE(RING,RADELEMINDEX,REFPTS)
35 
36 ENVELOPE is a structure with fields
37 Sigma   - [SIGMA(1); SIGMA(2)] - RMS size [m] along
38           the principal axis of a tilted ellips
39           Assuming normal distribution exp(-(Z^2)/(2*SIGMA))
40 Tilt    - Tilt angle of the XY ellips [rad]
41           Positive Tilt corresponds to Corkscrew (right)
42           rotatiom of XY plane around s-axis
43 R       - 6-by-6 equilibrium envelope matrix R
44
45 RMSDP   - RMS momentum spread
46 RMSBL   - RMS bunch length[m]</pre></div>
47
48<!-- crossreference -->
49<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
50This function calls:
51<ul style="list-style-image:url(../../matlabicon.gif)">
52<li><a href="findelemm66.html" class="code" title="function M66 = findelemm66(ELEM, MethodName, orbit_in);">findelemm66</a>     FINDELEMM66 numerically finds the 6x6 transfer matrix of an element</li><li><a href="findm66.html" class="code" title="function [M66, varargout] = findm66(RING, varargin);">findm66</a>        FINDM66 numerically finds the 6x6 transfer matrix of an accelerator lattice</li><li><a href="findmpoleraddiffmatrix.html" class="code" title="">findmpoleraddiffmatrix</a>      FINDMPOLERADDIFFMATIRX calculates radiation diffusion matrix of a multipole element</li></ul>
53This function is called by:
54<ul style="list-style-image:url(../../matlabicon.gif)">
55</ul>
56<!-- crossreference -->
57
58
59<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
60<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [ENVELOPE, RMSDP, RMSBL] = ohmienvelope(RING,RADELEMINDEX,varargin)</a>
610002 <span class="comment">%OHMIENVELOPE calculates equilibrium beam envelope in a</span>
620003 <span class="comment">% circular accelerator using Ohmi's beam envelope formalism [1]</span>
630004 <span class="comment">% [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994)</span>
640005 <span class="comment">%</span>
650006 <span class="comment">% [ENVELOPE, RMSDP, RMSBL] = ONMIENVELOPE(RING,RADELEMINDEX,REFPTS)</span>
660007 <span class="comment">%</span>
670008 <span class="comment">% ENVELOPE is a structure with fields</span>
680009 <span class="comment">% Sigma   - [SIGMA(1); SIGMA(2)] - RMS size [m] along</span>
690010 <span class="comment">%           the principal axis of a tilted ellips</span>
700011 <span class="comment">%           Assuming normal distribution exp(-(Z^2)/(2*SIGMA))</span>
710012 <span class="comment">% Tilt    - Tilt angle of the XY ellips [rad]</span>
720013 <span class="comment">%           Positive Tilt corresponds to Corkscrew (right)</span>
730014 <span class="comment">%           rotatiom of XY plane around s-axis</span>
740015 <span class="comment">% R       - 6-by-6 equilibrium envelope matrix R</span>
750016 <span class="comment">%</span>
760017 <span class="comment">% RMSDP   - RMS momentum spread</span>
770018 <span class="comment">% RMSBL   - RMS bunch length[m]</span>
780019
790020 NumElements = length(RING);
800021 
810022 [MRING, MS, orbit] = <a href="findm66.html" class="code" title="function [M66, varargout] = findm66(RING, varargin);">findm66</a>(RING,1:NumElements+1);
820023
830024 B = cell(1,NumElements); <span class="comment">% B{i} is the diffusion matrix of the i-th element</span>
840025 BATEXIT = B;             <span class="comment">% BATEXIT{i} is the cumulative diffusion matrix from</span>
850026                          <span class="comment">% element 0 to the end of the i-th element</span>
860027 M = B;                   <span class="comment">% 6-by-6 transfer matrix of the i-th element</span>
870028
880029 <span class="comment">% calculate Radiation-Diffusion matrix B for elements with radiation</span>
890030 <span class="keyword">for</span> i = RADELEMINDEX
900031    B{i} = <a href="findmpoleraddiffmatrix.html" class="code" title="">findmpoleraddiffmatrix</a>(RING{i},orbit(:,i));
910032 <span class="keyword">end</span>
920033
930034 <span class="comment">% Calculate 6-by-6 linear transfer matrix in each element</span>
940035 <span class="comment">% near the equilibrium orbit</span>
950036 <span class="keyword">for</span> i = 1:NumElements
960037    ELEM = RING{i};
970038    M{i} = <a href="findelemm66.html" class="code" title="function M66 = findelemm66(ELEM, MethodName, orbit_in);">findelemm66</a>(ELEM,ELEM.PassMethod,orbit(:,i));
980039    <span class="comment">% Set Radiation-Diffusion matrix B to 0 in elements without radiation</span>
990040    <span class="keyword">if</span> isempty(B{i})
1000041       B{i} = zeros(6,6);
1010042    <span class="keyword">end</span>
1020043 <span class="keyword">end</span>
1030044 <span class="comment">% Calculate cumulative Radiation-Diffusion matrix for the ring</span>
1040045 BCUM = zeros(6,6); <span class="comment">% Cumulative diffusion matrix of the entire ring</span>
1050046
1060047 <span class="keyword">for</span> i = 1:NumElements
1070048    BCUM = M{i}*BCUM*M{i}' + B{i};
1080049    BATEXIT{i} = BCUM;
1090050 <span class="keyword">end</span>
1100051 <span class="comment">% ------------------------------------------------------------------------</span>
1110052 <span class="comment">% Equation for the moment matrix R is</span>
1120053 <span class="comment">%         R = MRING*R*MRING' + BCUM;</span>
1130054 <span class="comment">% We rewrite it in the form of Lyapunov equation to use MATLAB's LYAP function</span>
1140055 <span class="comment">%            AA*R + R*BB = -CC</span>
1150056 <span class="comment">% where</span>
1160057 <span class="comment">%                AA =  inv(MRING)</span>
1170058 <span class="comment">%                BB = -MRING'</span>
1180059 <span class="comment">%                CC = -inv(MRING)*BCUM</span>
1190060 <span class="comment">% -----------------------------------------------------------------------</span>
1200061 AA =  inv(MRING);
1210062 BB = -MRING';
1220063 CC = -inv(MRING)*BCUM;
1230064 
1240065 R = lyap(AA,BB,CC);     <span class="comment">% Envelope matrix at the rinng entrance</span>
1250066
1260067 RMSDP = sqrt(R(5,5));   <span class="comment">% R.M.S. energy spread</span>
1270068 RMSBL = sqrt(R(6,6));   <span class="comment">% R.M.S. bunch lenght</span>
1280069
1290070 <span class="keyword">if</span> nargin == 2 <span class="comment">% no REFPTS</span>
1300071     ENVELOPE.R = R;
1310072 <span class="keyword">elseif</span> nargin == 3
1320073     REFPTS = varargin{1};
1330074     
1340075     REFPTSX = REFPTS;
1350076     <span class="comment">% complete REFPTS with 1 and NumElements+1 if necessary</span>
1360077     <span class="keyword">if</span> REFPTS(1)~=1
1370078         REFPTSX = [1 REFPTS];
1380079     <span class="keyword">end</span>
1390080     <span class="keyword">if</span> REFPTS(end)~= NumElements+1
1400081         REFPTSX = [REFPTSX NumElements+1];
1410082     <span class="keyword">end</span>
1420083     <span class="comment">% Now REFPTS has at least 2 ponts and the first one is the ring entrance</span>
1430084     
1440085     NRX = length(REFPTSX);
1450086     ENVELOPE = struct(<span class="string">'Sigma'</span>,num2cell(zeros(2,NRX),1),<span class="keyword">...</span>
1460087         <span class="string">'Tilt'</span>,0,<span class="string">'R'</span>,zeros(6));
1470088     
1480089     ENVELOPE(1).R = R;
1490090
1500091     <span class="keyword">for</span> i=2:NRX
1510092         ELEM = REFPTSX(i);
1520093         ENVELOPE(i).R = MS(:,:,ELEM)*R*MS(:,:,ELEM)'+BATEXIT{ELEM-1};
1530094     <span class="keyword">end</span>
1540095     
1550096   
1560097     <span class="keyword">if</span> REFPTS(1)~=1
1570098         ENVELOPE = ENVELOPE(2:end);
1580099     <span class="keyword">end</span>
1590100     <span class="keyword">if</span> REFPTS(end)~= NumElements+1
1600101         ENVELOPE = ENVELOPE(1:end-1);
1610102     <span class="keyword">end</span>
1620103
1630104 <span class="keyword">else</span>
1640105     error(<span class="string">'Too many input arguments'</span>);
1650106 <span class="keyword">end</span>
1660107
1670108 <span class="keyword">for</span> i=1:length(ENVELOPE)
1680109     [U,DR] = eig(ENVELOPE(i).R([1 3],[1 3]));
1690110     ENVELOPE(i).Tilt = asin((U(2,1)-U(1,2))/2);
1700111     ENVELOPE(i).Sigma(1) = sqrt(DR(1,1));
1710112     ENVELOPE(i).Sigma(2) = sqrt(DR(2,2));
1720113 <span class="keyword">end</span></pre></div>
173<hr><address>Generated on Mon 21-May-2007 15:26:45 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
174</body>
175</html>
Note: See TracBrowser for help on using the repository browser.