findorbit4

FINDORBIT4 finds closed orbit in the 4-d transverse phase
 space 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 the CONSTANT MOMENTUM constraint, dP2 = dP1 = dP and
    there is NO constraint on the 6-th coordinate CT

 IMPORTANT!!! FINDORBIT4 imposes a constraint on dP and relaxes
    the constraint on the revolution frequency. A physical storage
    ring does exactly the opposite: the momentum deviation of a
    particle on the closed orbit settles at the value
    such that the revolution is synchronous with the RF cavity

                 HarmNumber*Frev = Frf

    To impose this artificial constraint in FINDORBIT4
    PassMethod used for any elemen SHOULD NOT
    1. change the longitudinal momentum dP (cavities , magnets with radiation)
    2. have any time dependence (localized impedance, fast kickers etc)

 FINDORBIT4(RING,dP) is 4x1 vector - fixed point at the
    entrance of the 1-st element of the RING (x,px,y,py)

 FINDORBIT4(RING,dP,REFPTS) is 4-by-Length(REFPTS)
     array of column vectors - fixed points (x,px,y,py)
     at the entrance of each element indexed  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

 FINDORBIT4(RING,dP,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] = FINDORBIT4( ... )
     The optional second return parameter is
     a 6-by-1 vector of initial conditions
     on the closed orbit at the entrance to the RING.

 See also FINDSYNCORBIT, FINDORBIT6.

0001 function [orbit,fixedpoint] = findorbit4(RING,dP,varargin)
0002 %FINDORBIT4 finds closed orbit in the 4-d transverse phase
0003 % space by numerically solving  for a fixed point of the one turn
0004 % map M calculated with LINEPASS
0005 %
0006 %         (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
0007 %
0008 %    under the CONSTANT MOMENTUM constraint, dP2 = dP1 = dP and
0009 %    there is NO constraint on the 6-th coordinate CT
0010 %
0011 % IMPORTANT!!! FINDORBIT4 imposes a constraint on dP and relaxes
0012 %    the constraint on the revolution frequency. A physical storage
0013 %    ring does exactly the opposite: the momentum deviation of a
0014 %    particle on the closed orbit settles at the value
0015 %    such that the revolution is synchronous with the RF cavity
0016 %
0017 %                 HarmNumber*Frev = Frf
0018 %
0019 %    To impose this artificial constraint in FINDORBIT4
0020 %    PassMethod used for any elemen SHOULD NOT
0021 %    1. change the longitudinal momentum dP (cavities , magnets with radiation)
0022 %    2. have any time dependence (localized impedance, fast kickers etc)
0023 %
0024 % FINDORBIT4(RING,dP) is 4x1 vector - fixed point at the
0025 %    entrance of the 1-st element of the RING (x,px,y,py)
0026 %
0027 % FINDORBIT4(RING,dP,REFPTS) is 4-by-Length(REFPTS)
0028 %     array of column vectors - fixed points (x,px,y,py)
0029 %     at the entrance of each element indexed  REFPTS array.
0030 %     REFPTS is an array of increasing indexes that  select elements
0031 %     from the range 1 to length(RING)+1.
0032 %     See further explanation of REFPTS in the 'help' for FINDSPOS
0033 %
0034 % FINDORBIT4(RING,dP,REFPTS,GUESS) - same as above but the search
0035 %     for the fixed point starts at the initial condition GUESS
0036 %     Otherwise the search starts from [0 0 0 0 0 0]'.
0037 %     GUESS must be a 6-by-1 vector;
0038 %
0039 % [ORBIT, FIXEDPOINT] = FINDORBIT4( ... )
0040 %     The optional second return parameter is
0041 %     a 6-by-1 vector of initial conditions
0042 %     on the closed orbit at the entrance to the RING.
0043 %
0044 % See also FINDSYNCORBIT, FINDORBIT6.
0045 
0046 if ~iscell(RING)
0047     error('First argument must be a cell array');
0048 end
0049 
0050 d = 1e-6;       % step size for numerical differentiation
0051 dps = 1e-12;    % convergence threshold
0052 %dps=eps;       % convergence threshold
0053 max_iterations = 20;
0054 
0055 if nargin >= 4    % Check if guess argument was supplied
0056     if isnumeric(varargin{2}) && all(eq(size(varargin{2}),[6,1]))
0057         Ri=varargin{2};
0058     else
0059         error('The last argument GUESS must be a 6-by-1 vector');
0060     end
0061 else
0062     Ri = zeros(6,1);
0063 end
0064 % Set the momentum component of Ri to the specified dP
0065 Ri(5) = dP;
0066 D = [d*eye(4) zeros(4,1);zeros(2,5)];
0067 %D = [0.5*d*eye(4) -0.5*d*eye(4) zeros(4,1);zeros(2,9)];
0068 
0069 args={};
0070 change=Inf;
0071 itercount = 0;
0072 while (change > dps) && (itercount < max_iterations)
0073     RMATi = Ri(:,ones(1,5)) + D;
0074     %RMATi = Ri(:,ones(1,9)) + D;
0075     RMATf = linepass(RING,RMATi,args{:});
0076     Rf = RMATf(:,end);
0077     % compute the transverse part of the Jacobian
0078     J4 = (RMATf(1:4,1:4)-RMATf(1:4,5*ones(1,4)))/d;
0079     %J4 = (RMATf(1:4,1:4)-RMATf(1:4,5:8))/d;
0080     Ri_next = Ri +  [(eye(4) - J4)\(Rf(1:4)-Ri(1:4)); 0; 0];
0081     change = norm(Ri_next - Ri);
0082     Ri = Ri_next;
0083     itercount = itercount+1;
0084     args={'KeepLattice'};
0085 end
0086 
0087 if itercount == max_iterations
0088     warning('Maximum number of iterations reached. Possible non-convergence')
0089 end
0090 
0091 if (nargin<3) || (isscalar(varargin{1}) && (varargin{1}==(length(RING)+1)))
0092     % return only the fixed point at the entrance of RING{1}
0093     orbit=Ri(1:4,1);
0094 else    % 3-rd input argument - vector of reference points along the RING
0095     % is supplied - return orbit
0096     orb6 = linepass(RING,Ri,varargin{1},'KeepLattice');
0097     orbit = orb6(1:4,:);
0098 end
0099 
0100 if nargout >= 2
0101     fixedpoint=Ri;
0102 end