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