locoresponsematrix

LOCORESPONSEMATRIX - Calculate the BPM response matrix and dispersion function
 M = LOCORESPONCEMATRIX(RINGData, BPMData, CMData)

 Accelerator Toolbox implementation of generic LOCO function

 RINGData - must have fields 'Lattice', 'CavityFrequency', 'CavityHarmNumber'
            RINGData.Lattice - AT lattice cell arrary
            RINGData.CavityFrequency  [Hz]
            RINGData.CavityHarmNumber [Hz]

 CMData -   must have fields: 'HCMIndex', 'VCMIndex', 'HCMKicks', 'VCMKicks', 'HCMCoupling', 'VCMCoupling'
            CMData.HCMIndex    - indexes in the AT lattice of elements used as horizontal correctors
            CMData.VCMIndex    - indexes in the AT lattice of elements used as vertical correctors
                                 Elements used as correctors in both planes should be included in both lists
            CMData.HCMKicks    - kick  size [radians] of horizontal correctors in the horizontal plane
            CMData.VCMKicks    - kick  size [radians] of vertical correctors in the vertical plane
            CMData.HCMCoupling - corrector coupling coefficient into another plane:
                                 0.01 coupling means that for 1e-3 kick in the horizontal direction there
                                 is a 1e-5 rad kick in the vertical direction
            CMData.VCMCoupling - corrector coupling coefficient into another plane:
                                 0.01 coupling means that for 1e-3 kick in the vertical direction there
                                 is a 1e-5 rad kick in the horizontal direction

 BPMData -  must have field 'BPMIndex'
            CMData.BPMIndex - indexes of all BPMs or observation points in the AT lattice
                              All BPS and observation points (single plane too)
                              are included in CMData.BPMIndex.

 Return value: a matrix with number of rows equal to 2*length(CMData.BPMIndex) and the number of columns
               equal length(CMData.HCMIndex)+length(CMData.VCMIndex)

 Additional string flags (in any order)

 LOCORESPONSEMATRIX(...,ClosedOrbitType,...)
       ClosedOrbitType is 'fixedmomentum',  'fixedpathlength' (default)

 LOCORESPONCEMATRIX(..., 'linear') calculates M using linear approximation  !!! including the dispersion terms

 LOCORESPONCEMATRIX(..., 'RF', DeltaRF) - 'RF' switch must be followed by the value of DeltaRF [Hz]

 LOCORESPONCEMATRIX(..., 'ResponseMatrixMeasurement', 'oneway') - 'oneway' switch is used
                        when the response matrix was measured only kicking the i-th corrector
                        to +KicksCoupled(i) one way, default: ResponseMatrixMeasurement = 'bidirectional'

 LOCORESPONCEMATRIX(..., 'DispersionMeasurement', 'oneway') - 'oneway' switch is used
                        when the dispersion was measured only by varying the
                        RF frequency in one direction, default: DispersionMeasurement = 'bidirectional'

 Or a Flags structure can be an input argument:
 LOCORESPONCEMATRIX(..., Flags)
    Flags.ResponseMatrixMeasurement = 'oneway' or {'bi-directional'}
    Flags.DispersionMeasurement     = 'oneway' or {'bi-directional'}
    Flags.ResponseMatrixCalculator  = {'linear'} or 'full'
    Flags.ClosedOrbitType           = 'fixedmomentum' or {'fixedpathlength'}

  NOTE
  1. Flag names are not case sensitive
  2. JR - 12/7/05 added vectorized linear rm calculation (NewVectorizedMethod = 1)

0001 function RM = locoresponsematrix(RINGData, BPMData, CMData, varargin)
0002 %LOCORESPONSEMATRIX - Calculate the BPM response matrix and dispersion function
0003 % M = LOCORESPONCEMATRIX(RINGData, BPMData, CMData)
0004 %
0005 % Accelerator Toolbox implementation of generic LOCO function
0006 %
0007 % RINGData - must have fields 'Lattice', 'CavityFrequency', 'CavityHarmNumber'
0008 %            RINGData.Lattice - AT lattice cell arrary
0009 %            RINGData.CavityFrequency  [Hz]
0010 %            RINGData.CavityHarmNumber [Hz]
0011 %
0012 % CMData -   must have fields: 'HCMIndex', 'VCMIndex', 'HCMKicks', 'VCMKicks', 'HCMCoupling', 'VCMCoupling'
0013 %            CMData.HCMIndex    - indexes in the AT lattice of elements used as horizontal correctors
0014 %            CMData.VCMIndex    - indexes in the AT lattice of elements used as vertical correctors
0015 %                                 Elements used as correctors in both planes should be included in both lists
0016 %            CMData.HCMKicks    - kick  size [radians] of horizontal correctors in the horizontal plane
0017 %            CMData.VCMKicks    - kick  size [radians] of vertical correctors in the vertical plane
0018 %            CMData.HCMCoupling - corrector coupling coefficient into another plane:
0019 %                                 0.01 coupling means that for 1e-3 kick in the horizontal direction there
0020 %                                 is a 1e-5 rad kick in the vertical direction
0021 %            CMData.VCMCoupling - corrector coupling coefficient into another plane:
0022 %                                 0.01 coupling means that for 1e-3 kick in the vertical direction there
0023 %                                 is a 1e-5 rad kick in the horizontal direction
0024 %
0025 % BPMData -  must have field 'BPMIndex'
0026 %            CMData.BPMIndex - indexes of all BPMs or observation points in the AT lattice
0027 %                              All BPS and observation points (single plane too)
0028 %                              are included in CMData.BPMIndex.
0029 %
0030 % Return value: a matrix with number of rows equal to 2*length(CMData.BPMIndex) and the number of columns
0031 %               equal length(CMData.HCMIndex)+length(CMData.VCMIndex)
0032 %
0033 % Additional string flags (in any order)
0034 %
0035 % LOCORESPONSEMATRIX(...,ClosedOrbitType,...)
0036 %       ClosedOrbitType is 'fixedmomentum',  'fixedpathlength' (default)
0037 %
0038 % LOCORESPONCEMATRIX(..., 'linear') calculates M using linear approximation  !!! including the dispersion terms
0039 %
0040 % LOCORESPONCEMATRIX(..., 'RF', DeltaRF) - 'RF' switch must be followed by the value of DeltaRF [Hz]
0041 %
0042 % LOCORESPONCEMATRIX(..., 'ResponseMatrixMeasurement', 'oneway') - 'oneway' switch is used
0043 %                        when the response matrix was measured only kicking the i-th corrector
0044 %                        to +KicksCoupled(i) one way, default: ResponseMatrixMeasurement = 'bidirectional'
0045 %
0046 % LOCORESPONCEMATRIX(..., 'DispersionMeasurement', 'oneway') - 'oneway' switch is used
0047 %                        when the dispersion was measured only by varying the
0048 %                        RF frequency in one direction, default: DispersionMeasurement = 'bidirectional'
0049 %
0050 % Or a Flags structure can be an input argument:
0051 % LOCORESPONCEMATRIX(..., Flags)
0052 %    Flags.ResponseMatrixMeasurement = 'oneway' or {'bi-directional'}
0053 %    Flags.DispersionMeasurement     = 'oneway' or {'bi-directional'}
0054 %    Flags.ResponseMatrixCalculator  = {'linear'} or 'full'
0055 %    Flags.ClosedOrbitType           = 'fixedmomentum' or {'fixedpathlength'}
0056 %
0057 %  NOTE
0058 %  1. Flag names are not case sensitive
0059 %  2. JR - 12/7/05 added vectorized linear rm calculation (NewVectorizedMethod = 1)
0060 
0061 NewVectorizedMethod = 1;
0062 
0063 C = 2.99792458e8;
0064 
0065 % Defaults
0066 ResponseMatrixMeasurement = 'bidirectional';
0067 DispersionMeasurement     = 'bidirectional';
0068 ResponseMatrixCalculator  = 'linear';
0069 ClosedOrbitType           = 'fixedpathlength';
0070 MachineType               = 'StorageRing';
0071 
0072 RFFLAG = 0;
0073 DeltaRF = [];
0074 
0075 N = nargin-3;
0076 i = 0;
0077 while i < N
0078     i = i + 1;
0079     if isstruct(varargin{i})
0080         Flags = varargin{i};
0081         if isfield(Flags,'ResponseMatrixCalculator')
0082             ResponseMatrixCalculator = Flags.ResponseMatrixCalculator;
0083         end
0084         if isfield(Flags,'ClosedOrbitType')
0085             ClosedOrbitType = Flags.ClosedOrbitType;
0086         end
0087         if isfield(Flags,'ResponseMatrixMeasurement')
0088             ResponseMatrixMeasurement = Flags.ResponseMatrixMeasurement;
0089         end
0090         if isfield(Flags,'DispersionMeasurement')
0091             DispersionMeasurement = Flags.DispersionMeasurement;
0092         end
0093         if isfield(Flags,'MachineType')
0094             MachineType = Flags.MachineType;
0095         end
0096     elseif ischar(varargin{i})
0097         switch lower(varargin{i})
0098             case 'linear'
0099                 ResponseMatrixCalculator = 'linear';
0100             case 'full'
0101                 ResponseMatrixCalculator = 'full';
0102             case 'fixedmomentum'
0103                 ClosedOrbitType = 'fixedmomentum';
0104             case 'fixedpathlength'
0105                 ClosedOrbitType = 'fixedpathlength';
0106             case 'rf'
0107                 if (i+4<=nargin) && isnumeric(varargin{i+1})
0108                     RFFLAG = 1;
0109                     DeltaRF = varargin{i+1};
0110                     i = i + 1;
0111                 else
0112                     error('''RF'' flag must be followed by a numeric value of delta RF [Hz]');
0113                 end
0114             case 'dispersionmeasurement'
0115                 if (i+4<=nargin) && ischar(varargin{i+1})
0116                     DispersionMeasurement = varargin{i+1};
0117                     i = i + 1;
0118                 else
0119                     error('''DispersionMeasurement'' flag must be followed by ''oneway'' or ''bidirectional''');
0120                 end
0121             case 'responsematrixmeasurement'
0122                 if (i+4<=nargin) && ischar(varargin{i+1})
0123                     ResponseMatrixMeasurement = varargin{i+1};
0124                     i = i + 1;
0125                 else
0126                     error('''ResponseMatrixMeasurement'' flag must be followed by ''oneway'' or ''bidirectional''');
0127                 end
0128             case {'storagering','booster','boosterring'}
0129                 MachineType = 'StorageRing';
0130             case {'transport','linac'}
0131                 MachineType = 'Transport';
0132             otherwise
0133                 warning('Unknown switch ignored.');
0134         end
0135     else
0136         warning('Unknown switch ignored.');
0137     end
0138 end
0139 
0140 
0141 % Input checks
0142 if ~strcmpi(ResponseMatrixMeasurement, 'bidirectional') && ~strcmpi(ResponseMatrixMeasurement, 'oneway')
0143     error('Unknown ResponseMatrixMeasurement type');
0144 end
0145 if ~strcmpi(DispersionMeasurement, 'bidirectional') && ~strcmpi(DispersionMeasurement, 'oneway')
0146     error('Unknown DispersionMeasurement type');
0147 end
0148 if ~strcmpi(ResponseMatrixCalculator, 'linear') && ~strcmpi(ResponseMatrixCalculator, 'full')
0149     error('Unknown ResponseMatrixCalculator method');
0150 end
0151 if ~strcmpi(ClosedOrbitType, 'fixedpathlength') && ~strcmpi(ClosedOrbitType, 'fixedmomentum')
0152     error('Unknown ClosedOrbitType method');
0153 end
0154 if ~strcmpi(MachineType, 'StorageRing') && ~strcmpi(MachineType, 'Transport')
0155     error('Unknown MachineType');
0156 end
0157 
0158 
0159 % Initialize
0160 ATRING = RINGData.Lattice;
0161 NHC = length(CMData.HCMIndex);
0162 NVC = length(CMData.VCMIndex);
0163 NBPM = length(BPMData.BPMIndex);
0164 
0165 
0166 if strcmpi(MachineType, 'StorageRing')
0167 
0168     if RFFLAG
0169         % Add an extra column for the orbit responce to the RF frequency change
0170         RM = zeros(2*NBPM,NHC+NVC+1);
0171     else
0172         RM = zeros(2*NBPM,NHC+NVC);
0173     end
0174 
0175     NE = length(ATRING);
0176     if strcmpi(ResponseMatrixCalculator, 'Linear')
0177         % Calculate linear optics and chromatic finctions for the model
0178         [M44,T,ClosedOrbit] = findm44(ATRING,0,1:NE+1);
0179         DP = 0.00001;
0180         ClosedOrbitDP = findorbit4(ATRING,DP,1:NE+1);
0181         Dispersion = (ClosedOrbitDP-ClosedOrbit)/DP;
0182         L0 = findspos(ATRING,NE+1);
0183 
0184         %%X1(6)/(DP*L0); % is not the same as mcf(ATRING)???  G. Portmann
0185         %X1 = ringpass(ATRING,[ClosedOrbitDP(:,1);DP;0]);
0186         %MCF = X1(6)/(DP*L0);
0187         MCF = mcf(ATRING);
0188 
0189         % Transfer matrixes through individual correctors
0190         M44HCOR = cell(1,NHC);
0191         M44VCOR = cell(1,NVC);
0192         for i=1:NHC
0193             M44HCOR{i}=findelemm44(ATRING{CMData.HCMIndex(i)},ATRING{CMData.HCMIndex(i)}.PassMethod,[ClosedOrbit(:,CMData.HCMIndex(i));0;0]);
0194         end
0195         for i=1:NVC
0196             match = find(CMData.VCMIndex(i)==CMData.HCMIndex);
0197             if match
0198                 M44VCOR{i}=M44HCOR{match};
0199             else
0200                 M44VCOR{i}=findelemm44(ATRING{CMData.VCMIndex(i)},ATRING{CMData.VCMIndex(i)}.PassMethod,[ClosedOrbit(:,CMData.VCMIndex(i));0;0]);
0201             end
0202         end
0203 
0204 
0205         % Assemble arrays of corrector kicks including coupling
0206         HCORTheta = zeros(4,NHC);
0207         VCORTheta = zeros(4,NVC);
0208 
0209         HCORTheta(2,:) = CMData.HCMKicks(:)';
0210         HCORTheta(4,:) = CMData.HCMCoupling(:)' .* CMData.HCMKicks(:)';
0211         VCORTheta(2,:) = CMData.VCMCoupling(:)' .* CMData.VCMKicks(:)';
0212         VCORTheta(4,:) = CMData.VCMKicks(:)';
0213 
0214 
0215         % Calculate closed orbit at the exit of each corrector magnet WITH applied kick
0216         for i=1:NHC
0217             CI = CMData.HCMIndex(i);
0218             InverseT = inv(T(:,:,CI));
0219             OrbitEntrance = (inv(eye(4)-T(:,:,CI)*M44*InverseT)*...
0220                 T(:,:,CI)*M44*InverseT*(eye(4)+inv(M44HCOR{i}))*HCORTheta(:,i)/2);
0221 
0222             OrbitExit = HCORTheta(:,i)/2+M44HCOR{i}*(OrbitEntrance+HCORTheta(:,i)/2);
0223 
0224 
0225             R0 = inv(T(:,:,CI+1))*OrbitExit;
0226 
0227             if NewVectorizedMethod
0228 
0229                 % very slow loop - use vector operations instead
0230                 vectind = BPMData.BPMIndex(1:NBPM);
0231 
0232                 % convert multidimensional loop to single matrix product
0233                 T3 = T([1, 3],:,vectind);
0234                 T2 = reshape(permute(T3, [1 3 2]), NBPM*2, 4);
0235 
0236                 % vector comparison
0237                 bgtc = find(vectind > CMData.HCMIndex(i));
0238                 bltc = find(vectind <= CMData.HCMIndex(i));
0239                 bgtc = [(bgtc - 1) .* 2 + 1 (bgtc - 1) .* 2 + 2];
0240                 bltc = [(bltc - 1) .* 2 + 1 (bltc - 1) .* 2 + 2];
0241 
0242                 % conditionally multiply
0243                 Tout1 = T2 * R0;
0244                 Tout2 = T2 * M44 * R0;
0245                 Tout = zeros(size(Tout1));
0246                 Tout(bgtc,:) = Tout1(bgtc,:);
0247                 Tout(bltc,:) = Tout2(bltc,:);
0248 
0249                 % interleave the writes
0250                 jjj = zeros(2, NBPM);
0251                 jjj(1,:) = 1:NBPM;
0252                 jjj(2,:) = NBPM+1:NBPM*2;
0253 
0254                 RM(jjj(:), i) = Tout;
0255 
0256             else
0257 
0258                 for j=1:NBPM
0259                     if BPMData.BPMIndex(j)>CMData.HCMIndex(i)
0260                         RM([j, j+NBPM],i) = T([1, 3],:,BPMData.BPMIndex(j))*R0;
0261                     else
0262                         RM([j, j+NBPM],i) = T([1, 3],:,BPMData.BPMIndex(j))*M44*R0;
0263                     end
0264                 end
0265 
0266             end
0267 
0268             if strcmpi(ClosedOrbitType, 'FixedPathLength')
0269                 % Use the average value of the dispersion at entrance and exit
0270                 D = HCORTheta(2,i) * ...
0271                     (Dispersion(1,CMData.HCMIndex(i))+Dispersion(1,CMData.HCMIndex(i)+1)) * ...
0272                     Dispersion([1 3],BPMData.BPMIndex) /L0/MCF/2;
0273 
0274                 RM(1:NBPM,i) = RM(1:NBPM,i) - D(1,:)';
0275                 RM(NBPM+1:end,i) = RM(NBPM+1:end,i) - D(2,:)';
0276             end
0277 
0278         end
0279 
0280         for i=1:NVC
0281             CI = CMData.VCMIndex(i);
0282 
0283             InverseT = inv(T(:,:,CI));
0284             OrbitEntrance = (inv(eye(4)-T(:,:,CI)*M44*InverseT) * T(:,:,CI) * M44 * ...
0285                 InverseT * (eye(4)+inv(M44VCOR{i}))*VCORTheta(:,i)/2);
0286             OrbitExit = VCORTheta(:,i)/2+M44VCOR{i}*(OrbitEntrance+VCORTheta(:,i)/2);
0287 
0288             R0 = inv(T(:,:,CI+1))*OrbitExit;
0289 
0290             if NewVectorizedMethod
0291 
0292                 % very slow loop - use vector operations instead
0293                 vectind = BPMData.BPMIndex(1:NBPM);
0294 
0295                 % convert multidimensional loop to single matrix product
0296                 T3 = T([1, 3],:,vectind);
0297                 T2 = reshape(permute(T3, [1 3 2]), NBPM*2, 4);
0298 
0299                 % vector comparison
0300                 bgtc = find(vectind > CMData.VCMIndex(i));
0301                 bltc = find(vectind <= CMData.VCMIndex(i));
0302                 bgtc = [(bgtc - 1) .* 2 + 1 (bgtc - 1) .* 2 + 2];
0303                 bltc = [(bltc - 1) .* 2 + 1 (bltc - 1) .* 2 + 2];
0304 
0305                 % conditionally multiply
0306                 Tout1 = T2 * R0;
0307                 Tout2 = T2 * M44 * R0;
0308                 Tout = zeros(size(Tout1));
0309                 Tout(bgtc,:) = Tout1(bgtc,:);
0310                 Tout(bltc,:) = Tout2(bltc,:);
0311 
0312                 % interleave the writes
0313                 jjj = zeros(2, NBPM);
0314                 jjj(1,:) = 1:NBPM;
0315                 jjj(2,:) = NBPM+1:NBPM*2;
0316 
0317                 RM(jjj(:), i+NHC) = Tout;
0318 
0319             else
0320 
0321                 for j=1:NBPM
0322                     if BPMData.BPMIndex(j)>CMData.VCMIndex(i)
0323                         RM([j, j+NBPM],i+NHC) = T([1, 3],:,BPMData.BPMIndex(j))*R0;
0324                     else
0325                         RM([j, j+NBPM],i+NHC) = T([1, 3],:,BPMData.BPMIndex(j))*M44*R0;
0326                     end
0327                 end
0328 
0329             end
0330 
0331             % Vertical correctors with coupling to X and non-zero horizontal dispersion
0332             if strcmpi(ClosedOrbitType, 'FixedPathLength')
0333                 % Use the average value of the dispersion at entrance and exit
0334                 D = VCORTheta(2,i)*(Dispersion(1,CMData.VCMIndex(i))+Dispersion(1,CMData.VCMIndex(i)+1))*...
0335                     Dispersion([1 3],BPMData.BPMIndex)/L0/MCF/2;
0336                 RM(1:NBPM,NHC+i) = RM(1:NBPM,NHC+i) - D(1,:)';
0337                 RM(NBPM+1:end,NHC+i) = RM(NBPM+1:end,NHC+i) - D(2,:)';
0338             end
0339         end
0340 
0341         if RFFLAG
0342             if strcmpi(DispersionMeasurement, 'Bidirectional')
0343                 ORBITPLUS = findsyncorbit(RINGData.Lattice, (-C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
0344                 ORBIT0    = findsyncorbit(RINGData.Lattice, ( C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
0345             else
0346                 ORBITPLUS = findsyncorbit(RINGData.Lattice, -C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2, 1:length(RINGData.Lattice)+1);
0347                 ORBIT0    = findsyncorbit(RINGData.Lattice, 0, 1:length(RINGData.Lattice)+1);
0348             end
0349             D = ORBITPLUS([1 3],BPMData.BPMIndex) - ORBIT0([1 3],BPMData.BPMIndex);
0350             RM(:,end) = [D(1,:)'; D(2,:)'];
0351         end
0352 
0353 
0354     elseif strcmpi(ClosedOrbitType, 'FixedPathLength')
0355         % Exact calculation using FINDSYNCORBIT
0356 
0357         for i = 1:NHC
0358             switch ATRING{CMData.HCMIndex(i)}.PassMethod
0359                 case 'CorrectorPass'
0360 
0361                     KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
0362 
0363                     if strcmpi(ResponseMatrixMeasurement, 'bidirectional')
0364                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i)/2;
0365                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
0366                         ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0367 
0368                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.HCMKicks(i)/2;
0369                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
0370                         ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0371                     else
0372                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i);
0373                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i);
0374                         ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0375 
0376                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1);
0377                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2);
0378                         ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0379                     end
0380 
0381                     ATRING{CMData.HCMIndex(i)}.KickAngle = KickAngle0;
0382 
0383                     RM(:,i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
0384 
0385                 case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
0386                     error('Not implemented yet');
0387                 otherwise
0388                     error('Unknown pass method for corrector');
0389             end
0390         end
0391 
0392         for i = 1:NVC
0393             switch ATRING{CMData.VCMIndex(i)}.PassMethod
0394                 case 'CorrectorPass'
0395                     KickAngle0 = ATRING{CMData.VCMIndex(i)}.KickAngle;
0396 
0397                     if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
0398                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i)/2;
0399                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
0400                         ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0401 
0402                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.VCMKicks(i)/2;
0403                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
0404                         ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0405                     else
0406                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i);
0407                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i);
0408                         ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0409 
0410                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2);
0411                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1);
0412                         ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
0413                     end
0414 
0415                     ATRING{CMData.VCMIndex(i)}.KickAngle = KickAngle0;
0416 
0417                     RM(:,NHC+i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
0418 
0419                 case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
0420                     error('Not implemented yet');
0421                 otherwise
0422                     error('Unknown pass method for corrector')
0423             end
0424         end
0425 
0426         if RFFLAG
0427             if strcmpi(DispersionMeasurement, 'bidirectional')
0428                 ORBITPLUS = findsyncorbit(RINGData.Lattice, (-C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
0429                 ORBIT0    = findsyncorbit(RINGData.Lattice, ( C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
0430             else
0431                 ORBITPLUS = findsyncorbit(RINGData.Lattice, -C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2, 1:length(RINGData.Lattice)+1);
0432                 ORBIT0    = findsyncorbit(RINGData.Lattice, 0, 1:length(RINGData.Lattice)+1);
0433             end
0434 
0435             D = ORBITPLUS([1 3],BPMData.BPMIndex) - ORBIT0([1 3],BPMData.BPMIndex);
0436             RM(:,end) = [D(1,:)';D(2,:)'];
0437         end
0438 
0439     elseif strcmpi(ClosedOrbitType, 'FixedMomentum')
0440         % ClosedOrbitType = 'fixedmomentum' - Exact calculation using FINDORBIT4
0441         for i = 1:NHC
0442             switch ATRING{CMData.HCMIndex(i)}.PassMethod
0443                 case 'CorrectorPass'
0444                     KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
0445 
0446                     if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
0447                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i)/2;
0448                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
0449                         ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0450 
0451                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.HCMKicks(i)/2;
0452                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
0453                         ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0454                     else
0455                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i);
0456                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i);
0457                         ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0458 
0459                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1);
0460                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2);
0461                         ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0462                     end
0463 
0464                     ATRING{CMData.HCMIndex(i)}.KickAngle = KickAngle0;
0465 
0466                     RM(:,i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
0467 
0468                 case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
0469                     error('Not implemented yet');
0470                 otherwise
0471                     error('Unknown pass method for corrector');
0472             end
0473         end
0474 
0475         for i = 1:NVC
0476             switch ATRING{CMData.HCMIndex(i)}.PassMethod
0477                 case 'CorrectorPass'
0478 
0479                     KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
0480 
0481                     if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
0482                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i)/2;
0483                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
0484                         ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0485 
0486                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.VCMKicks(i)/2;
0487                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
0488                         ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0489                     else
0490                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i);
0491                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i);
0492                         ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0493 
0494                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2);
0495                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1);
0496                         ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
0497                     end
0498 
0499                     ATRING{CMData.VCMIndex(i)}.KickAngle = KickAngle0;
0500 
0501                     RM(:,NHC+i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
0502 
0503                 case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
0504                     error('Not implemented yet');
0505                 otherwise
0506                     error('Unknown pass method for corrector')
0507             end
0508         end
0509 
0510         if RFFLAG
0511             if strcmpi(DispersionMeasurement, 'Bidirectional')
0512                 ORBITPLUS = findsyncorbit(RINGData.Lattice, (-C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
0513                 ORBIT0    = findsyncorbit(RINGData.Lattice, ( C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
0514             else
0515                 ORBITPLUS = findsyncorbit(RINGData.Lattice, -C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2, 1:length(RINGData.Lattice)+1);
0516                 ORBIT0    = findsyncorbit(RINGData.Lattice, 0, 1:length(RINGData.Lattice)+1);
0517             end
0518 
0519             D = ORBITPLUS([1 3],BPMData.BPMIndex) - ORBIT0([1 3],BPMData.BPMIndex);
0520             RM(:,end) = [D(1,:)'; D(2,:)'];
0521         end
0522 
0523     else
0524         error('ClosedOrbitType method unknown');
0525     end
0526 
0527     % End of storage ring
0528     
0529     
0530 else
0531     
0532     
0533     % Transport line
0534     if RFFLAG
0535         RFFLAG = 0;
0536         fprintf('   RF flag ignored for transport lines.\n');
0537     end
0538 
0539     %if strcmpi(ResponseMatrixCalculator, 'Linear')
0540     %    ResponseMatrixCalculator = 'Full';
0541     %    fprintf('   ResponseMatrixCalculator set to ''Full'' for transport lines (linear is actualy slower for small transport lines).\n');
0542     %end
0543 
0544     RM = zeros(2*NBPM,NHC+NVC);
0545     
0546     NE = length(ATRING);
0547 
0548     % Find the response matrix about an initial orbit
0549     if isfield(ATRING{1}, 'TwissData')
0550         TwissData = ATRING{1}.TwissData;
0551         R0 = [TwissData.ClosedOrbit; TwissData.dP; TwissData.dL];
0552     else
0553         % Try the middlelayer
0554         try
0555             TwissData = getfamilydata('TwissData');
0556             R0 = [TwissData.ClosedOrbit; TwissData.dP; TwissData.dL];
0557             if isempty(TwissData)
0558                 R0 = [0 0 0 0 0 0]';
0559             end
0560         catch
0561             R0 = [0 0 0 0 0 0]';
0562         end
0563     end
0564 
0565     if strcmpi(ResponseMatrixCalculator, 'Linear')
0566                        
0567         % Calculate linear optics and chromatic finctions for the model
0568         [M44, T, ClosedOrbit1] = findm44(ATRING, 0, 1:NE+1);
0569 
0570         % Don't use the closed orbit that findm44 use.  Use linepass.
0571         ClosedOrbit = linepass(ATRING, R0, 1:length(ATRING)+1);
0572     
0573         % Transfer matrixes through individual correctors
0574         M44HCOR = cell(1,NHC);
0575         M44VCOR = cell(1,NVC);
0576 
0577         for i=1:NHC
0578             M44HCOR{i} = findelemm44(ATRING{CMData.HCMIndex(i)}, ATRING{CMData.HCMIndex(i)}.PassMethod, ClosedOrbit(:,CMData.HCMIndex(i)));
0579         end
0580         for i=1:NVC
0581             match = find(CMData.VCMIndex(i)==CMData.HCMIndex);
0582             if match
0583                 M44VCOR{i} = M44HCOR{match};
0584             else
0585                 M44VCOR{i} = findelemm44(ATRING{CMData.VCMIndex(i)}, ATRING{CMData.VCMIndex(i)}.PassMethod, ClosedOrbit(:,CMData.VCMIndex(i)));
0586             end
0587         end
0588 
0589 
0590         % Assemble arrays of corrector kicks including coupling
0591         HCORTheta = zeros(4,NHC);
0592         VCORTheta = zeros(4,NVC);
0593 
0594         HCORTheta(2,:) = CMData.HCMKicks(:)';
0595         HCORTheta(4,:) = CMData.HCMCoupling(:)' .* CMData.HCMKicks(:)';
0596         VCORTheta(2,:) = CMData.VCMCoupling(:)' .* CMData.VCMKicks(:)';
0597         VCORTheta(4,:) = CMData.VCMKicks(:)';
0598 
0599 
0600         % Calculate closed orbit at the exit of each corrector magnet WITH applied kick
0601         for i = 1:NHC
0602             CI = CMData.HCMIndex(i);
0603             InverseT = inv(T(:,:,CI));
0604             
0605             % Split the kick on either side of the corrector
0606             OrbitExit = M44HCOR{i} * HCORTheta(:,i)/2 + HCORTheta(:,i)/2;
0607             OrbitEntrance = inv(M44HCOR{i}) * OrbitExit;
0608             R0 = InverseT * OrbitEntrance(1:4);
0609             
0610             for j = 1:NBPM
0611                 if BPMData.BPMIndex(j) > CMData.HCMIndex(i)
0612                     RM([j, j+NBPM],i) = T([1 3],:,BPMData.BPMIndex(j)) * R0;
0613                 end
0614             end
0615         end
0616 
0617         for i=1:NVC
0618             CI = CMData.VCMIndex(i);
0619             InverseT = inv(T(:,:,CI));
0620 
0621             % Split the kick on either side of the corrector
0622             OrbitExit = M44VCOR{i} * VCORTheta(:,i)/2 + VCORTheta(:,i)/2;
0623             OrbitEntrance = inv(M44HCOR{i}) * OrbitExit;
0624             R0 = InverseT * OrbitEntrance(1:4);
0625             
0626             for j=1:NBPM
0627                 if BPMData.BPMIndex(j)>CMData.VCMIndex(i)
0628                     RM([j, j+NBPM],i+NHC) = T([1 3],:,BPMData.BPMIndex(j)) * R0;
0629                 end
0630             end
0631         end
0632         
0633     elseif strcmpi(ResponseMatrixCalculator, 'Full')
0634         
0635         % Exact calculation using linepass
0636         for i = 1:NHC
0637             switch ATRING{CMData.HCMIndex(i)}.PassMethod
0638                 case 'CorrectorPass'
0639 
0640                     KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
0641 
0642                     if strcmpi(ResponseMatrixMeasurement, 'bidirectional')
0643                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i)/2;
0644                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
0645                         ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
0646 
0647                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.HCMKicks(i)/2;
0648                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
0649                         ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
0650                     else
0651                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i);
0652                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i);
0653                         ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
0654 
0655                         ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1);
0656                         ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2);
0657                         ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
0658                     end
0659 
0660                     ATRING{CMData.HCMIndex(i)}.KickAngle = KickAngle0;
0661 
0662                     RM(:,i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
0663 
0664                 case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
0665                     error('Not implemented yet');
0666                 otherwise
0667                     error('Unknown pass method for corrector');
0668             end
0669         end
0670 
0671         for i = 1:NVC
0672             switch ATRING{CMData.VCMIndex(i)}.PassMethod
0673                 case 'CorrectorPass'
0674                     KickAngle0 = ATRING{CMData.VCMIndex(i)}.KickAngle;
0675 
0676                     if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
0677                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i)/2;
0678                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
0679                         ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
0680 
0681                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.VCMKicks(i)/2;
0682                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
0683                         ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
0684                     else
0685                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i);
0686                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i);
0687                         ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
0688 
0689                         ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2);
0690                         ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1);
0691                         ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
0692                     end
0693 
0694                     ATRING{CMData.VCMIndex(i)}.KickAngle = KickAngle0;
0695 
0696                     RM(:,NHC+i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
0697 
0698                 case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
0699                     error('Not implemented yet');
0700                 otherwise
0701                     error('Unknown pass method for corrector')
0702             end
0703         end
0704     else
0705         error('Unknown method for transfer lines.');
0706     end
0707 end