source: MML/trunk/mml/at/modeldisp.m @ 4

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

Initial import--MML version from SOLEIL@2013
File size: 20.6 KB
Line 
1function [Dx, Dy, Sx, Sy] = modeldisp(varargin)
2%MODELDISP - Returns the dispersion function of the model
3%  [Dx, Dy, Sx, Sy] = modeldisp(DeltaRF, Family1, DeviceList1, Family2, DeviceList2, ModulationMethod)
4%  [Dx, Dy, Sx, Sy] = modeldisp(DeltaRF, Family1, DeviceList1);  % Uses the same information for both planes
5%  [Dx, Dy, Sx, Sy] = modeldisp(DeltaRF, Family1, Family2);      % Assumes DeviceList1 and DeviceList2 are []
6%  [Dx, Dy, Sx, Sy] = modeldisp(Family1, Family2);               % Assumes DeviceList1 and DeviceList2 are []
7%  [Dx, Dy, Sx, Sy] = modeldisp(DeltaRF, Family1);               % Family2 = Family1
8%  [Dx, Dy, Sx, Sy] = modeldisp;                                 % DeltaRF, Family, DeviceList are all optional
9%
10%  INPUTS
11%  1. DeltaRF is the change in RF frequency {Default: .2% energy change}
12%  2. Family1 and Family2 are the family names for where to measure the horizontal/vertical dispersion
13%     A family name can be a middlelayer family, an AT family (FamName), or 'All' for all AT elements
14%     plus the end of the ring.  However, if using a non-middlelayer family, 'Physics' units must be used.
15%     {Default or []: 'All'}
16%  3. DeviceList1 and DeviceList2 are the device list corresponding to Family2 and Family2.
17%     Or a subindex array of the AT family.
18%     {Default or []: the entire list}
19%  4. ModulationMethod - Method for changing the ActuatorFamily
20%                       'bipolar' changes the RF by +/- DeltaRF/2 {Default}
21%                       'unipolar' changes the RF from 0 to DeltaRF
22%  5. Units - {Default: 'Physics'}
23%     'Physics' - Forces dispersion units of meters/(dp/p) add 'Physics' with an optional input of the
24%                 momentum compaction factor (mcf will be found from the AT function mcf if empty). 
25%     'Hardware' - Forces dispersion units of (Hardware BPM Units)/(Hardware RF units)
26%  6. 'Struct'  will return data structures instead of vectors
27%     'Numeric' will return vector outputs {Default}
28%  7. Optional display
29%     'Display'   - Plot the dispersion {Default if no outputs exist}
30%     'NoDisplay' - Dispersion will not be plotted {Default if outputs exist}
31%  8. 'NoArchive' - No file archive {Default}
32%     'Archive'   - Save a dispersion data structure to \<Directory.DispData>\<DispArchiveFile><Date><Time>.mat
33%                   To change the filename, included the filename after the 'Archive', '' to browse
34%
35%  OUTPUTS
36%  1. Dx and Dy - Horizontal and vertical dispersion function
37%                 Note: it does not matter what the input families are
38%                       Dx is always horizontal and Dy vertical
39%  2. Sx and Sy are longitudinal locations in the ring [meters].
40%
41%  For hardware units:
42%  Dx = Delta BPMx / Delta RF and Dy = Delta BPMy / Delta RF
43%       hence Dx and Dy are not quite the definition of dispersion
44%
45%               x2(RF0+DeltaRF/2) - x1(RF0-DeltaRF/2)
46%           D = -------------------------------------
47%                              DeltaRF
48%           
49%           where RF0 = is the present RF frequency
50%
51%  For physics units: DeltaRF is converted to change in energy, dp/p,
52%  hence dispersion is has units of [meters/(dp/p)]
53%   
54%  NOTE
55%  1. This function uses the model coordinate system, ie, no BPM rolls.  If hardware
56%     units are used, only a gain change will be applied.
57%  2. This function only call the AT model.
58%  3. If no output exist, the dispersion function will be plotted to the screen.
59%  4. Family2 and DeviceList1 can be any family.  For instance, if Family2='VCM'
60%     and DeviceList1=[], then Dx is the horizontal dispersion function at the
61%     vertical corrector magnets (similarly for Family2 and DeviceList2).
62%
63%  Written by Greg Portmann
64
65
66global THERING
67if isempty(THERING)
68    error('Simulator variable is not setup properly');
69end
70
71DeltaRFDefault = .1;  % Hz
72Family1 = 'ALL';
73Family2 = 'ALL';
74%Family1 = 'BPMx';
75%Family2 = 'BPMy';
76DeviceList1 = [];   
77DeviceList2 = [];   
78ModulationMethod = 'bipolar';
79StructOutputFlag = 0;
80NumericOutputFlag = 0;
81ArchiveFlag = 0;
82FileName = -1;
83if nargout == 0
84    DisplayFlag = 1;
85else
86    DisplayFlag = 0;
87end
88ModeFlag = {};  % model, online, manual, or '' for default mode
89UnitsFlag = ''; % hardware, physics, or '' for default units
90MCF = [];
91
92InputFlags = {};
93for i = length(varargin):-1:1
94    if isstruct(varargin{i})
95        % Ignor structures
96    elseif iscell(varargin{i})
97        % Ignor cells
98    elseif strcmpi(varargin{i},'struct')
99        StructOutputFlag = 1;
100        varargin(i) = [];
101    elseif strcmpi(varargin{i},'numeric')
102        NumericOutputFlag = 1;
103        StructOutputFlag = 0;
104        varargin(i) = [];
105    elseif strcmpi(varargin{i},'archive')
106        ArchiveFlag = 1;
107        if length(varargin) > i
108            % Look for a filename as the next input
109            if ischar(varargin{i+1})
110                FileName = varargin{i+1};
111                varargin(i+1) = [];
112            end
113        end
114        varargin(i) = [];
115    elseif strcmpi(varargin{i},'noarchive')
116        ArchiveFlag = 0;
117        varargin(i) = [];
118    elseif strcmpi(varargin{i},'Display')
119        DisplayFlag = 1;
120        varargin(i) = [];
121    elseif strcmpi(varargin{i},'NoDisplay')
122        DisplayFlag = 0;
123        varargin(i) = [];
124    elseif strcmpi(varargin{i},'bipolar')
125        ModulationMethod = 'bipolar';
126        varargin(i) = [];
127    elseif strcmpi(varargin{i},'unipolar')
128        ModulationMethod = 'unipolar';
129        varargin(i) = [];
130    elseif strcmpi(varargin{i},'eta') | strcmpi(varargin{i},'physics')
131        UnitsFlag = 'Physics';
132        varargin(i) = [];
133        if length(varargin) >= i
134            if isnumeric(varargin{i})
135                MCF = varargin{i};
136                varargin(i) = [];
137            end
138        end
139    elseif strcmpi(varargin{i},'hardware')
140        UnitsFlag = varargin{i};
141        varargin(i) = [];
142        if length(varargin) >= i
143            if isnumeric(varargin{i})
144                MCF = varargin{i};
145                varargin(i) = [];
146            end
147        end
148    elseif strcmpi(varargin{i},'simulator') | strcmpi(varargin{i},'model')
149        % Ignor
150        varargin(i) = [];
151    elseif strcmpi(varargin{i},'online')
152        % Ignor
153        varargin(i) = [];
154    elseif strcmpi(varargin{i},'manual')
155        % Ignor
156        varargin(i) = [];
157    end
158end
159
160
161% Look for DeltaRF input
162if length(varargin) >= 1
163    if isnumeric(varargin{1})
164        DeltaRF = varargin{1};
165        varargin(1) = [];
166    else
167        DeltaRF = [];
168    end
169else
170    DeltaRF = [];
171end
172
173
174% Look for BPMx family info
175if length(varargin) >= 1
176    if isstr(varargin{1})
177        Family1 = varargin{1};
178        varargin(1) = [];
179        if length(varargin) >= 1
180            if isnumeric(varargin{1})
181                DeviceList1 = varargin{1};
182                varargin(1) = [];
183            end
184        end
185    elseif isnumeric(varargin{1})
186        DeviceList1 = varargin{1};
187        varargin(1) = [];
188    elseif isstruct(varargin{1})
189        Family1 = varargin{1}.FamilyName;
190        DeviceList1 = varargin{1}.DeviceList;
191        if isempty(UnitsFlag)
192            UnitsFlag = varargin{1}.Units;
193        end
194        if ~NumericOutputFlag
195            % Only change StructOutputFlag if 'numeric' is not on the input line
196            StructOutputFlag = 1;
197        end
198        varargin(1) = [];     
199    end
200end
201
202% Look for BPMy family info
203if length(varargin) >= 1
204    if isstr(varargin{1})
205        Family2 = varargin{1};
206        varargin(1) = [];
207        if length(varargin) >= 1
208            if isnumeric(varargin{1})
209                DeviceList2 = varargin{1};
210                varargin(1) = [];
211            end
212        end
213    elseif isnumeric(varargin{1})
214        DeviceList2 = varargin{1};
215        varargin(1) = [];
216    elseif isstruct(varargin{1})
217        Family2 = varargin{1}.FamilyName;
218        DeviceList2 = varargin{1}.DeviceList;
219        if isempty(UnitsFlag)
220            UnitsFlag = varargin{1}.Units;
221        end
222        if ~NumericOutputFlag
223            % Only change StructOutputFlag if 'numeric' is not on the input line
224            StructOutputFlag = 1;
225        end
226        varargin(1) = [];     
227    end
228else
229    Family2 = Family1;
230    DeviceList2 = DeviceList1;
231end
232% End of input parsing
233
234
235% Archive data structure
236if ArchiveFlag
237    if isempty(FileName)
238        FileName = appendtimestamp(getfamilydata('Default', 'DispArchiveFile'));
239        DirectoryName = getfamilydata('Directory','DispData');
240        if isempty(DirectoryName)
241            DirectoryName = [getfamilydata('Directory','DataRoot') 'Dispersion', filesep];
242        else
243            % Make sure default directory exists
244            DirStart = pwd;
245            [DirectoryName, ErrorFlag] = gotodirectory(DirectoryName);
246            cd(DirStart);
247        end
248        [FileName, DirectoryName] = uiputfile('*.mat', 'Select Dispersion File', [DirectoryName FileName]);
249        if FileName == 0
250            ArchiveFlag = 0;
251            disp('   Dispersion measurement canceled.');
252            Dx=[]; Dy=[]; FileName='';
253            return
254        end
255        FileName = [DirectoryName, FileName];
256    elseif FileName == -1
257        FileName = appendtimestamp(getfamilydata('Default', 'DispArchiveFile'));
258        DirectoryName = getfamilydata('Directory','DispData');
259        if isempty(DirectoryName)
260            DirectoryName = [getfamilydata('Directory','DataRoot') 'Dispersion', filesep];
261        end
262        FileName = [DirectoryName, FileName];
263    end
264end
265
266
267% Get the input units
268if isempty(UnitsFlag)
269    UnitsFlag = 'Physics';
270%     if strcmpi(Family1, 'All') | strcmpi(Family2, 'All')
271%         UnitsFlag = 'Physics';
272%     else
273%         RFUnitsInput = getfamilydata('RF','Setpoint','Units');
274%         UnitsFlag = RFUnitsInput;
275%     end
276else
277    RFUnitsInput = UnitsFlag;
278end
279
280
281% Get DeltaRF in Hardware units
282RFHWUnits = getfamilydata('RF','Setpoint','HWUnits');
283
284
285% Make sure DeltaRF is in hardware units
286if isempty(DeltaRF) | ~isnumeric(DeltaRF)
287    % Get the default from the AD is in Hardware units
288    DeltaRF = getfamilydata('DeltaRFDisp');
289   
290    % If the default is not in the AD
291    if isempty(DeltaRF)
292        % Here is the second level default
293        DeltaRF = getrf('Model','Hardware') * modelmcf * .002;  % .2% energy change
294    end
295
296else
297    if strcmpi(UnitsFlag, 'Physics')
298        % Change to hardware
299        RFUnitsInput = UnitsFlag;
300        DeltaRF = physics2hw('RF', 'Setpoint', DeltaRF, [1 1]);
301    end
302end
303
304
305% DeltaRF are in hardware units at this point
306% Change to AT units
307DeltaRF = hw2physics('RF', 'Setpoint', DeltaRF, [1 1]);
308
309
310% Check DeltaRF for resonable values
311if DeltaRF > 15000;  % Hz
312    tmp = questdlg(sprintf('%f Hz is a large RF change.  Do you want to continue?', DeltaRF),'Dispersion Measurement','YES','NO','YES');
313    if strcmp(tmp,'NO')
314        Dx=[];  Dy=[];
315        return
316    end
317end
318
319if isempty(MCF)
320    MCF = modelmcf;
321end
322
323
324%%%%%%%%%%%%%%%%%%%%%%%%%%%%
325% Get the model dispersion %
326%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327
328% Find which orbit method to use
329[CavityState, PassMethod, iCavity] = getcavity;
330if isempty(CavityState)
331    % No cavity in AT model
332    if isradiationon
333        fprintf('   Turning radiation off since there is no RF cavity.\n');
334        setradiation off;
335    end
336    ATMethod = 'findorbit4';
337else
338    if strcmpi(deblank(CavityState(1,:)), 'On') | isradiationon
339        % findorbit6 recommended when cavity or radiation is on
340        if strcmpi(deblank(CavityState(1,:)), 'Off')
341            % Turn off cavity in AT model
342            fprintf('   Turning the RF cavity on, since radiation is on.\n');
343            setcavity('On');
344        end
345        ATMethod = 'findorbit6';
346    else
347        % Cavity is off in AT model
348        ATMethod = 'findsyncorbit';
349    end
350end
351
352
353if strcmpi(ATMethod, 'findsyncorbit')     
354    C = 2.99792458e8;
355    CavityFrequency  = THERING{iCavity(1)}.Frequency;
356    CavityHarmNumber = THERING{iCavity(1)}.HarmNumber;
357    L = findspos(THERING,length(THERING)+1)';  % getfamilydata('Circumference')
358    f0 = C * CavityHarmNumber / L;
359    RFChange = CavityFrequency - f0;   % To account for all the changes made to the RF [Hz]
360   
361    if strcmpi(ModulationMethod, 'bipolar')
362        ORBITPLUS = findsyncorbit(THERING, (-C*( DeltaRF+RFChange)*CavityHarmNumber/CavityFrequency^2)/2, 1:length(THERING)+1);
363        ORBIT0    = findsyncorbit(THERING, (-C*(-DeltaRF+RFChange)*CavityHarmNumber/CavityFrequency^2)/2, 1:length(THERING)+1);
364    else
365        ORBITPLUS = findsyncorbit(THERING, -C*(DeltaRF+RFChange)*CavityHarmNumber/CavityFrequency^2, 1:length(THERING)+1);
366        ORBIT0    = findsyncorbit(THERING, -C*(        RFChange)*CavityHarmNumber/CavityFrequency^2, 1:length(THERING)+1);
367    end   
368elseif strcmpi(ATMethod, 'findorbit4')
369    % Use findorbit4
370    % CavityFrequency is needed to for convert DeltaRF to dP
371    % Since on cavity is in the model, use the frequency calculated from the harmonic number
372    C = 2.99792458e8;
373    CavityHarmNumber = getfamilydata('HarmonicNumber');
374    if isempty(CavityHarmNumber)
375        error('Add HarmonicNumber to the ML so that the RF can be inferred when no cavity is in the model');
376    end
377    L = findspos(THERING,length(THERING)+1)';  % getfamilydata('Circumference')
378    CavityFrequency = C * CavityHarmNumber / L;
379    dP = DeltaRF / MCF/ CavityFrequency;
380    %DeltaRF = CavityFrequency * MCF * dP;
381    if strcmpi(ModulationMethod, 'bipolar')
382        ORBITPLUS = findorbit4(THERING, -dP/2, 1:length(THERING)+1);
383        ORBIT0    = findorbit4(THERING,  dP/2, 1:length(THERING)+1);
384    else
385        ORBITPLUS = findorbit4(THERING, -dP, 1:length(THERING)+1);
386        ORBIT0    = findorbit4(THERING,   0, 1:length(THERING)+1);
387    end   
388   
389elseif strcmpi(ATMethod, 'findorbit6')
390    % Use findorbit6
391    CavityFrequency = THERING{iCavity(1)}.Frequency;
392    if strcmpi(ModulationMethod, 'bipolar')
393        for kk = 1:length(iCavity)
394            THERING{iCavity(kk)}.Frequency = CavityFrequency + DeltaRF/2;
395        end               
396        ORBITPLUS = findorbit6(THERING, 1:length(THERING)+1);
397    else
398        for kk = 1:length(iCavity)
399            THERING{iCavity(kk)}.Frequency = CavityFrequency + DeltaRF;
400        end               
401        ORBITPLUS = findorbit6(THERING, 1:length(THERING)+1);
402    end   
403    if strcmpi(ModulationMethod, 'bipolar')
404        for kk = 1:length(iCavity)
405            THERING{iCavity(kk)}.Frequency = CavityFrequency - DeltaRF/2;
406        end               
407        ORBIT0 = findorbit6(THERING, 1:length(THERING)+1);
408        for kk = 1:length(iCavity)
409            THERING{iCavity(kk)}.Frequency = CavityFrequency;
410        end               
411    else
412        for kk = 1:length(iCavity)
413            THERING{iCavity(kk)}.Frequency = CavityFrequency;
414        end               
415        ORBIT0 = findorbit6(THERING, 1:length(THERING)+1);
416    end   
417end
418
419
420% Horizontal plane
421if strcmpi(Family1,'All')
422    Index1 = 1:length(THERING)+1;
423    if ~isempty(DeviceList1)
424        Index1 = Index1(DeviceList1);
425    end
426elseif isfamily(Family1)
427    Index1 = family2atindex(Family1, DeviceList1);
428else
429    Index1 = findcells(THERING, 'FamName', Family1);
430    if ~isempty(DeviceList1)
431        Index1 = Index1(DeviceList1);
432    end
433end
434if isempty(Index1)
435    error(sprintf('Family1=%s could not be found in the AO or AT deck',Family1));
436else
437    Index1 = Index1(:)';    % Row vector
438end
439
440D = ORBITPLUS([1 3], Index1) - ORBIT0([1 3], Index1);
441Sx = findspos(THERING, Index1);
442Sx = Sx(:)';
443Dx = D(1,:)' / DeltaRF;
444
445% Vertical plane
446if strcmpi(Family2,'All')
447    Index2 = 1:length(THERING)+1;
448    if ~isempty(DeviceList2)
449        Index2 = Index2(DeviceList2);
450    end
451elseif isfamily(Family2)
452    Index2 = family2atindex(Family2, DeviceList2);
453else
454    Index2 = findcells(THERING, 'FamName', Family2);
455    if ~isempty(DeviceList2)
456        Index2 = Index2(DeviceList2);
457    end
458end
459if isempty(Index2)
460    error(sprintf('Family2=%s could not be found in the AO or AT deck',Family2));
461else
462    Index2 = Index2(:)';    % Row vector
463end
464
465D = ORBITPLUS([1 3], Index2) - ORBIT0([1 3], Index2);
466Sy = findspos(THERING, Index2);
467Sy = Sy(:)';
468Dy = D(2,:)' / DeltaRF;
469
470% All elements
471D = ORBITPLUS([1 3], :) - ORBIT0([1 3], :);
472DxAll = D(1,:)' / DeltaRF;
473DyAll = D(2,:)' / DeltaRF;
474SAll = findspos(THERING, 1:length(THERING)+1);
475
476
477% Dispersion has units meters/Hz at this point
478% Change to true physics units meters/(dp/p)
479Dx = -CavityFrequency * MCF * Dx;
480Dy = -CavityFrequency * MCF * Dy;
481DxAll = -CavityFrequency * MCF * DxAll;
482DyAll = -CavityFrequency * MCF * DyAll;
483
484% Fill the dispersion structure (response matrix structure + some fields)
485d(1).Data = Dx;
486d(2).Data = Dy;
487d(1).FamilyName = 'DispersionX';
488d(2).FamilyName = 'DispersionY';
489
490d(1).Actuator = getsp('RF','Physics','Struct');
491d(2).Actuator = d(1).Actuator;
492% d(1).Actuator.Data = CavityFrequency;
493% d(2).Actuator.Data = CavityFrequency;
494% d(1).Actuator.Units = 'Physics';
495% d(2).Actuator.Units = 'Physics';
496% d(1).Actuator.UnitsString = 'Hz';
497% d(2).Actuator.UnitsString = 'Hz';
498d(1).ActuatorDelta = DeltaRF;
499d(2).ActuatorDelta = DeltaRF;
500
501if isfamily(Family1) & ismemberof(Family1, 'BPM')
502    d(1).Monitor = family2datastruct(Family1, 'Monitor', DeviceList1, 'Physics');
503else
504    d(1).Monitor.Data = ORBIT0(1, Index1)';
505    d(1).Monitor.FamilyName = ['x@', Family1];
506    d(1).Monitor.Field = 'Monitor';
507    d(1).Monitor.DeviceList = DeviceList1;
508    d(1).Monitor.Units = 'Physics';
509    d(1).Monitor.UnitsString = 'meters';
510    d(1).Monitor.TimeStamp = clock;
511    d(1).Monitor.t = 0;
512    d(1).Monitor.tout = 0;
513end
514if isfamily(Family2) & ismemberof(Family2, 'BPM')
515    d(2).Monitor = family2datastruct(Family2, 'Monitor', DeviceList2, 'Physics');
516else
517    d(2).Monitor.Data = ORBIT0(2, Index2)';
518    d(2).Monitor.FamilyName = ['y@', Family2];
519    d(2).Monitor.Field = 'Monitor';
520    d(2).Monitor.DeviceList = DeviceList2;
521    d(2).Monitor.Units = 'Physics';
522    d(2).Monitor.UnitsString = 'meters';
523    d(2).Monitor.TimeStamp = clock;
524    d(2).Monitor.t = 0;
525    d(2).Monitor.tout = 0;
526end
527
528d(1).TimeStamp = clock;
529d(2).TimeStamp = d(1).TimeStamp;
530
531for i=1:2
532    d(i).Mode = 'Model';
533    d(i).Units = 'Physics';
534    d(i).UnitsString = 'm/(dp/p)';
535    d(i).ModulationMethod = ModulationMethod;
536    d(i).OperationalMode = getfamilydata('OperationalMode');
537    d(i).DataDescriptor = 'Dispersion';
538    d(i).CreatedBy = 'modeldisp';
539    d(i).GeV = getenergy('Model');
540    d(i).WaitFlag = -2;
541    d(i).MCF = MCF;
542    d(i).dp = -DeltaRF / (CavityFrequency*MCF);
543end
544
545
546% Final units conversion
547if strcmpi(UnitsFlag, 'Hardware')
548    if ~ismemberof(Family1, 'BPM') | ~ismemberof(Family2, 'BPM')
549        error('Hardware units cannot be used when not using a BPM family');
550    end
551    d = physics2hw(d);
552end
553
554
555% Output
556% Plot if no output
557if DisplayFlag
558    % Plot dispersion
559   
560    % plotdisp cannot handle a AT family name or 'All'
561    % so a new dispersion plot was written below
562    %plotdisp(d);
563
564    UnitsFlag = d(1).Units;
565    UnitsString = d(1).UnitsString;
566
567    if strcmp(UnitsString, 'mm/MHz')
568        TitleString = sprintf('Change in Orbit / Change in RF (Delta RF=%g Hz)', DeltaRF);
569    elseif strcmp(UnitsString, 'meters/(dp/p)') | strcmp(UnitsString, 'm/(dp/p)')
570        TitleString = sprintf('Dispersion Function: %s  (\\alpha=%.2e, f_{RF}=%f Hz, \\Delta f_{RF}=%g Hz)', texlabel('-alpha f_{RF} {Delta}Orbit / {Delta}f_{RF}'), MCF, CavityFrequency, DeltaRF);
571    else
572        TitleString = 'Dispersion Function';
573    end
574   
575    h1 = subplot(5,1,[1 2]);
576    plot(SAll, DxAll, '-b');
577    if strcmpi(Family1,'All')
578        xlabel('Position [meters]');
579    else
580        hold on;
581        plot(Sx, Dx, '.b');
582        hold off;
583        xlabel(sprintf('%s Position [meters]', Family1));
584    end
585    xlim([0 SAll(end)]);
586    ylabel(sprintf('Horizontal [%s]', UnitsString));
587    title(TitleString, 'Interpreter', 'Tex');
588   
589    h2 = subplot(5,1,3);
590    drawlattice;
591   
592    h3 = subplot(5,1,[4 5]);
593    plot(SAll, DyAll, '-b');
594    if strcmpi(Family2,'All')
595        xlabel('Position [meters]');
596    else
597        hold on;
598        plot(Sy, Dy, '.b');
599        hold off;
600        xlabel(sprintf('%s Position [meters]', Family2));
601    end
602    ylabel(sprintf('Vertical [%s]', UnitsString));
603   
604    linkaxes([h1 h2 h3],'x')
605    set([h1 h2 h3],'XGrid','On','YGrid','On');
606end
607
608% Output
609if StructOutputFlag
610    Dx = d(1);
611    Dy = d(2);
612else
613    Dx = d(1).Data;
614    Dy = d(2).Data;
615end
616
617
618% Archive data structure
619if ArchiveFlag
620    % If the filename contains a directory then make sure it exists
621    [DirectoryName, FileName, Ext] = fileparts(FileName);
622    DirStart = pwd;
623    [DirectoryName, ErrorFlag] = gotodirectory(DirectoryName);
624    BPMxDisp = d(1);
625    BPMyDisp = d(2);
626    save(FileName, 'BPMxDisp', 'BPMyDisp');
627    if DisplayFlag
628        fprintf('   Model dispersion data saved to %s.mat\n', [DirectoryName FileName]);
629        if ErrorFlag
630            fprintf('   Warning: %s was not the desired directory\n', DirectoryName);
631        end
632    end
633    cd(DirStart);
634    FileName = [DirectoryName FileName];
635end
636if FileName == -1
637    FileName = '';
638end
639
640
Note: See TracBrowser for help on using the repository browser.