findelemm44

FINDELEMM44 numerically finds the 4x4 transfer matrix of an element
  FINDELEM66(ELEM, METHODNAME, ORBITIN)
     ELEM          - the element data structure
     METHODNAME    - name of the pass-method function
                   (default:  ELEM.PassMethod)
     ORBITIN       - 6-by-1 phase space coordinates at the entrance
                   (default: zeros(6,1))
                   The transverse matrix is momentum-dependent,
                   the 5-th component of ORBITIN is used as the DP value

 See also FINDELEMM66

0001 function M44 = findelemm44(ELEM, MethodName, R0)
0002 %FINDELEMM44 numerically finds the 4x4 transfer matrix of an element
0003 %  FINDELEM66(ELEM, METHODNAME, ORBITIN)
0004 %     ELEM          - the element data structure
0005 %     METHODNAME    - name of the pass-method function
0006 %                   (default:  ELEM.PassMethod)
0007 %     ORBITIN       - 6-by-1 phase space coordinates at the entrance
0008 %                   (default: zeros(6,1))
0009 %                   The transverse matrix is momentum-dependent,
0010 %                   the 5-th component of ORBITIN is used as the DP value
0011 %
0012 % See also FINDELEMM66
0013 
0014 if (nargin < 3) || isempty(R0),  R0 = zeros(6,1); end
0015 if (nargin < 2) || isempty(MethodName),  MethodName=ELEM.PassMethod; end
0016 
0017 % Determine step size to use for numerical differentiation
0018 global NUMDIFPARAMS
0019 % Transverse
0020 if isfield(NUMDIFPARAMS,'XYStep')
0021     dt = NUMDIFPARAMS.XYStep';
0022 else
0023     % optimal differentiation step - Numerical Recipes
0024     dt =  6.055454452393343e-006;
0025 end
0026 
0027 % Build a diagonal matrix of initial conditions
0028 D4 = [0.5*dt*eye(4);zeros(2,4)];
0029 % Add to the orbit_in
0030 RIN = R0(:,ones(1,8)) + [D4 -D4];
0031 % Propagate through the element
0032 ROUT = feval(MethodName,ELEM,RIN);
0033 % Calculate numerical derivative
0034 M44 = ((ROUT(1:4,1:4)-ROUT(1:4,5:8))./dt);