source: MML/trunk/applications/loco/locoresponsematrix.m @ 4

Last change on this file since 4 was 4, checked in by zhangj, 10 years ago

Initial import--MML version from SOLEIL@2013
File size: 30.9 KB
Line 
1function RM = locoresponsematrix(RINGData, BPMData, CMData, varargin)
2%LOCORESPONSEMATRIX - Calculate the BPM response matrix and dispersion function
3% M = LOCORESPONCEMATRIX(RINGData, BPMData, CMData)
4%
5% Accelerator Toolbox implementation of generic LOCO function
6%
7% RINGData - must have fields 'Lattice', 'CavityFrequency', 'CavityHarmNumber'
8%            RINGData.Lattice - AT lattice cell arrary
9%            RINGData.CavityFrequency  [Hz]
10%            RINGData.CavityHarmNumber [Hz]
11%
12% CMData -   must have fields: 'HCMIndex', 'VCMIndex', 'HCMKicks', 'VCMKicks', 'HCMCoupling', 'VCMCoupling'
13%            CMData.HCMIndex    - indexes in the AT lattice of elements used as horizontal correctors
14%            CMData.VCMIndex    - indexes in the AT lattice of elements used as vertical correctors
15%                                 Elements used as correctors in both planes should be included in both lists
16%            CMData.HCMKicks    - kick  size [radians] of horizontal correctors in the horizontal plane
17%            CMData.VCMKicks    - kick  size [radians] of vertical correctors in the vertical plane
18%            CMData.HCMCoupling - corrector coupling coefficient into another plane:
19%                                 0.01 coupling means that for 1e-3 kick in the horizontal direction there
20%                                 is a 1e-5 rad kick in the vertical direction
21%            CMData.VCMCoupling - corrector coupling coefficient into another plane:
22%                                 0.01 coupling means that for 1e-3 kick in the vertical direction there
23%                                 is a 1e-5 rad kick in the horizontal direction
24%
25% BPMData -  must have field 'BPMIndex'
26%            CMData.BPMIndex - indexes of all BPMs or observation points in the AT lattice
27%                              All BPS and observation points (single plane too)
28%                              are included in CMData.BPMIndex.
29%
30% Return value: a matrix with number of rows equal to 2*length(CMData.BPMIndex) and the number of columns
31%               equal length(CMData.HCMIndex)+length(CMData.VCMIndex)
32%
33% Additional string flags (in any order)
34%
35% LOCORESPONSEMATRIX(...,ClosedOrbitType,...)
36%       ClosedOrbitType is 'fixedmomentum',  'fixedpathlength' (default)
37%
38% LOCORESPONCEMATRIX(..., 'linear') calculates M using linear approximation  !!! including the dispersion terms
39%
40% LOCORESPONCEMATRIX(..., 'RF', DeltaRF) - 'RF' switch must be followed by the value of DeltaRF [Hz]
41%
42% LOCORESPONCEMATRIX(..., 'ResponseMatrixMeasurement', 'oneway') - 'oneway' switch is used
43%                        when the response matrix was measured only kicking the i-th corrector
44%                        to +KicksCoupled(i) one way, default: ResponseMatrixMeasurement = 'bidirectional'
45%
46% LOCORESPONCEMATRIX(..., 'DispersionMeasurement', 'oneway') - 'oneway' switch is used
47%                        when the dispersion was measured only by varying the
48%                        RF frequency in one direction, default: DispersionMeasurement = 'bidirectional'
49%
50% Or a Flags structure can be an input argument:
51% LOCORESPONCEMATRIX(..., Flags)
52%    Flags.ResponseMatrixMeasurement = 'oneway' or {'bi-directional'}
53%    Flags.DispersionMeasurement     = 'oneway' or {'bi-directional'}
54%    Flags.ResponseMatrixCalculator  = {'linear'} or 'full'
55%    Flags.ClosedOrbitType           = 'fixedmomentum' or {'fixedpathlength'}
56%
57%  NOTE
58%  1. Flag names are not case sensitive
59%  2. JR - 12/7/05 added vectorized linear rm calculation (NewVectorizedMethod = 1)
60
61NewVectorizedMethod = 1;
62
63C = 2.99792458e8;
64
65% Defaults
66ResponseMatrixMeasurement = 'bidirectional';
67DispersionMeasurement     = 'bidirectional';
68ResponseMatrixCalculator  = 'linear';
69ClosedOrbitType           = 'fixedpathlength';
70MachineType               = 'StorageRing';
71
72RFFLAG = 0;
73DeltaRF = [];
74
75N = nargin-3;
76i = 0;
77while i < N
78    i = i + 1;
79    if isstruct(varargin{i})
80        Flags = varargin{i};
81        if isfield(Flags,'ResponseMatrixCalculator')
82            ResponseMatrixCalculator = Flags.ResponseMatrixCalculator;
83        end
84        if isfield(Flags,'ClosedOrbitType')
85            ClosedOrbitType = Flags.ClosedOrbitType;
86        end
87        if isfield(Flags,'ResponseMatrixMeasurement')
88            ResponseMatrixMeasurement = Flags.ResponseMatrixMeasurement;
89        end
90        if isfield(Flags,'DispersionMeasurement')
91            DispersionMeasurement = Flags.DispersionMeasurement;
92        end
93        if isfield(Flags,'MachineType')
94            MachineType = Flags.MachineType;
95        end
96    elseif ischar(varargin{i})
97        switch lower(varargin{i})
98            case 'linear'
99                ResponseMatrixCalculator = 'linear';
100            case 'full'
101                ResponseMatrixCalculator = 'full';
102            case 'fixedmomentum'
103                ClosedOrbitType = 'fixedmomentum';
104            case 'fixedpathlength'
105                ClosedOrbitType = 'fixedpathlength';
106            case 'rf'
107                if (i+4<=nargin) && isnumeric(varargin{i+1})
108                    RFFLAG = 1;
109                    DeltaRF = varargin{i+1};
110                    i = i + 1;
111                else
112                    error('''RF'' flag must be followed by a numeric value of delta RF [Hz]');
113                end
114            case 'dispersionmeasurement'
115                if (i+4<=nargin) && ischar(varargin{i+1})
116                    DispersionMeasurement = varargin{i+1};
117                    i = i + 1;
118                else
119                    error('''DispersionMeasurement'' flag must be followed by ''oneway'' or ''bidirectional''');
120                end
121            case 'responsematrixmeasurement'
122                if (i+4<=nargin) && ischar(varargin{i+1})
123                    ResponseMatrixMeasurement = varargin{i+1};
124                    i = i + 1;
125                else
126                    error('''ResponseMatrixMeasurement'' flag must be followed by ''oneway'' or ''bidirectional''');
127                end
128            case {'storagering','booster','boosterring'}
129                MachineType = 'StorageRing';
130            case {'transport','linac'}
131                MachineType = 'Transport';
132            otherwise
133                warning('Unknown switch ignored.');
134        end
135    else
136        warning('Unknown switch ignored.');
137    end
138end
139
140
141% Input checks
142if ~strcmpi(ResponseMatrixMeasurement, 'bidirectional') && ~strcmpi(ResponseMatrixMeasurement, 'oneway')
143    error('Unknown ResponseMatrixMeasurement type');
144end
145if ~strcmpi(DispersionMeasurement, 'bidirectional') && ~strcmpi(DispersionMeasurement, 'oneway')
146    error('Unknown DispersionMeasurement type');
147end
148if ~strcmpi(ResponseMatrixCalculator, 'linear') && ~strcmpi(ResponseMatrixCalculator, 'full')
149    error('Unknown ResponseMatrixCalculator method');
150end
151if ~strcmpi(ClosedOrbitType, 'fixedpathlength') && ~strcmpi(ClosedOrbitType, 'fixedmomentum')
152    error('Unknown ClosedOrbitType method');
153end
154if ~strcmpi(MachineType, 'StorageRing') && ~strcmpi(MachineType, 'Transport')
155    error('Unknown MachineType');
156end
157
158
159% Initialize
160ATRING = RINGData.Lattice;
161NHC = length(CMData.HCMIndex);
162NVC = length(CMData.VCMIndex);
163NBPM = length(BPMData.BPMIndex);
164
165
166if strcmpi(MachineType, 'StorageRing')
167
168    if RFFLAG
169        % Add an extra column for the orbit responce to the RF frequency change
170        RM = zeros(2*NBPM,NHC+NVC+1);
171    else
172        RM = zeros(2*NBPM,NHC+NVC);
173    end
174
175    NE = length(ATRING);
176    if strcmpi(ResponseMatrixCalculator, 'Linear')
177        % Calculate linear optics and chromatic finctions for the model
178        [M44,T,ClosedOrbit] = findm44(ATRING,0,1:NE+1);
179        DP = 0.00001;
180        ClosedOrbitDP = findorbit4(ATRING,DP,1:NE+1);
181        Dispersion = (ClosedOrbitDP-ClosedOrbit)/DP;
182        L0 = findspos(ATRING,NE+1);
183
184        %%X1(6)/(DP*L0); % is not the same as mcf(ATRING)???  G. Portmann
185        %X1 = ringpass(ATRING,[ClosedOrbitDP(:,1);DP;0]);
186        %MCF = X1(6)/(DP*L0);
187        MCF = mcf(ATRING);
188
189        % Transfer matrixes through individual correctors
190        M44HCOR = cell(1,NHC);
191        M44VCOR = cell(1,NVC);
192        for i=1:NHC
193            M44HCOR{i}=findelemm44(ATRING{CMData.HCMIndex(i)},ATRING{CMData.HCMIndex(i)}.PassMethod,[ClosedOrbit(:,CMData.HCMIndex(i));0;0]);
194        end
195        for i=1:NVC
196            match = find(CMData.VCMIndex(i)==CMData.HCMIndex);
197            if match
198                M44VCOR{i}=M44HCOR{match};
199            else
200                M44VCOR{i}=findelemm44(ATRING{CMData.VCMIndex(i)},ATRING{CMData.VCMIndex(i)}.PassMethod,[ClosedOrbit(:,CMData.VCMIndex(i));0;0]);
201            end
202        end
203
204
205        % Assemble arrays of corrector kicks including coupling
206        HCORTheta = zeros(4,NHC);
207        VCORTheta = zeros(4,NVC);
208
209        HCORTheta(2,:) = CMData.HCMKicks(:)';
210        HCORTheta(4,:) = CMData.HCMCoupling(:)' .* CMData.HCMKicks(:)';
211        VCORTheta(2,:) = CMData.VCMCoupling(:)' .* CMData.VCMKicks(:)';
212        VCORTheta(4,:) = CMData.VCMKicks(:)';
213
214
215        % Calculate closed orbit at the exit of each corrector magnet WITH applied kick
216        for i=1:NHC
217            CI = CMData.HCMIndex(i);
218            InverseT = inv(T(:,:,CI));
219            OrbitEntrance = (inv(eye(4)-T(:,:,CI)*M44*InverseT)*...
220                T(:,:,CI)*M44*InverseT*(eye(4)+inv(M44HCOR{i}))*HCORTheta(:,i)/2);
221
222            OrbitExit = HCORTheta(:,i)/2+M44HCOR{i}*(OrbitEntrance+HCORTheta(:,i)/2);
223
224
225            R0 = inv(T(:,:,CI+1))*OrbitExit;
226
227            if NewVectorizedMethod
228
229                % very slow loop - use vector operations instead
230                vectind = BPMData.BPMIndex(1:NBPM);
231
232                % convert multidimensional loop to single matrix product
233                T3 = T([1, 3],:,vectind);
234                T2 = reshape(permute(T3, [1 3 2]), NBPM*2, 4);
235
236                % vector comparison
237                bgtc = find(vectind > CMData.HCMIndex(i));
238                bltc = find(vectind <= CMData.HCMIndex(i));
239                bgtc = [(bgtc - 1) .* 2 + 1 (bgtc - 1) .* 2 + 2];
240                bltc = [(bltc - 1) .* 2 + 1 (bltc - 1) .* 2 + 2];
241
242                % conditionally multiply
243                Tout1 = T2 * R0;
244                Tout2 = T2 * M44 * R0;
245                Tout = zeros(size(Tout1));
246                Tout(bgtc,:) = Tout1(bgtc,:);
247                Tout(bltc,:) = Tout2(bltc,:);
248
249                % interleave the writes
250                jjj = zeros(2, NBPM);
251                jjj(1,:) = 1:NBPM;
252                jjj(2,:) = NBPM+1:NBPM*2;
253
254                RM(jjj(:), i) = Tout;
255
256            else
257
258                for j=1:NBPM
259                    if BPMData.BPMIndex(j)>CMData.HCMIndex(i)
260                        RM([j, j+NBPM],i) = T([1, 3],:,BPMData.BPMIndex(j))*R0;
261                    else
262                        RM([j, j+NBPM],i) = T([1, 3],:,BPMData.BPMIndex(j))*M44*R0;
263                    end
264                end
265
266            end
267
268            if strcmpi(ClosedOrbitType, 'FixedPathLength')
269                % Use the average value of the dispersion at entrance and exit
270                D = HCORTheta(2,i) * ...
271                    (Dispersion(1,CMData.HCMIndex(i))+Dispersion(1,CMData.HCMIndex(i)+1)) * ...
272                    Dispersion([1 3],BPMData.BPMIndex) /L0/MCF/2;
273
274                RM(1:NBPM,i) = RM(1:NBPM,i) - D(1,:)';
275                RM(NBPM+1:end,i) = RM(NBPM+1:end,i) - D(2,:)';
276            end
277
278        end
279
280        for i=1:NVC
281            CI = CMData.VCMIndex(i);
282
283            InverseT = inv(T(:,:,CI));
284            OrbitEntrance = (inv(eye(4)-T(:,:,CI)*M44*InverseT) * T(:,:,CI) * M44 * ...
285                InverseT * (eye(4)+inv(M44VCOR{i}))*VCORTheta(:,i)/2);
286            OrbitExit = VCORTheta(:,i)/2+M44VCOR{i}*(OrbitEntrance+VCORTheta(:,i)/2);
287
288            R0 = inv(T(:,:,CI+1))*OrbitExit;
289
290            if NewVectorizedMethod
291
292                % very slow loop - use vector operations instead
293                vectind = BPMData.BPMIndex(1:NBPM);
294
295                % convert multidimensional loop to single matrix product
296                T3 = T([1, 3],:,vectind);
297                T2 = reshape(permute(T3, [1 3 2]), NBPM*2, 4);
298
299                % vector comparison
300                bgtc = find(vectind > CMData.VCMIndex(i));
301                bltc = find(vectind <= CMData.VCMIndex(i));
302                bgtc = [(bgtc - 1) .* 2 + 1 (bgtc - 1) .* 2 + 2];
303                bltc = [(bltc - 1) .* 2 + 1 (bltc - 1) .* 2 + 2];
304
305                % conditionally multiply
306                Tout1 = T2 * R0;
307                Tout2 = T2 * M44 * R0;
308                Tout = zeros(size(Tout1));
309                Tout(bgtc,:) = Tout1(bgtc,:);
310                Tout(bltc,:) = Tout2(bltc,:);
311
312                % interleave the writes
313                jjj = zeros(2, NBPM);
314                jjj(1,:) = 1:NBPM;
315                jjj(2,:) = NBPM+1:NBPM*2;
316
317                RM(jjj(:), i+NHC) = Tout;
318
319            else
320
321                for j=1:NBPM
322                    if BPMData.BPMIndex(j)>CMData.VCMIndex(i)
323                        RM([j, j+NBPM],i+NHC) = T([1, 3],:,BPMData.BPMIndex(j))*R0;
324                    else
325                        RM([j, j+NBPM],i+NHC) = T([1, 3],:,BPMData.BPMIndex(j))*M44*R0;
326                    end
327                end
328
329            end
330
331            % Vertical correctors with coupling to X and non-zero horizontal dispersion
332            if strcmpi(ClosedOrbitType, 'FixedPathLength')
333                % Use the average value of the dispersion at entrance and exit
334                D = VCORTheta(2,i)*(Dispersion(1,CMData.VCMIndex(i))+Dispersion(1,CMData.VCMIndex(i)+1))*...
335                    Dispersion([1 3],BPMData.BPMIndex)/L0/MCF/2;
336                RM(1:NBPM,NHC+i) = RM(1:NBPM,NHC+i) - D(1,:)';
337                RM(NBPM+1:end,NHC+i) = RM(NBPM+1:end,NHC+i) - D(2,:)';
338            end
339        end
340
341        if RFFLAG
342            if strcmpi(DispersionMeasurement, 'Bidirectional')
343                ORBITPLUS = findsyncorbit(RINGData.Lattice, (-C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
344                ORBIT0    = findsyncorbit(RINGData.Lattice, ( C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
345            else
346                ORBITPLUS = findsyncorbit(RINGData.Lattice, -C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2, 1:length(RINGData.Lattice)+1);
347                ORBIT0    = findsyncorbit(RINGData.Lattice, 0, 1:length(RINGData.Lattice)+1);
348            end
349            D = ORBITPLUS([1 3],BPMData.BPMIndex) - ORBIT0([1 3],BPMData.BPMIndex);
350            RM(:,end) = [D(1,:)'; D(2,:)'];
351        end
352
353
354    elseif strcmpi(ClosedOrbitType, 'FixedPathLength')
355        % Exact calculation using FINDSYNCORBIT
356
357        for i = 1:NHC
358            switch ATRING{CMData.HCMIndex(i)}.PassMethod
359                case 'CorrectorPass'
360
361                    KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
362
363                    if strcmpi(ResponseMatrixMeasurement, 'bidirectional')
364                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i)/2;
365                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
366                        ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
367
368                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.HCMKicks(i)/2;
369                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
370                        ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
371                    else
372                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i);
373                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i);
374                        ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
375
376                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1);
377                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2);
378                        ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
379                    end
380
381                    ATRING{CMData.HCMIndex(i)}.KickAngle = KickAngle0;
382
383                    RM(:,i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
384
385                case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
386                    error('Not implemented yet');
387                otherwise
388                    error('Unknown pass method for corrector');
389            end
390        end
391
392        for i = 1:NVC
393            switch ATRING{CMData.VCMIndex(i)}.PassMethod
394                case 'CorrectorPass'
395                    KickAngle0 = ATRING{CMData.VCMIndex(i)}.KickAngle;
396
397                    if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
398                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i)/2;
399                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
400                        ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
401
402                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.VCMKicks(i)/2;
403                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
404                        ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
405                    else
406                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i);
407                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i);
408                        ORBITPLUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
409
410                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2);
411                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1);
412                        ORBITMINUS = findsyncorbit(ATRING,0,BPMData.BPMIndex);
413                    end
414
415                    ATRING{CMData.VCMIndex(i)}.KickAngle = KickAngle0;
416
417                    RM(:,NHC+i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
418
419                case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
420                    error('Not implemented yet');
421                otherwise
422                    error('Unknown pass method for corrector')
423            end
424        end
425
426        if RFFLAG
427            if strcmpi(DispersionMeasurement, 'bidirectional')
428                ORBITPLUS = findsyncorbit(RINGData.Lattice, (-C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
429                ORBIT0    = findsyncorbit(RINGData.Lattice, ( C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
430            else
431                ORBITPLUS = findsyncorbit(RINGData.Lattice, -C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2, 1:length(RINGData.Lattice)+1);
432                ORBIT0    = findsyncorbit(RINGData.Lattice, 0, 1:length(RINGData.Lattice)+1);
433            end
434
435            D = ORBITPLUS([1 3],BPMData.BPMIndex) - ORBIT0([1 3],BPMData.BPMIndex);
436            RM(:,end) = [D(1,:)';D(2,:)'];
437        end
438
439    elseif strcmpi(ClosedOrbitType, 'FixedMomentum')
440        % ClosedOrbitType = 'fixedmomentum' - Exact calculation using FINDORBIT4
441        for i = 1:NHC
442            switch ATRING{CMData.HCMIndex(i)}.PassMethod
443                case 'CorrectorPass'
444                    KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
445
446                    if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
447                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i)/2;
448                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
449                        ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
450
451                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.HCMKicks(i)/2;
452                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
453                        ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
454                    else
455                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i);
456                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i);
457                        ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
458
459                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1);
460                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2);
461                        ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
462                    end
463
464                    ATRING{CMData.HCMIndex(i)}.KickAngle = KickAngle0;
465
466                    RM(:,i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
467
468                case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
469                    error('Not implemented yet');
470                otherwise
471                    error('Unknown pass method for corrector');
472            end
473        end
474
475        for i = 1:NVC
476            switch ATRING{CMData.HCMIndex(i)}.PassMethod
477                case 'CorrectorPass'
478
479                    KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
480
481                    if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
482                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i)/2;
483                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
484                        ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
485
486                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.VCMKicks(i)/2;
487                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
488                        ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
489                    else
490                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i);
491                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i);
492                        ORBITPLUS = findorbit4(ATRING,0,BPMData.BPMIndex);
493
494                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2);
495                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1);
496                        ORBITMINUS = findorbit4(ATRING,0,BPMData.BPMIndex);
497                    end
498
499                    ATRING{CMData.VCMIndex(i)}.KickAngle = KickAngle0;
500
501                    RM(:,NHC+i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
502
503                case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
504                    error('Not implemented yet');
505                otherwise
506                    error('Unknown pass method for corrector')
507            end
508        end
509
510        if RFFLAG
511            if strcmpi(DispersionMeasurement, 'Bidirectional')
512                ORBITPLUS = findsyncorbit(RINGData.Lattice, (-C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
513                ORBIT0    = findsyncorbit(RINGData.Lattice, ( C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2)/2, 1:length(RINGData.Lattice)+1);
514            else
515                ORBITPLUS = findsyncorbit(RINGData.Lattice, -C*DeltaRF*RINGData.CavityHarmNumber/RINGData.CavityFrequency^2, 1:length(RINGData.Lattice)+1);
516                ORBIT0    = findsyncorbit(RINGData.Lattice, 0, 1:length(RINGData.Lattice)+1);
517            end
518
519            D = ORBITPLUS([1 3],BPMData.BPMIndex) - ORBIT0([1 3],BPMData.BPMIndex);
520            RM(:,end) = [D(1,:)'; D(2,:)'];
521        end
522
523    else
524        error('ClosedOrbitType method unknown');
525    end
526
527    % End of storage ring
528   
529   
530else
531   
532   
533    % Transport line
534    if RFFLAG
535        RFFLAG = 0;
536        fprintf('   RF flag ignored for transport lines.\n');
537    end
538
539    %if strcmpi(ResponseMatrixCalculator, 'Linear')
540    %    ResponseMatrixCalculator = 'Full';
541    %    fprintf('   ResponseMatrixCalculator set to ''Full'' for transport lines (linear is actualy slower for small transport lines).\n');
542    %end
543
544    RM = zeros(2*NBPM,NHC+NVC);
545   
546    NE = length(ATRING);
547
548    % Find the response matrix about an initial orbit
549    if isfield(ATRING{1}, 'TwissData')
550        TwissData = ATRING{1}.TwissData;
551        R0 = [TwissData.ClosedOrbit; TwissData.dP; TwissData.dL];
552    else
553        % Try the middlelayer
554        try
555            TwissData = getfamilydata('TwissData');
556            R0 = [TwissData.ClosedOrbit; TwissData.dP; TwissData.dL];
557            if isempty(TwissData)
558                R0 = [0 0 0 0 0 0]';
559            end
560        catch
561            R0 = [0 0 0 0 0 0]';
562        end
563    end
564
565    if strcmpi(ResponseMatrixCalculator, 'Linear')
566                       
567        % Calculate linear optics and chromatic finctions for the model
568        [M44, T, ClosedOrbit1] = findm44(ATRING, 0, 1:NE+1);
569
570        % Don't use the closed orbit that findm44 use.  Use linepass.
571        ClosedOrbit = linepass(ATRING, R0, 1:length(ATRING)+1);
572   
573        % Transfer matrixes through individual correctors
574        M44HCOR = cell(1,NHC);
575        M44VCOR = cell(1,NVC);
576
577        for i=1:NHC
578            M44HCOR{i} = findelemm44(ATRING{CMData.HCMIndex(i)}, ATRING{CMData.HCMIndex(i)}.PassMethod, ClosedOrbit(:,CMData.HCMIndex(i)));
579        end
580        for i=1:NVC
581            match = find(CMData.VCMIndex(i)==CMData.HCMIndex);
582            if match
583                M44VCOR{i} = M44HCOR{match};
584            else
585                M44VCOR{i} = findelemm44(ATRING{CMData.VCMIndex(i)}, ATRING{CMData.VCMIndex(i)}.PassMethod, ClosedOrbit(:,CMData.VCMIndex(i)));
586            end
587        end
588
589
590        % Assemble arrays of corrector kicks including coupling
591        HCORTheta = zeros(4,NHC);
592        VCORTheta = zeros(4,NVC);
593
594        HCORTheta(2,:) = CMData.HCMKicks(:)';
595        HCORTheta(4,:) = CMData.HCMCoupling(:)' .* CMData.HCMKicks(:)';
596        VCORTheta(2,:) = CMData.VCMCoupling(:)' .* CMData.VCMKicks(:)';
597        VCORTheta(4,:) = CMData.VCMKicks(:)';
598
599
600        % Calculate closed orbit at the exit of each corrector magnet WITH applied kick
601        for i = 1:NHC
602            CI = CMData.HCMIndex(i);
603            InverseT = inv(T(:,:,CI));
604           
605            % Split the kick on either side of the corrector
606            OrbitExit = M44HCOR{i} * HCORTheta(:,i)/2 + HCORTheta(:,i)/2;
607            OrbitEntrance = inv(M44HCOR{i}) * OrbitExit;
608            R0 = InverseT * OrbitEntrance(1:4);
609           
610            for j = 1:NBPM
611                if BPMData.BPMIndex(j) > CMData.HCMIndex(i)
612                    RM([j, j+NBPM],i) = T([1 3],:,BPMData.BPMIndex(j)) * R0;
613                end
614            end
615        end
616
617        for i=1:NVC
618            CI = CMData.VCMIndex(i);
619            InverseT = inv(T(:,:,CI));
620
621            % Split the kick on either side of the corrector
622            OrbitExit = M44VCOR{i} * VCORTheta(:,i)/2 + VCORTheta(:,i)/2;
623            OrbitEntrance = inv(M44HCOR{i}) * OrbitExit;
624            R0 = InverseT * OrbitEntrance(1:4);
625           
626            for j=1:NBPM
627                if BPMData.BPMIndex(j)>CMData.VCMIndex(i)
628                    RM([j, j+NBPM],i+NHC) = T([1 3],:,BPMData.BPMIndex(j)) * R0;
629                end
630            end
631        end
632       
633    elseif strcmpi(ResponseMatrixCalculator, 'Full')
634       
635        % Exact calculation using linepass
636        for i = 1:NHC
637            switch ATRING{CMData.HCMIndex(i)}.PassMethod
638                case 'CorrectorPass'
639
640                    KickAngle0 = ATRING{CMData.HCMIndex(i)}.KickAngle;
641
642                    if strcmpi(ResponseMatrixMeasurement, 'bidirectional')
643                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i)/2;
644                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
645                        ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
646
647                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.HCMKicks(i)/2;
648                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.HCMKicks(i)*CMData.HCMCoupling(i)/2;
649                        ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
650                    else
651                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.HCMKicks(i);
652                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.HCMKicks(i)*CMData.HCMCoupling(i);
653                        ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
654
655                        ATRING{CMData.HCMIndex(i)}.KickAngle(1) = KickAngle0(1);
656                        ATRING{CMData.HCMIndex(i)}.KickAngle(2) = KickAngle0(2);
657                        ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
658                    end
659
660                    ATRING{CMData.HCMIndex(i)}.KickAngle = KickAngle0;
661
662                    RM(:,i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
663
664                case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
665                    error('Not implemented yet');
666                otherwise
667                    error('Unknown pass method for corrector');
668            end
669        end
670
671        for i = 1:NVC
672            switch ATRING{CMData.VCMIndex(i)}.PassMethod
673                case 'CorrectorPass'
674                    KickAngle0 = ATRING{CMData.VCMIndex(i)}.KickAngle;
675
676                    if strcmpi(ResponseMatrixMeasurement, 'Bidirectional')
677                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i)/2;
678                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
679                        ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
680
681                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) - CMData.VCMKicks(i)/2;
682                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) - CMData.VCMKicks(i)*CMData.VCMCoupling(i)/2;
683                        ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
684                    else
685                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2) + CMData.VCMKicks(i);
686                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1) + CMData.VCMKicks(i)*CMData.VCMCoupling(i);
687                        ORBITPLUS = linepass(ATRING,R0,BPMData.BPMIndex);
688
689                        ATRING{CMData.VCMIndex(i)}.KickAngle(2) = KickAngle0(2);
690                        ATRING{CMData.VCMIndex(i)}.KickAngle(1) = KickAngle0(1);
691                        ORBITMINUS = linepass(ATRING,R0,BPMData.BPMIndex);
692                    end
693
694                    ATRING{CMData.VCMIndex(i)}.KickAngle = KickAngle0;
695
696                    RM(:,NHC+i) = [ORBITPLUS(1,:)-ORBITMINUS(1,:),ORBITPLUS(3,:)-ORBITMINUS(3,:)]';
697
698                case {'StrMPoleSymplectic4Pass','BndMPoleSymplectic4Pass'}
699                    error('Not implemented yet');
700                otherwise
701                    error('Unknown pass method for corrector')
702            end
703        end
704    else
705        error('Unknown method for transfer lines.');
706    end
707end
Note: See TracBrowser for help on using the repository browser.