1 | function [orbit, varargout] = findsyncorbit(RING,dCT,varargin); |
---|
2 | %FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity |
---|
3 | % and momentum deviation dP (first 5 components of the phase space vector) |
---|
4 | % by numerically solving for a fixed point |
---|
5 | % of the one turn map M calculated with LINEPASS |
---|
6 | % |
---|
7 | % (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1) |
---|
8 | % |
---|
9 | % under constraints CT2 - CT1 = dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where |
---|
10 | % Frev0 = Frf0/HarmNumber is the design revolution frequency |
---|
11 | % Frev = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency |
---|
12 | % |
---|
13 | % IMPORTANT!!! FINDSYNCORBIT imposes a constraint (CT2 - CT1) and |
---|
14 | % dP2 = dP1 but no constraint on the value of dP1, dP2 |
---|
15 | % The algorithm assumes time-independent fixed-momentum ring |
---|
16 | % to reduce the dimensionality of the problem. |
---|
17 | % To impose this artificial constraint in FINDSYNCORBIT |
---|
18 | % PassMethod used for any element SHOULD NOT |
---|
19 | % 1. change the longitudinal momentum dP (cavities , magnets with radiation) |
---|
20 | % 2. have any time dependence (localized impedance, fast kickers etc). |
---|
21 | % |
---|
22 | % |
---|
23 | % FINDSYNCORBIT(RING,dCT) is 5x1 vector - fixed point at the |
---|
24 | % entrance of the 1-st element of the RING (x,px,y,py,dP) |
---|
25 | % |
---|
26 | % FINDSYNCORBIT(RING,dCT,REFPTS) is 5-by-Length(REFPTS) |
---|
27 | % array of column vectors - fixed points (x,px,y,py,dP) |
---|
28 | % at the entrance of each element indexed in REFPTS array. |
---|
29 | % REFPTS is an array of increasing indexes that select elements |
---|
30 | % from the range 1 to length(RING)+1. |
---|
31 | % See further explanation of REFPTS in the 'help' for FINDSPOS |
---|
32 | % |
---|
33 | % FINDSYNCORBIT(RING,dCT,REFPTS,GUESS) - same as above but the search |
---|
34 | % for the fixed point starts at the initial condition GUESS |
---|
35 | % Otherwise the search starts from [0 0 0 0 0 0]'. |
---|
36 | % GUESS must be a 6-by-1 vector; |
---|
37 | % |
---|
38 | % [ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... ) |
---|
39 | % The optional second return parameter is |
---|
40 | % a 6-by-1 vector of initial conditions |
---|
41 | % on the close orbit at the entrance to the RING. |
---|
42 | % |
---|
43 | % See also FINDORBIT, FINDORBIT4, FINDORBIT6. |
---|
44 | % |
---|
45 | if ~iscell(RING) |
---|
46 | error('First argument must be a cell array'); |
---|
47 | end |
---|
48 | |
---|
49 | d = 1e-9; % step size for numerical differentiation |
---|
50 | max_iterations = 20; |
---|
51 | J = zeros(5); |
---|
52 | |
---|
53 | if nargin==4 |
---|
54 | if isnumeric(varargin{2}) & all(eq(size(varargin{2}),[6,1])) |
---|
55 | Ri=varargin{2}; |
---|
56 | else |
---|
57 | error('The last argument GUESS must be a 6-by-1 vector'); |
---|
58 | end |
---|
59 | else |
---|
60 | Ri = zeros(6,1); |
---|
61 | end |
---|
62 | |
---|
63 | D5 = d*eye(5); |
---|
64 | %B5 = zeros(5,5); |
---|
65 | RMATi = zeros(6,6); |
---|
66 | theta5 = [0 0 0 0 dCT]'; |
---|
67 | |
---|
68 | |
---|
69 | |
---|
70 | %Fist iteration |
---|
71 | RMATi= Ri*ones(1,6) + [D5 zeros(5,1); zeros(1,6)]; |
---|
72 | RMATf = linepass(RING,RMATi); |
---|
73 | Rf = RMATf(:,6); |
---|
74 | % compute the transverse part of the Jacobian |
---|
75 | J5 = [RMATf([1:4,6],1:5)-RMATf([1:4,6],6)*ones(1,5)]/d; |
---|
76 | |
---|
77 | % Replace matrix inversion with \ |
---|
78 | %B5 = inv(diag([1 1 1 1 0]) - J5); |
---|
79 | % Ri_next = Ri + [B5* (Rf([1:4,6])-[Ri(1:4);0]-theta5) ; 0]; |
---|
80 | Ri_next = Ri + [(diag([1 1 1 1 0]) - J5)\(Rf([1:4,6])-[Ri(1:4);0]-theta5) ; 0]; |
---|
81 | change = norm(Ri_next - Ri); |
---|
82 | Ri = Ri_next; |
---|
83 | |
---|
84 | itercount = 1; |
---|
85 | |
---|
86 | |
---|
87 | while (change>eps) & (itercount < max_iterations) |
---|
88 | |
---|
89 | RMATi= Ri*ones(1,6) + [D5 zeros(5,1); zeros(1,6)]; |
---|
90 | RMATf = linepass(RING,RMATi,'reuse'); |
---|
91 | |
---|
92 | Rf = RMATf(:,6); |
---|
93 | % compute the transverse part of the Jacobian |
---|
94 | J5 = [RMATf([1:4,6],1:5)-RMATf([1:4,6],6)*ones(1,5)]/d; |
---|
95 | % Replace matrix inversion with \ |
---|
96 | %B5 = inv(diag([1 1 1 1 0]) - J5); |
---|
97 | %Ri_next = Ri + [B5*(Rf([1:4,6])-[Ri(1:4);0]-theta5); 0]; |
---|
98 | Ri_next = Ri + [(diag([1 1 1 1 0]) - J5)\(Rf([1:4,6])-[Ri(1:4);0]-theta5); 0]; |
---|
99 | change = norm(Ri_next - Ri); |
---|
100 | Ri = Ri_next; |
---|
101 | itercount = itercount+1; |
---|
102 | |
---|
103 | end; |
---|
104 | |
---|
105 | |
---|
106 | if(nargin<3) | (varargin{1}==(length(RING)+1)) |
---|
107 | % return only the fixed point at the entrance of RING{1} |
---|
108 | orbit=Ri(1:5,1); |
---|
109 | else % 3-rd input argument - vector of reference points along the Ring |
---|
110 | % is supplied - return orbit |
---|
111 | orb6 = linepass(RING,Ri,varargin{1},'reuse'); |
---|
112 | orbit = orb6(1:5,:); |
---|
113 | end |
---|
114 | |
---|
115 | if nargout==2 |
---|
116 | varargout{1}=Ri; |
---|
117 | end |
---|