source: MML/trunk/at/doc_html/at/atphysics/findsyncorbit.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: 11.4 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 findsyncorbit</title>
6  <meta name="keywords" content="findsyncorbit">
7  <meta name="description" content="FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity">
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; findsyncorbit.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>findsyncorbit
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>FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity</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 [orbit, varargout] = findsyncorbit(RING,dCT,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">FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity
31 and momentum deviation dP (first 5 components of the phase space vector)
32 by numerically solving  for a fixed point
33 of the one turn map M calculated with LINEPASS
34
35       (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
36   
37    under constraints CT2 - CT1 =  dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where
38    Frev0 = Frf0/HarmNumber is the design revolution frequency
39    Frev  = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency
40
41 IMPORTANT!!!  FINDSYNCORBIT imposes a constraint (CT2 - CT1) and
42    dP2 = dP1 but no constraint on the value of dP1, dP2
43    The algorithm assumes time-independent fixed-momentum ring
44    to reduce the dimensionality of the problem.
45    To impose this artificial constraint in FINDSYNCORBIT
46    PassMethod used for any element SHOULD NOT
47    1. change the longitudinal momentum dP (cavities , magnets with radiation)
48    2. have any time dependence (localized impedance, fast kickers etc).
49
50
51 FINDSYNCORBIT(RING,dCT) is 5x1 vector - fixed point at the
52        entrance of the 1-st element of the RING (x,px,y,py,dP)
53
54 FINDSYNCORBIT(RING,dCT,REFPTS) is 5-by-Length(REFPTS)
55     array of column vectors - fixed points (x,px,y,py,dP)
56     at the entrance of each element indexed in REFPTS array.
57     REFPTS is an array of increasing indexes that  select elements
58     from the range 1 to length(RING)+1.
59     See further explanation of REFPTS in the 'help' for FINDSPOS
60
61 FINDSYNCORBIT(RING,dCT,REFPTS,GUESS) - same as above but the search
62     for the fixed point starts at the initial condition GUESS
63     Otherwise the search starts from [0 0 0 0 0 0]'.
64     GUESS must be a 6-by-1 vector;
65
66 [ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... )
67     The optional second return parameter is
68     a 6-by-1 vector of initial conditions
69     on the close orbit at the entrance to the RING. 
70
71 See also <a href="findorbit.html" class="code" title="function [orbit, varargout]  = findorbit(RING,D, varargin);">FINDORBIT</a>, <a href="findorbit4.html" class="code" title="function orbit = findorbit4(RING,dP,varargin);">FINDORBIT4</a>, <a href="findorbit6.html" class="code" title="function orbit = findorbit6(RING,varargin);">FINDORBIT6</a>.</pre></div>
72
73<!-- crossreference -->
74<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
75This function calls:
76<ul style="list-style-image:url(../../matlabicon.gif)">
77</ul>
78This function is called by:
79<ul style="list-style-image:url(../../matlabicon.gif)">
80<li><a href="findrespm.html" class="code" title="function C = findrespm(RING, OBSINDEX, PERTURB, PVALUE, varargin)">findrespm</a>       FINDRESPM computes the change in the closed orbit due to parameter perturbations</li></ul>
81<!-- crossreference -->
82
83
84<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
85<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [orbit, varargout] = findsyncorbit(RING,dCT,varargin);</a>
860002 <span class="comment">%FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity</span>
870003 <span class="comment">% and momentum deviation dP (first 5 components of the phase space vector)</span>
880004 <span class="comment">% by numerically solving  for a fixed point</span>
890005 <span class="comment">% of the one turn map M calculated with LINEPASS</span>
900006 <span class="comment">%</span>
910007 <span class="comment">%       (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)</span>
920008 <span class="comment">%</span>
930009 <span class="comment">%    under constraints CT2 - CT1 =  dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where</span>
940010 <span class="comment">%    Frev0 = Frf0/HarmNumber is the design revolution frequency</span>
950011 <span class="comment">%    Frev  = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency</span>
960012 <span class="comment">%</span>
970013 <span class="comment">% IMPORTANT!!!  FINDSYNCORBIT imposes a constraint (CT2 - CT1) and</span>
980014 <span class="comment">%    dP2 = dP1 but no constraint on the value of dP1, dP2</span>
990015 <span class="comment">%    The algorithm assumes time-independent fixed-momentum ring</span>
1000016 <span class="comment">%    to reduce the dimensionality of the problem.</span>
1010017 <span class="comment">%    To impose this artificial constraint in FINDSYNCORBIT</span>
1020018 <span class="comment">%    PassMethod used for any element SHOULD NOT</span>
1030019 <span class="comment">%    1. change the longitudinal momentum dP (cavities , magnets with radiation)</span>
1040020 <span class="comment">%    2. have any time dependence (localized impedance, fast kickers etc).</span>
1050021 <span class="comment">%</span>
1060022 <span class="comment">%</span>
1070023 <span class="comment">% FINDSYNCORBIT(RING,dCT) is 5x1 vector - fixed point at the</span>
1080024 <span class="comment">%        entrance of the 1-st element of the RING (x,px,y,py,dP)</span>
1090025 <span class="comment">%</span>
1100026 <span class="comment">% FINDSYNCORBIT(RING,dCT,REFPTS) is 5-by-Length(REFPTS)</span>
1110027 <span class="comment">%     array of column vectors - fixed points (x,px,y,py,dP)</span>
1120028 <span class="comment">%     at the entrance of each element indexed in REFPTS array.</span>
1130029 <span class="comment">%     REFPTS is an array of increasing indexes that  select elements</span>
1140030 <span class="comment">%     from the range 1 to length(RING)+1.</span>
1150031 <span class="comment">%     See further explanation of REFPTS in the 'help' for FINDSPOS</span>
1160032 <span class="comment">%</span>
1170033 <span class="comment">% FINDSYNCORBIT(RING,dCT,REFPTS,GUESS) - same as above but the search</span>
1180034 <span class="comment">%     for the fixed point starts at the initial condition GUESS</span>
1190035 <span class="comment">%     Otherwise the search starts from [0 0 0 0 0 0]'.</span>
1200036 <span class="comment">%     GUESS must be a 6-by-1 vector;</span>
1210037 <span class="comment">%</span>
1220038 <span class="comment">% [ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... )</span>
1230039 <span class="comment">%     The optional second return parameter is</span>
1240040 <span class="comment">%     a 6-by-1 vector of initial conditions</span>
1250041 <span class="comment">%     on the close orbit at the entrance to the RING.</span>
1260042 <span class="comment">%</span>
1270043 <span class="comment">% See also FINDORBIT, FINDORBIT4, FINDORBIT6.</span>
1280044 <span class="comment">%</span>
1290045 <span class="keyword">if</span> ~iscell(RING)
1300046    error(<span class="string">'First argument must be a cell array'</span>);
1310047 <span class="keyword">end</span>
1320048
1330049 d = 1e-9;    <span class="comment">% step size for numerical differentiation</span>
1340050 max_iterations = 20;
1350051 J = zeros(5);
1360052
1370053 <span class="keyword">if</span> nargin==4
1380054     <span class="keyword">if</span> isnumeric(varargin{2}) &amp; all(eq(size(varargin{2}),[6,1]))
1390055        Ri=varargin{2};
1400056    <span class="keyword">else</span>
1410057        error(<span class="string">'The last argument GUESS must be a 6-by-1 vector'</span>);
1420058    <span class="keyword">end</span>
1430059 <span class="keyword">else</span>
1440060     Ri = zeros(6,1);
1450061 <span class="keyword">end</span>
1460062
1470063 D5 = d*eye(5); 
1480064 <span class="comment">%B5 = zeros(5,5);</span>
1490065 RMATi = zeros(6,6);
1500066 theta5 = [0 0 0 0  dCT]';
1510067
1520068
1530069
1540070 <span class="comment">%Fist iteration</span>
1550071 RMATi= Ri*ones(1,6) + [D5 zeros(5,1); zeros(1,6)];
1560072 RMATf = linepass(RING,RMATi);
1570073 Rf = RMATf(:,6);
1580074 <span class="comment">% compute the transverse part of the Jacobian</span>
1590075 J5 =  [RMATf([1:4,6],1:5)-RMATf([1:4,6],6)*ones(1,5)]/d;
1600076
1610077 <span class="comment">% Replace matrix inversion with \</span>
1620078 <span class="comment">%B5 = inv(diag([1 1 1 1 0]) - J5);</span>
1630079 <span class="comment">% Ri_next = Ri +  [B5* (Rf([1:4,6])-[Ri(1:4);0]-theta5) ; 0];</span>
1640080 Ri_next = Ri +  [(diag([1 1 1 1 0]) - J5)\(Rf([1:4,6])-[Ri(1:4);0]-theta5) ; 0];
1650081 change = norm(Ri_next - Ri);
1660082 Ri = Ri_next;
1670083
1680084 itercount = 1;
1690085
1700086
1710087 <span class="keyword">while</span> (change&gt;eps) &amp; (itercount &lt; max_iterations)
1720088   
1730089    RMATi= Ri*ones(1,6) + [D5 zeros(5,1); zeros(1,6)];   
1740090    RMATf = linepass(RING,RMATi,<span class="string">'reuse'</span>);
1750091
1760092    Rf = RMATf(:,6);
1770093    <span class="comment">% compute the transverse part of the Jacobian</span>
1780094    J5 =  [RMATf([1:4,6],1:5)-RMATf([1:4,6],6)*ones(1,5)]/d;
1790095    <span class="comment">% Replace matrix inversion with \</span>
1800096    <span class="comment">%B5 = inv(diag([1 1 1 1 0]) - J5);</span>
1810097    <span class="comment">%Ri_next = Ri +  [B5*(Rf([1:4,6])-[Ri(1:4);0]-theta5); 0];</span>
1820098    Ri_next = Ri +  [(diag([1 1 1 1 0]) - J5)\(Rf([1:4,6])-[Ri(1:4);0]-theta5); 0];
1830099    change = norm(Ri_next - Ri);
1840100    Ri = Ri_next;
1850101    itercount = itercount+1;
1860102   
1870103 <span class="keyword">end</span>;
1880104
1890105
1900106 <span class="keyword">if</span>(nargin&lt;3) | (varargin{1}==(length(RING)+1))
1910107     <span class="comment">% return only the fixed point at the entrance of RING{1}</span>
1920108     orbit=Ri(1:5,1);
1930109 <span class="keyword">else</span>            <span class="comment">% 3-rd input argument - vector of reference points along the Ring</span>
1940110                 <span class="comment">% is supplied - return orbit</span>
1950111    orb6 = linepass(RING,Ri,varargin{1},<span class="string">'reuse'</span>);
1960112    orbit = orb6(1:5,:);
1970113 <span class="keyword">end</span>
1980114
1990115 <span class="keyword">if</span> nargout==2
2000116     varargout{1}=Ri;
2010117 <span class="keyword">end</span></pre></div>
202<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>
203</body>
204</html>
Note: See TracBrowser for help on using the repository browser.