findm66

FINDM66 numerically finds the 6x6 transfer matrix of an accelerator lattice
  by differentiation of LINEPASS near the closed orbit
  FINDM66 uses FINDORBIT6 to search for the closed orbit in 6-d
  In order for this to work the ring MUST have a CAVITY element

 M66 = FINDM66(RING)finds full one-turn 6-by-6
    matrix at the entrance of the first element

 [M66,T] = FINDM66(RING,REFPTS) in addition to M finds
    6-by-6 transfer matrixes  between entrances of
    the first element and each element indexed by REFPTS.
    T is 6-by-6-by-length(REFPTS) 3 dimentional 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

    Note:
    When REFPTS= [ 1 2 .. ] the fist point is the entrance of the first element
    and T(:,:,1) - identity matrix
    When REFPTS= [  .. length(RING)+1] the last point is the exit of the last element
    and the entrance of the first element after 1 turn: T(:,:, ) = M

 [M66,T] = FINDM66(RING,REFPTS,ORBITIN) - Does not search for closed orbit.
   This syntax is useful to avoid recomputing the closed orbit if it is
   already known.

 [M66,T,orbit] = FINDM66(RING, REFPTS) in addition returns the closed orbit
    found in the process of lenearization

 See also FINDM44, FINDORBIT6

0001 function [M66, varargout] = findm66(LATTICE, varargin)
0002 %FINDM66 numerically finds the 6x6 transfer matrix of an accelerator lattice
0003 %  by differentiation of LINEPASS near the closed orbit
0004 %  FINDM66 uses FINDORBIT6 to search for the closed orbit in 6-d
0005 %  In order for this to work the ring MUST have a CAVITY element
0006 %
0007 % M66 = FINDM66(RING)finds full one-turn 6-by-6
0008 %    matrix at the entrance of the first element
0009 %
0010 % [M66,T] = FINDM66(RING,REFPTS) in addition to M finds
0011 %    6-by-6 transfer matrixes  between entrances of
0012 %    the first element and each element indexed by REFPTS.
0013 %    T is 6-by-6-by-length(REFPTS) 3 dimentional array.
0014 %
0015 %    REFPTS is an array of increasing indexes that  select elements
0016 %    from the range 1 to length(RING)+1.
0017 %    See further explanation of REFPTS in the 'help' for FINDSPOS
0018 %
0019 %    Note:
0020 %    When REFPTS= [ 1 2 .. ] the fist point is the entrance of the first element
0021 %    and T(:,:,1) - identity matrix
0022 %    When REFPTS= [  .. length(RING)+1] the last point is the exit of the last element
0023 %    and the entrance of the first element after 1 turn: T(:,:, ) = M
0024 %
0025 % [M66,T] = FINDM66(RING,REFPTS,ORBITIN) - Does not search for closed orbit.
0026 %   This syntax is useful to avoid recomputing the closed orbit if it is
0027 %   already known.
0028 %
0029 % [M66,T,orbit] = FINDM66(RING, REFPTS) in addition returns the closed orbit
0030 %    found in the process of lenearization
0031 %
0032 % See also FINDM44, FINDORBIT6
0033 
0034 if ~iscell(LATTICE)
0035     error('First argument must be a cell array');
0036 end
0037 
0038 NE = length(LATTICE);
0039 
0040 if nargin >= 3 && ~isempty(varargin{2}) % FINDM66(RING,REFPTS,ORBITIN)
0041     R0 = varargin{2};
0042 else
0043     R0 = findorbit6(LATTICE);
0044 end
0045 if nargin >= 2  % FINDM66(RING,REFPTS)
0046     if islogical(varargin{1})
0047         REFPTS=varargin{1};
0048         REFPTS(end+1:NE+1)=false;
0049     elseif isnumeric(varargin{1})
0050         REFPTS=setelems(false(1,NE+1),varargin{1});
0051     else
0052         error('REFPTS must be numeric or logical');
0053     end
0054 else
0055     REFPTS=false(1,NE+1);
0056 end
0057 refs=setelems(REFPTS,NE+1);
0058 reqs=REFPTS(refs);
0059 
0060 % Determine step size to use for numerical differentiation
0061 global NUMDIFPARAMS
0062 % Transverse
0063 if isfield(NUMDIFPARAMS,'XYStep')
0064     dt = NUMDIFPARAMS.XYStep';
0065 else
0066     % optimal differentiation step - Numerical Recipes
0067     dt =  6.055454452393343e-006;
0068 end
0069 % Longitudinal
0070 if isfield(NUMDIFPARAMS,'DPStep')
0071     dl = NUMDIFPARAMS.DPStep';
0072 else
0073     % optimal differentiation step - Numerical Recipes
0074     dl =  6.055454452393343e-006;
0075 end
0076 
0077 % Build a diagonal matrix of initial conditions
0078 D6 = [0.5*dt*eye(4),zeros(4,2);zeros(2,4),0.5*dl*eye(2)];
0079 % Add to the orbit_in. First 12 columns for derivative
0080 % 13-th column is for closed orbit
0081 RIN = R0(:,ones(1,13)) + [D6 -D6 zeros(6,1)];
0082 ROUT = linepass(LATTICE,RIN,refs);
0083 TMAT3 = reshape(ROUT,6,13,[]);
0084 M66 = [(TMAT3(:,1:4,end)-TMAT3(:,7:10,end))./dt,...
0085     (TMAT3(:,5:6,end)-TMAT3(:,11:12,end))./dl];
0086 
0087 if nargout >= 2 % Calculate matrixes at all REFPTS.
0088     varargout{1} = [(TMAT3(:,1:4,reqs)-TMAT3(:,7:10,reqs))./dt,...
0089         (TMAT3(:,5:6,reqs)-TMAT3(:,11:12,reqs))./dl];
0090     % Return closed orbit if requested
0091     if nargout >= 3
0092         varargout{2}=squeeze(TMAT3(:,13,reqs));
0093     end
0094 end
0095 
0096     function mask=setelems(mask,idx)
0097         mask(idx)=true;
0098     end
0099 
0100 end