FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity
and momentum deviation dP (first 5 components of the phase space vector)
by numerically solving for a fixed point
of the one turn map M calculated with LINEPASS
(X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
under constraints CT2 - CT1 = dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where
Frev0 = Frf0/HarmNumber is the design revolution frequency
Frev = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency
IMPORTANT!!! FINDSYNCORBIT imposes a constraint (CT2 - CT1) and
dP2 = dP1 but no constraint on the value of dP1, dP2
The algorithm assumes time-independent fixed-momentum ring
to reduce the dimensionality of the problem.
To impose this artificial constraint in FINDSYNCORBIT
PassMethod used for any element SHOULD NOT
1. change the longitudinal momentum dP (cavities , magnets with radiation)
2. have any time dependence (localized impedance, fast kickers etc).
FINDSYNCORBIT(RING,dCT) is 5x1 vector - fixed point at the
entrance of the 1-st element of the RING (x,px,y,py,dP)
FINDSYNCORBIT(RING,dCT,REFPTS) is 5-by-Length(REFPTS)
array of column vectors - fixed points (x,px,y,py,dP)
at the entrance of each element indexed in REFPTS array.
REFPTS is an array of increasing indexes that select elements
from the range 1 to length(RING)+1.
See further explanation of REFPTS in the 'help' for FINDSPOS
FINDSYNCORBIT(RING,dCT,REFPTS,GUESS) - same as above but the search
for the fixed point starts at the initial condition GUESS
Otherwise the search starts from [0 0 0 0 0 0]'.
GUESS must be a 6-by-1 vector;
[ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... )
The optional second return parameter is
a 6-by-1 vector of initial conditions
on the close orbit at the entrance to the RING.
See also FINDORBIT, FINDORBIT4, FINDORBIT6.0001 function [orbit, varargout] = findsyncorbit(RING,dCT,varargin); 0002 %FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity 0003 % and momentum deviation dP (first 5 components of the phase space vector) 0004 % by numerically solving for a fixed point 0005 % of the one turn map M calculated with LINEPASS 0006 % 0007 % (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1) 0008 % 0009 % under constraints CT2 - CT1 = dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where 0010 % Frev0 = Frf0/HarmNumber is the design revolution frequency 0011 % Frev = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency 0012 % 0013 % IMPORTANT!!! FINDSYNCORBIT imposes a constraint (CT2 - CT1) and 0014 % dP2 = dP1 but no constraint on the value of dP1, dP2 0015 % The algorithm assumes time-independent fixed-momentum ring 0016 % to reduce the dimensionality of the problem. 0017 % To impose this artificial constraint in FINDSYNCORBIT 0018 % PassMethod used for any element SHOULD NOT 0019 % 1. change the longitudinal momentum dP (cavities , magnets with radiation) 0020 % 2. have any time dependence (localized impedance, fast kickers etc). 0021 % 0022 % 0023 % FINDSYNCORBIT(RING,dCT) is 5x1 vector - fixed point at the 0024 % entrance of the 1-st element of the RING (x,px,y,py,dP) 0025 % 0026 % FINDSYNCORBIT(RING,dCT,REFPTS) is 5-by-Length(REFPTS) 0027 % array of column vectors - fixed points (x,px,y,py,dP) 0028 % at the entrance of each element indexed in REFPTS array. 0029 % REFPTS is an array of increasing indexes that select elements 0030 % from the range 1 to length(RING)+1. 0031 % See further explanation of REFPTS in the 'help' for FINDSPOS 0032 % 0033 % FINDSYNCORBIT(RING,dCT,REFPTS,GUESS) - same as above but the search 0034 % for the fixed point starts at the initial condition GUESS 0035 % Otherwise the search starts from [0 0 0 0 0 0]'. 0036 % GUESS must be a 6-by-1 vector; 0037 % 0038 % [ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... ) 0039 % The optional second return parameter is 0040 % a 6-by-1 vector of initial conditions 0041 % on the close orbit at the entrance to the RING. 0042 % 0043 % See also FINDORBIT, FINDORBIT4, FINDORBIT6. 0044 % 0045 if ~iscell(RING) 0046 error('First argument must be a cell array'); 0047 end 0048 0049 d = 1e-9; % step size for numerical differentiation 0050 max_iterations = 20; 0051 J = zeros(5); 0052 0053 if nargin==4 0054 if isnumeric(varargin{2}) & all(eq(size(varargin{2}),[6,1])) 0055 Ri=varargin{2}; 0056 else 0057 error('The last argument GUESS must be a 6-by-1 vector'); 0058 end 0059 else 0060 Ri = zeros(6,1); 0061 end 0062 0063 D5 = d*eye(5); 0064 %B5 = zeros(5,5); 0065 RMATi = zeros(6,6); 0066 theta5 = [0 0 0 0 dCT]'; 0067 0068 0069 0070 %Fist iteration 0071 RMATi= Ri*ones(1,6) + [D5 zeros(5,1); zeros(1,6)]; 0072 RMATf = linepass(RING,RMATi); 0073 Rf = RMATf(:,6); 0074 % compute the transverse part of the Jacobian 0075 J5 = [RMATf([1:4,6],1:5)-RMATf([1:4,6],6)*ones(1,5)]/d; 0076 0077 % Replace matrix inversion with \ 0078 %B5 = inv(diag([1 1 1 1 0]) - J5); 0079 % Ri_next = Ri + [B5* (Rf([1:4,6])-[Ri(1:4);0]-theta5) ; 0]; 0080 Ri_next = Ri + [(diag([1 1 1 1 0]) - J5)\(Rf([1:4,6])-[Ri(1:4);0]-theta5) ; 0]; 0081 change = norm(Ri_next - Ri); 0082 Ri = Ri_next; 0083 0084 itercount = 1; 0085 0086 0087 while (change>eps) & (itercount < max_iterations) 0088 0089 RMATi= Ri*ones(1,6) + [D5 zeros(5,1); zeros(1,6)]; 0090 RMATf = linepass(RING,RMATi,'reuse'); 0091 0092 Rf = RMATf(:,6); 0093 % compute the transverse part of the Jacobian 0094 J5 = [RMATf([1:4,6],1:5)-RMATf([1:4,6],6)*ones(1,5)]/d; 0095 % Replace matrix inversion with \ 0096 %B5 = inv(diag([1 1 1 1 0]) - J5); 0097 %Ri_next = Ri + [B5*(Rf([1:4,6])-[Ri(1:4);0]-theta5); 0]; 0098 Ri_next = Ri + [(diag([1 1 1 1 0]) - J5)\(Rf([1:4,6])-[Ri(1:4);0]-theta5); 0]; 0099 change = norm(Ri_next - Ri); 0100 Ri = Ri_next; 0101 itercount = itercount+1; 0102 0103 end; 0104 0105 0106 if(nargin<3) | (varargin{1}==(length(RING)+1)) 0107 % return only the fixed point at the entrance of RING{1} 0108 orbit=Ri(1:5,1); 0109 else % 3-rd input argument - vector of reference points along the Ring 0110 % is supplied - return orbit 0111 orb6 = linepass(RING,Ri,varargin{1},'reuse'); 0112 orbit = orb6(1:5,:); 0113 end 0114 0115 if nargout==2 0116 varargout{1}=Ri; 0117 end