source: MML/trunk/mml/doc_html/mml/quadplotphysics.html @ 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: 40.5 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2                "http://www.w3.org/TR/REC-html40/loose.dtd">
3<html>
4<head>
5  <title>Description of quadplotphysics</title>
6  <meta name="keywords" content="quadplotphysics">
7  <meta name="description" content="QUADPLOTPHYSICS - Plots quadrupole centering data in physics units">
8  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
9  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
10  <meta name="robots" content="index, follow">
11  <link type="text/css" rel="stylesheet" href="../m2html.css">
12</head>
13<body>
14<a name="_top"></a>
15<div><a href="../index.html">Home</a> &gt;  <a href="index.html">mml</a> &gt; quadplotphysics.m</div>
16
17<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
18<td align="right"><a href="index.html">Index for mml&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->
19
20<h1>quadplotphysics
21</h1>
22
23<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
24<div class="box"><strong>QUADPLOTPHYSICS - Plots quadrupole centering data in physics units</strong></div>
25
26<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
27<div class="box"><strong>function [QMS, WarningString] = quadplotphysics(Input1, FigureHandle, sigmaBPM) </strong></div>
28
29<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
30<div class="fragment"><pre class="comment">QUADPLOTPHYSICS - Plots quadrupole centering data in physics units
31  [QMS, WarningString] = quadplot(Input1, Handle, sigmaBPM)
32
33  INPUTS
34  1. Input1 can be a filename for the data or a QMS structure (see help quadcenter for details).
35     If empty or zero inputs, then a dialog box will be provided to select a file.
36  2. Handle can be a figure handle or a vector of 4 axes handles
37     If Handle=0, no results are plotted
38  3. Standard deviation of the BPMs (scalar if all BPMs are the same)
39     These should be in the data file, but this provides an override if not
40     found then the default is inf (ie, not used).
41
42  OUTPUTS
43  1. For details of the QMS data structure see help quadcenter
44     This function added:
45     QMS.Offset - Offset computed at each BPM
46     QMS.FitParameters - Fit parameter at each BPM
47     QMS.FitParametersStd - Sigma of the fit parameter at each BPM
48     QMS.BPMStd - BPM sigma at each BPM
49     QMS.OffsetSTDMontiCarlo - Monti Carlo estimate of the sigma of the offset (optional)
50
51  2. WarningString = string with warning message if you occurred
52
53  Written by Greg Portmann</pre></div>
54
55<!-- crossreference -->
56<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
57This function calls:
58<ul style="list-style-image:url(../matlabicon.gif)">
59<li><a href="getcrunch.html" class="code" title="function Data = getcrunch(varargin)">getcrunch</a>     GETCRUNCH- Returns the crunch values for a family (radians)</li><li><a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>  GETFAMILYDATA - Gets data associated with the accelerator control</li><li><a href="getgain.html" class="code" title="function Data = getgain(varargin)">getgain</a>     GETGAIN - Returns the gain for a family</li><li><a href="getroll.html" class="code" title="function Data = getroll(varargin)">getroll</a>       GETROLL - Returns the roll values for a family (radians)</li><li><a href="quaderrors.html" class="code" title="function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag)">quaderrors</a>  QUADERRORS - Computes the error bars for quadrupole center measurement</li></ul>
60This function is called by:
61<ul style="list-style-image:url(../matlabicon.gif)">
62</ul>
63<!-- crossreference -->
64
65
66<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
67<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [QMS, WarningString] = quadplotphysics(Input1, FigureHandle, sigmaBPM)</a>
680002 <span class="comment">%QUADPLOTPHYSICS - Plots quadrupole centering data in physics units</span>
690003 <span class="comment">%  [QMS, WarningString] = quadplot(Input1, Handle, sigmaBPM)</span>
700004 <span class="comment">%</span>
710005 <span class="comment">%  INPUTS</span>
720006 <span class="comment">%  1. Input1 can be a filename for the data or a QMS structure (see help quadcenter for details).</span>
730007 <span class="comment">%     If empty or zero inputs, then a dialog box will be provided to select a file.</span>
740008 <span class="comment">%  2. Handle can be a figure handle or a vector of 4 axes handles</span>
750009 <span class="comment">%     If Handle=0, no results are plotted</span>
760010 <span class="comment">%  3. Standard deviation of the BPMs (scalar if all BPMs are the same)</span>
770011 <span class="comment">%     These should be in the data file, but this provides an override if not</span>
780012 <span class="comment">%     found then the default is inf (ie, not used).</span>
790013 <span class="comment">%</span>
800014 <span class="comment">%  OUTPUTS</span>
810015 <span class="comment">%  1. For details of the QMS data structure see help quadcenter</span>
820016 <span class="comment">%     This function added:</span>
830017 <span class="comment">%     QMS.Offset - Offset computed at each BPM</span>
840018 <span class="comment">%     QMS.FitParameters - Fit parameter at each BPM</span>
850019 <span class="comment">%     QMS.FitParametersStd - Sigma of the fit parameter at each BPM</span>
860020 <span class="comment">%     QMS.BPMStd - BPM sigma at each BPM</span>
870021 <span class="comment">%     QMS.OffsetSTDMontiCarlo - Monti Carlo estimate of the sigma of the offset (optional)</span>
880022 <span class="comment">%</span>
890023 <span class="comment">%  2. WarningString = string with warning message if you occurred</span>
900024 <span class="comment">%</span>
910025 <span class="comment">%  Written by Greg Portmann</span>
920026
930027
940028 <span class="comment">% To Do:</span>
950029 <span class="comment">% 1. It wouldn't be to difficult to added a LS weight based on slope or even the ideal weighting of std(center).</span>
960030 <span class="comment">%    I haven't done it yet because the BPM errors are usually roughly equal at most accelerators.</span>
970031
980032
990033 <span class="comment">% Remove BPM if it's slope less than MinSlopeFraction * (the maximum slope)</span>
1000034 MinSlopeFraction = .25;
1010035
1020036 <span class="comment">% # of STD of the center calculation allowed</span>
1030037 CenterOutlierFactor = 1;
1040038
1050039
1060040 QMS = [];
1070041 WarningString = <span class="string">''</span>;
1080042
1090043
1100044 <span class="comment">% Inputs</span>
1110045 <span class="keyword">try</span>
1120046     <span class="keyword">if</span> nargin == 0
1130047         [FileName, PathName] = uigetfile(<span class="string">'*.mat'</span>, <span class="string">'Select a Quadrupole Center File'</span>, [<a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Directory'</span>,<span class="string">'DataRoot'</span>), <span class="string">'QMS'</span>, filesep]);
1140048         <span class="keyword">if</span> ~isstr(FileName)
1150049             <span class="keyword">return</span>
1160050         <span class="keyword">else</span>
1170051             load([PathName,FileName]);
1180052         <span class="keyword">end</span>
1190053     <span class="keyword">else</span>
1200054         <span class="keyword">if</span> isempty(Input1)
1210055             [FileName, PathName] = uigetfile(<span class="string">'*.mat'</span>, <span class="string">'Select a Quadrupole Center File'</span>, [<a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Directory'</span>,<span class="string">'DataRoot'</span>), <span class="string">'QMS'</span>, filesep]);
1220056             <span class="keyword">if</span> ~isstr(FileName)
1230057                 <span class="keyword">return</span>
1240058             <span class="keyword">else</span>
1250059                 load([PathName,FileName]);
1260060             <span class="keyword">end</span>
1270061         <span class="keyword">elseif</span> isstr(Input1)
1280062             FileName = Input1;
1290063             load(FileName);
1300064         <span class="keyword">else</span>
1310065             QMS = Input1;
1320066             FileName = [];
1330067         <span class="keyword">end</span>
1340068     <span class="keyword">end</span>
1350069 <span class="keyword">catch</span>
1360070     error(<span class="string">'Problem getting input data'</span>);
1370071 <span class="keyword">end</span>
1380072 <span class="keyword">if</span> nargin &lt; 2
1390073     FigureHandle = [];
1400074 <span class="keyword">end</span>
1410075
1420076
1430077  QMS.QuadraticFit = 0;
1440078 
1450079 QuadraticFit = QMS.QuadraticFit;       <span class="comment">% 0 = linear fit, else quadratic fit</span>
1460080 OutlierFactor = QMS.OutlierFactor;     <span class="comment">% if abs(data - fit) &gt; OutlierFactor, then remove that BPM</span>
1470081
1480082 <span class="comment">% Get BPM standard deviation</span>
1490083 <span class="keyword">if</span> nargin &lt; 3
1500084     <span class="comment">% Get from the data file</span>
1510085     <span class="keyword">if</span> isfield(QMS, <span class="string">'BPMSTD'</span>)
1520086         sigmaBPM = QMS.BPMSTD;
1530087     <span class="keyword">else</span>
1540088         sigmaBPM = inf;
1550089     <span class="keyword">end</span>
1560090 <span class="keyword">end</span>
1570091 <span class="keyword">if</span> isempty(sigmaBPM)
1580092     sigmaBPM = inf;
1590093 <span class="keyword">end</span>
1600094 <span class="keyword">if</span> isnan(sigmaBPM) | isinf(sigmaBPM)
1610095     sigmaBPM = inf;
1620096     fprintf(<span class="string">'   WARNING: BPM standard deviation is unknown, hence there is no BPM outlier condition.\n'</span>);
1630097 <span class="keyword">end</span>
1640098 sigmaBPM = sigmaBPM(:);
1650099 QMS.BPMSTD = sigmaBPM;
1660100
1670101 <span class="comment">% Get figure handle</span>
1680102 <span class="keyword">if</span> all(FigureHandle ~= 0)
1690103     <span class="keyword">if</span> isempty(FigureHandle)
1700104         FigureHandle = figure;
1710105         clf reset
1720106         AxesHandles(1) = subplot(3,1,1);
1730107         AxesHandles(2) = subplot(3,1,2);
1740108         AxesHandles(3) = subplot(3,1,3);
1750109         <span class="comment">%AxesHandles(4) = subplot(4,1,4);</span>
1760110     <span class="keyword">else</span>
1770111         <span class="keyword">if</span> length(FigureHandle) == 1
1780112             FigureHandle = figure(FigureHandle);
1790113             clf reset
1800114             AxesHandles(1) = subplot(3,1,1);
1810115             AxesHandles(2) = subplot(3,1,2);
1820116             AxesHandles(3) = subplot(3,1,3);
1830117             <span class="comment">%AxesHandles(4) = subplot(4,1,4);</span>
1840118         <span class="keyword">elseif</span> length(FigureHandle) == 4
1850119             FigureHandle = figure;
1860120             clf reset
1870121             AxesHandles = FigureHandle;
1880122         <span class="keyword">else</span>
1890123             error(<span class="string">'Improper size of input FigureHandle'</span>);
1900124         <span class="keyword">end</span>
1910125     <span class="keyword">end</span>
1920126 <span class="keyword">end</span>
1930127
1940128 Buffer = .01;
1950129 HeightBuffer = .08;
1960130 <span class="keyword">if</span> QMS.QuadPlane == 1   
1970131     x1 = QMS.x1;
1980132     x2 = QMS.x2;
1990133     
2000134     <span class="comment">% Plot setup</span>
2010135     <span class="keyword">if</span> all(FigureHandle ~= 0)
2020136         set(FigureHandle,<span class="string">'units'</span>,<span class="string">'normal'</span>,<span class="string">'position'</span>,[.0+Buffer .27+Buffer .5-2*Buffer .72-2*Buffer-HeightBuffer]);
2030137     <span class="keyword">end</span>
2040138 <span class="keyword">else</span>
2050139     x1 = QMS.y1;
2060140     x2 = QMS.y2;
2070141     
2080142     <span class="comment">% Plot setup</span>
2090143     <span class="keyword">if</span> all(FigureHandle ~= 0)
2100144         set(FigureHandle,<span class="string">'units'</span>,<span class="string">'normal'</span>,<span class="string">'position'</span>,[.5+Buffer .27+Buffer .5-2*Buffer .72-2*Buffer-HeightBuffer]);
2110145     <span class="keyword">end</span>
2120146 <span class="keyword">end</span>
2130147
2140148
2150149 [BPMelem1, iNotFound] = findrowindex(QMS.BPMDev, QMS.BPMDevList);
2160150 <span class="keyword">if</span> ~isempty(iNotFound)
2170151     error(<span class="string">'BPM at the quadrupole not found in the BPM device list'</span>);
2180152 <span class="keyword">end</span>
2190153
2200154 <span class="comment">% Expand sigmaBPM is necessary</span>
2210155 <span class="keyword">if</span> length(sigmaBPM) == 1
2220156     sigmaBPM = ones(size(x1,1),1) * sigmaBPM;
2230157 <span class="keyword">end</span>
2240158
2250159 N = size(x1,2);
2260160
2270161 <span class="comment">% Change the number of points</span>
2280162 <span class="comment">% if 0</span>
2290163 <span class="comment">%     Ndiv2 = floor(size(x1,2)/2);</span>
2300164 <span class="comment">%     Npoint1 = Ndiv2;</span>
2310165 <span class="comment">%     Npoint2 = Ndiv2+2;</span>
2320166 <span class="comment">%     fprintf('  Using %d points (%d to %d, total %d).', Npoint2-Npoint1+1, Npoint1, Npoint2, N)</span>
2330167 <span class="comment">%     x1 = x1(:,Npoint1:Npoint2);</span>
2340168 <span class="comment">%     x2 = x2(:,Npoint1:Npoint2);</span>
2350169 <span class="comment">%     y1 = y1(:,Npoint1:Npoint2);</span>
2360170 <span class="comment">%     y2 = y2(:,Npoint1:Npoint2);</span>
2370171 <span class="comment">%     N = size(x1,2);</span>
2380172 <span class="comment">%     Ndiv2 = floor(size(x1,2)/2);</span>
2390173 <span class="comment">% end</span>
2400174
2410175
2420176
2430177
2440178 <span class="comment">%</span>
2450179 <span class="comment">% QUAD0 = getquad(QMS, 'Model');</span>
2460180 <span class="comment">% CM0 = getsp(QMS.CorrFamily, QMS.CorrDevList, 'Model');</span>
2470181 <span class="comment">%</span>
2480182 <span class="comment">%</span>
2490183 <span class="comment">% % Start the corrector a little lower first for hysteresis reasons</span>
2500184 <span class="comment">% CorrStep = 2 * QMS.CorrDelta / (N-1);</span>
2510185 <span class="comment">% stepsp(QMS.CorrFamily, -QMS.CorrDelta, QMS.CorrDevList, -1, 'Model');</span>
2520186 <span class="comment">%</span>
2530187 <span class="comment">% %XstartModel = getam(BPMxFamily, BPMxDev)</span>
2540188 <span class="comment">% for i = 1:N</span>
2550189 <span class="comment">%     % Step the vertical orbit</span>
2560190 <span class="comment">%     if i ~= 1</span>
2570191 <span class="comment">%         stepsp(QMS.CorrFamily, CorrStep, QMS.CorrDevList, -1, 'Model');</span>
2580192 <span class="comment">%     end</span>
2590193 <span class="comment">%</span>
2600194 <span class="comment">% %    fprintf('   %d. %s(%d,%d) = %+5.2f, %s(%d,%d) = %+.5f %s\n', i, QMS.CorrFamily, QMS.CorrDevList(1,:), getsp(QMS.CorrFamily, QMS.CorrDevList(1,:),'Model'),  BPMyFamily, BPMyDev, getam(BPMyFamily, BPMyDev,'Model'), QMS_Vertical.Orbit0.UnitsString); pause(0);</span>
2610195 <span class="comment">%</span>
2620196 <span class="comment">%     %OrbitCorrection(XstartModel, BPMxFamily, BPMxDev, HCMFamily, HCMDev, 2, 'Model');</span>
2630197 <span class="comment">%</span>
2640198 <span class="comment">%     if strcmpi(lower(QMS.ModulationMethod), 'sweep')</span>
2650199 <span class="comment">%         % One dimensional sweep of the quadrupole</span>
2660200 <span class="comment">%         xm1(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');</span>
2670201 <span class="comment">%         xm0(:,i) = xm1(:,i);</span>
2680202 <span class="comment">%         setquad(QMS, i*QMS.QuadDelta+QUAD0, -1, 'Model');</span>
2690203 <span class="comment">%         xm2(:,i) = getam(QMS.BPMFamily, QMS.BPMDevList, 'Model');</span>
2700204 <span class="comment">%     elseif strcmpi(lower(QMS.ModulationMethod), 'bipolar')</span>
2710205 <span class="comment">%         % Modulate the quadrupole</span>
2720206 <span class="comment">%         xm0(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');</span>
2730207 <span class="comment">%         [xq0(:,i), yq0(:,i)] = modeltwiss('x', QMS.QuadFamily, QMS.QuadDev);</span>
2740208 <span class="comment">%</span>
2750209 <span class="comment">%         setquad(QMS, QMS.QuadDelta+QUAD0, -1, 'Model');</span>
2760210 <span class="comment">%         xm1(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');</span>
2770211 <span class="comment">%         [xq1(:,i), yq1(:,i)] = modeltwiss('x', QMS.QuadFamily, QMS.QuadDev);</span>
2780212 <span class="comment">%</span>
2790213 <span class="comment">%</span>
2800214 <span class="comment">%         setquad(QMS,-QMS.QuadDelta+QUAD0, -1, 'Model');</span>
2810215 <span class="comment">%         xm2(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');</span>
2820216 <span class="comment">%         [xq2(:,i), yq2(:,i)] = modeltwiss('x', QMS.QuadFamily, QMS.QuadDev);</span>
2830217 <span class="comment">%</span>
2840218 <span class="comment">%         setquad(QMS, QUAD0, -1, 'Model');</span>
2850219 <span class="comment">%     elseif strcmpi(lower(QMS.ModulationMethod), 'unipolar')</span>
2860220 <span class="comment">%         % Modulate the quadrupole</span>
2870221 <span class="comment">%         xm1(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');</span>
2880222 <span class="comment">%         xm0(:,i) = x1(:,i);</span>
2890223 <span class="comment">%         setquad(QMS, QMS.QuadDelta+QUAD0, -1, 'Model');</span>
2900224 <span class="comment">%         xm2(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');</span>
2910225 <span class="comment">%         setquad(QMS, QUAD0, -1, 'Model');</span>
2920226 <span class="comment">%     end</span>
2930227 <span class="comment">% end</span>
2940228 <span class="comment">%</span>
2950229 <span class="comment">% setquad(QMS, QUAD0, -1, 'Model');</span>
2960230 <span class="comment">% setsp(QMS.CorrFamily, CM0, QMS.CorrDevList, -1, 'Model');</span>
2970231 <span class="comment">%</span>
2980232 <span class="comment">% xq0 = 1000*xq0;</span>
2990233 <span class="comment">% xq1 = 1000*xq1;</span>
3000234 <span class="comment">% xq2 = 1000*xq2;</span>
3010235 <span class="comment">%</span>
3020236 <span class="comment">% yq0 = 1000*yq0;</span>
3030237 <span class="comment">% yq1 = 1000*yq1;</span>
3040238 <span class="comment">% yq2 = 1000*yq2;</span>
3050239 <span class="comment">%</span>
3060240 <span class="comment">% x = xm0 + yq0-(xm0-xm0(3));</span>
3070241 <span class="comment">% x = xm0(3) + yq0-yq0(3);</span>
3080242 <span class="comment">% x = x';</span>
3090243
3100244
3110245 <span class="comment">% Convert to physics units</span>
3120246 <span class="comment">%x0 = hw2physics('BPMx', 'Monitor', x0, QMS.BPMDevList);</span>
3130247 <span class="comment">%x1 = hw2physics('BPMx', 'Monitor', x1, QMS.BPMDevList);</span>
3140248 <span class="comment">%x2 = hw2physics('BPMx', 'Monitor', x2, QMS.BPMDevList);</span>
3150249
3160250
3170251 Gx = <a href="getgain.html" class="code" title="function Data = getgain(varargin)">getgain</a>(<span class="string">'BPMx'</span>, QMS.BPMDevList);
3180252 Gy = <a href="getgain.html" class="code" title="function Data = getgain(varargin)">getgain</a>(<span class="string">'BPMy'</span>, QMS.BPMDevList);
3190253 C = <a href="getcrunch.html" class="code" title="function Data = getcrunch(varargin)">getcrunch</a>(<span class="string">'BPMx'</span>, QMS.BPMDevList);
3200254 R = <a href="getroll.html" class="code" title="function Data = getroll(varargin)">getroll</a>(<span class="string">'BPMx'</span>, QMS.BPMDevList);
3210255
3220256 <span class="keyword">for</span> i = 1:length(Gx)
3230257     m = [cos(R(i)) -sin(R(i)); sin(R(i)) cos(R(i))] * [1 C(i); C(i) 1] * [Gx(i) 0;0 Gy(i)] / sqrt(1-C(i)^2);
3240258
3250259     <span class="keyword">for</span> j = 1:size(QMS.x1,2)
3260260         <span class="keyword">if</span> j == 1
3270261             <span class="keyword">if</span> isfield(QMS, <span class="string">'x0'</span>)
3280262                 a = m * [QMS.x0(i,j); QMS.y0(i,j)];
3290263                 QMS.x0(i,j) = a(1);
3300264                 QMS.y0(i,j) = a(2);
3310265             <span class="keyword">end</span>
3320266         <span class="keyword">end</span>
3330267
3340268         a = m * [QMS.x1(i,j); QMS.y1(i,j)];
3350269         QMS.x1(i,j) = a(1);
3360270         QMS.y1(i,j) = a(2);
3370271         
3380272         a = m * [QMS.x2(i,j); QMS.y2(i,j)];
3390273         QMS.x2(i,j) = a(1);
3400274         QMS.y2(i,j) = a(2);
3410275     <span class="keyword">end</span>
3420276 <span class="keyword">end</span>
3430277
3440278 x1 = QMS.y1;
3450279 x2 = QMS.y2;
3460280
3470281
3480282
3490283 <span class="keyword">if</span> isfield(QMS, <span class="string">'OldCenter'</span>)
3500284  <span class="comment">%   QMS.OldCenter = hw2physics(QMS.BPMFamily, 'Monitor', QMS.OldCenter, QMS.BPMDev);</span>
3510285 <span class="keyword">end</span>
3520286
3530287 <span class="keyword">if</span> isfield(QMS, <span class="string">'Orbit0'</span>)
3540288   <span class="comment">%  QMS.Orbit0 = hw2physics(QMS.Orbit0);</span>
3550289     BPMUnitsString = QMS.Orbit0.UnitsString;
3560290 <span class="keyword">else</span>
3570291     BPMUnitsString = <span class="string">'m'</span>;
3580292 <span class="keyword">end</span>
3590293
3600294
3610295 <span class="comment">% Fit verses the position at the BPM next to the quadrupole</span>
3620296 <span class="keyword">if</span> strcmpi(QMS.ModulationMethod, <span class="string">'sweep'</span>)
3630297     <span class="comment">% One dimensional sweep of the quadrupole</span>
3640298     <span class="comment">%x = x1(BPMelem1,:)';</span>
3650299     x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2;     
3660300 <span class="keyword">elseif</span> strcmpi(QMS.ModulationMethod, <span class="string">'bipolar'</span>)
3670301     <span class="comment">% Modulation of the quadrupole</span>
3680302     x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2;     
3690303 <span class="keyword">elseif</span> strcmpi(QMS.ModulationMethod, <span class="string">'unipolar'</span>)
3700304     <span class="comment">% Unipolar modulation of the quadrupole</span>
3710305     <span class="comment">%x = x1(BPMelem1,:)';</span>
3720306     x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2;     
3730307 <span class="keyword">end</span>
3740308
3750309
3760310 <span class="comment">% Figure #1</span>
3770311 <span class="keyword">if</span> all(FigureHandle ~= 0)
3780312     axes(AxesHandles(1));
3790313     plot(x, x2-x1, <span class="string">'-x'</span>);
3800314     <span class="comment">%plot(linspace(-DelHCM,DelHCM,3), x2-x1);</span>
3810315     xlabel(sprintf(<span class="string">'%s(%d,%d), raw values [%s]'</span>, QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
3820316     ylabel(sprintf(<span class="string">'%s [%s]'</span>, QMS.BPMFamily, BPMUnitsString));
3830317     <span class="keyword">if</span> isempty(FileName)
3840318         title(sprintf(<span class="string">'Center for %s(%d,%d) %s(%d,%d)'</span>, QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), QMS.QuadFamily, QMS.QuadDev), <span class="string">'interpreter'</span>, <span class="string">'none'</span>);
3850319     <span class="keyword">else</span>
3860320         title(sprintf(<span class="string">'Center for %s(%d,%d) %s(%d,%d) (%s)'</span>, QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), QMS.QuadFamily, QMS.QuadDev, FileName), <span class="string">'interpreter'</span>, <span class="string">'none'</span>);
3870321     <span class="keyword">end</span>
3880322     grid on
3890323     axis tight
3900324 <span class="keyword">end</span>
3910325
3920326 <span class="keyword">if</span> isempty(FileName)
3930327     fprintf(<span class="string">'   Calculating the center of %s(%d,%d) using %s(%d,%d)\n'</span>, QMS.QuadFamily, QMS.QuadDev, QMS.BPMFamily, QMS.BPMDev);
3940328 <span class="keyword">else</span>
3950329     fprintf(<span class="string">'   Calculating the center of %s(%d,%d) using %s(%d,%d) (Data file: %s)\n'</span>, QMS.QuadFamily, QMS.QuadDev, QMS.BPMFamily, QMS.BPMDev, FileName);
3960330 <span class="keyword">end</span>
3970331 fprintf(<span class="string">'   Quadrupole modulation delta = %.3f amps, max. corrector step = %.3f amps\n'</span>, QMS.QuadDelta, QMS.CorrDelta);
3980332
3990333
4000334 <span class="comment">% Least squares fit</span>
4010335 merit = x2 - x1;
4020336 <span class="keyword">if</span> QuadraticFit
4030337     X = [ones(size(x)) x x.^2];
4040338     fprintf(<span class="string">'   %d point parabolic least squares fit\n'</span>, N);
4050339 <span class="keyword">else</span>
4060340     X = [ones(size(x)) x];
4070341     fprintf(<span class="string">'   %d point linear least squares fit\n'</span>, N);
4080342 <span class="keyword">end</span>
4090343
4100344
4110345 <span class="comment">% Axes #2</span>
4120346 <span class="keyword">if</span> all(FigureHandle ~= 0)
4130347     axes(AxesHandles(2));
4140348     xx = linspace(x(1), x(end), 200);
4150349 <span class="keyword">end</span>
4160350
4170351 iBPMOutlier = [];
4180352 invXX   = inv(X'*X);
4190353 invXX_X = invXX * X';
4200354
4210355 <span class="keyword">for</span> i = 1:size(x1,1)
4220356     <span class="comment">% least-square fit: m = slope and b = Y-intercept</span>
4230357     b = invXX_X * merit(i,:)';
4240358     bhat(i,:) = b';
4250359     
4260360     <span class="comment">% Should equal</span>
4270361     <span class="comment">%b = X \merit(i,:)';</span>
4280362     <span class="comment">%bhat1(i,:) = b';</span>
4290363     
4300364     <span class="comment">% Standard deviation</span>
4310365     bstd = sigmaBPM(i) * invXX;
4320366     bhatstd(i,:) = diag(bstd)';  <span class="comment">% hopefully cross-correlation terms are small</span>
4330367     
4340368     <span class="keyword">if</span> all(FigureHandle ~= 0)   
4350369         <span class="keyword">if</span> QuadraticFit
4360370             y = b(3)*xx.^2 + b(2)*xx + b(1);
4370371         <span class="keyword">else</span>
4380372             y = b(2)*xx + b(1);
4390373         <span class="keyword">end</span>
4400374 <span class="comment">%        plot(xx, y); hold on</span>
4410375     <span class="keyword">end</span>
4420376     
4430377     <span class="comment">% Outlier condition: remove if the error between the fit and the data is greater than OutlierFactor</span>
4440378     <span class="keyword">if</span> QuadraticFit
4450379         y = b(3)*x.^2 + b(2)*x + b(1);
4460380     <span class="keyword">else</span>
4470381         y = b(2)*x + b(1);
4480382     <span class="keyword">end</span>
4490383     <span class="keyword">if</span> max(abs(y - merit(i,:)')) &gt; OutlierFactor * sigmaBPM(i)    <span class="comment">% OutlierFactor was absolute max(abs(y - merit(i,:)')) &gt; OutlierFactor</span>
4500384         iBPMOutlier = [iBPMOutlier;i];
4510385     <span class="keyword">end</span>
4520386     
4530387     <span class="keyword">if</span> QuadraticFit
4540388         <span class="comment">% Quadratic fit</span>
4550389         <span class="keyword">if</span> b(2) &gt; 0
4560390             offset(i,1) = (-b(2) + sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3));
4570391         <span class="keyword">else</span>
4580392             offset(i,1) = (-b(2) - sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3));
4590393         <span class="keyword">end</span>
4600394         <span class="keyword">if</span> ~isreal(offset(i,1))
4610395             <span class="comment">% (b^2-4ac) can be negative but it will only happen if the slope is very small.  The offset</span>
4620396             <span class="comment">% should just get thrown out later as an outlier but change the solution to the minimum of the parabola.</span>
4630397             offset(i,1) = -b(2) / b(1) / 2;
4640398         <span class="keyword">end</span>
4650399     <span class="keyword">else</span>
4660400         <span class="comment">% Linear fit</span>
4670401         offset(i,1) = -b(1)/b(2);
4680402     <span class="keyword">end</span>
4690403 <span class="keyword">end</span>
4700404
4710405 QMS.Offset = offset;
4720406 QMS.FitParameters = bhat;
4730407 QMS.FitParametersStd = bhatstd;
4740408 QMS.BPMStd = sigmaBPM;
4750409
4760410
4770411 <span class="comment">% % Label axes #2</span>
4780412 <span class="comment">% if all(FigureHandle ~= 0)</span>
4790413 <span class="comment">%     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));</span>
4800414 <span class="comment">%     ylabel(sprintf('BPM LS Fit [%s]', BPMUnitsString));</span>
4810415 <span class="comment">%     grid on</span>
4820416 <span class="comment">%     axis tight</span>
4830417 <span class="comment">% end</span>
4840418
4850419
4860420 <span class="comment">% Compute offset for different conditions</span>
4870421 fprintf(<span class="string">'   %d total BPMs\n'</span>, length(offset));
4880422 fprintf(<span class="string">'   BPMs are removed for 1. Bad Status, 2. BPM Outlier, 3. Small Slope, or 4. Center Outlier\n'</span>);
4890423 <span class="comment">%m1 = mean(offset);</span>
4900424 <span class="comment">%s1 = std(offset);</span>
4910425 <span class="comment">%fprintf('   0. Mean = %.5f %s, STD = %.5f %s, all %d BPMs\n', m1, BPMUnitsString, s1, BPMUnitsString, length(offset));</span>
4920426
4930427
4940428 <span class="comment">% Remove bad Status BPMs</span>
4950429 iStatus = find(QMS.BPMStatus==0);
4960430 iBad = iStatus;
4970431 <span class="keyword">if</span> length(iBad) == length(offset)
4980432     error(<span class="string">'All the BPMs have bad status'</span>);
4990433 <span class="keyword">end</span>
5000434 offset1 = offset;
5010435 offset1(iBad) = [];
5020436 m2 = mean(offset1);
5030437 s2 = std(offset1);
5040438 fprintf(<span class="string">'   1. Mean = %+.5f %s, STD = %.5f %s, %2d points with bad status\n'</span>, m2, BPMUnitsString, s2, BPMUnitsString, length(iBad));
5050439
5060440 <span class="comment">% Remove bad Status + Outliers</span>
5070441 iBad = unique([iStatus; iBPMOutlier]);
5080442 <span class="keyword">if</span> length(iBad) == length(offset)
5090443     error(<span class="string">'All BPMs either have bad status or failed the BPM outlier condition'</span>);
5100444 <span class="keyword">end</span>
5110445 offset1 = offset;
5120446 offset1(iBad) = [];
5130447 m2 = mean(offset1);
5140448 s2 = std(offset1);
5150449 fprintf(<span class="string">'   2. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1 and abs(fit - measured data) &gt; %.2f std(BPM) (BPM outlier)\n'</span>, m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), OutlierFactor);
5160450
5170451
5180452 <span class="comment">% Remove bad Status + Small slopes</span>
5190453 <span class="comment">%iSlope = find(abs(bhat(:,2)) &lt; max(abs(bhat(:,2)))*MinSlopeFraction);</span>
5200454
5210455 <span class="comment">% Look for slope outliers</span>
5220456 Slopes = abs(bhat(:,2));
5230457 [Slopes, i] = sort(Slopes);
5240458 Slopes = Slopes(round(end/2):end);  <span class="comment">% remove the first half</span>
5250459 <span class="keyword">if</span> length(Slopes) &gt; 5
5260460     SlopesMax = Slopes(end-4);
5270461 <span class="keyword">else</span>
5280462     SlopesMax = Slopes(end);
5290463 <span class="keyword">end</span>
5300464 <span class="comment">%i = find(abs(Slopes-mean(Slopes)) &gt; 2 * std(Slopes));</span>
5310465 <span class="comment">%Slopes(i) = [];</span>
5320466
5330467 iSlope = find(abs(bhat(:,2)) &lt; SlopesMax * MinSlopeFraction);
5340468 iBad = unique([iStatus; iBPMOutlier; iSlope]);
5350469 <span class="keyword">if</span> length(iBad) == length(offset)
5360470     error(<span class="string">'All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition'</span>);
5370471 <span class="keyword">end</span>
5380472 offset1 = offset;
5390473 offset1(iBad) = [];
5400474 m2 = mean(offset1);
5410475 s2 = std(offset1);
5420476 fprintf(<span class="string">'   3. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1, 2, and slope &lt; %.2f max(slope)\n'</span>, m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), MinSlopeFraction);
5430477
5440478
5450479 <span class="comment">% Offset outlier offsets-mean(offsets) greater than 1 std</span>
5460480 itotal = (1:length(offset))';
5470481 iok = itotal;
5480482
5490483 offset1 = offset;
5500484 offset1(iBad) = [];
5510485 iok(iBad) = [];
5520486
5530487 i = find(abs(offset1-mean(offset1)) &gt; CenterOutlierFactor * std(offset1));
5540488 iCenterOutlier = iok(i);
5550489 iBad = unique([iBad; iCenterOutlier]);
5560490 <span class="keyword">if</span> length(iBad) == length(offset)
5570491     error(<span class="string">'All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition, , or failed the center outlier condition'</span>);
5580492 <span class="keyword">end</span>
5590493 offset1(i) = [];
5600494 iok(i) = [];
5610495
5620496 m2 = mean(offset1);
5630497 s2 = std(offset1);
5640498 QMS.GoodIndex = iok;
5650499 fprintf(<span class="string">'   4. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1, 2, 3, and abs(center-mean(center)) &gt; %.2f std(center) (Center outlier)\n'</span>, m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), CenterOutlierFactor);
5660500
5670501
5680502 NN = length(offset);
5690503
5700504 <span class="comment">% Axes #4</span>
5710505 <span class="keyword">if</span> all(FigureHandle ~= 0)
5720506     axes(AxesHandles(3));
5730507     [xx1,yy1]=stairs(1:NN,offset);
5740508     offset1 = offset;
5750509     offset1(iBad) = NaN*ones(length(iBad),1);
5760510     [xx2, yy2] = stairs(1:NN, offset1);
5770511     plot(xx1,yy1,<span class="string">'r'</span>, xx2,yy2,<span class="string">'b'</span>);
5780512     xlabel(<span class="string">'BPM Number'</span>);
5790513     ylabel(sprintf(<span class="string">'BPM Center [%s]'</span>, BPMUnitsString));
5800514     grid on
5810515     axis tight
5820516     <span class="comment">%xaxis([0 NN+1]);</span>
5830517     axis([0 NN+1 min(offset1)-.1e-6 max(offset1)+.1e-6]);
5840518 <span class="keyword">end</span>
5850519
5860520
5870521 <span class="comment">% Axes #2</span>
5880522
5890523 <span class="keyword">if</span> all(FigureHandle ~= 0)
5900524     <span class="keyword">if</span> 0
5910525         <span class="comment">% Plot red line over the bad lines</span>
5920526         axes(AxesHandles(2));
5930527         <span class="keyword">for</span> j = 1:length(iBad)
5940528             i = iBad(j);
5950529             <span class="keyword">if</span> QuadraticFit
5960530                 y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
5970531             <span class="keyword">else</span>
5980532                 y = bhat(i,2)*xx + bhat(i,1);
5990533             <span class="keyword">end</span>
6000534             plot(xx, y,<span class="string">'r'</span>);
6010535         <span class="keyword">end</span>
6020536         hold off
6030537         axis tight
6040538     <span class="keyword">else</span>
6050539         <span class="comment">% Only plot the good data</span>
6060540         axes(AxesHandles(2));
6070541         yy = [];
6080542         <span class="keyword">for</span> i = 1:size(x1,1)
6090543             <span class="keyword">if</span> ~any(i == iBad)
6100544                 <span class="keyword">if</span> QuadraticFit
6110545                     y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
6120546                 <span class="keyword">else</span>
6130547                     y = bhat(i,2)*xx + bhat(i,1);
6140548                 <span class="keyword">end</span>
6150549                 yy = [yy;y];
6160550             <span class="keyword">end</span>
6170551         <span class="keyword">end</span>
6180552         plot(xx, yy);
6190553         hold off
6200554         grid on
6210555         axis tight
6220556     <span class="keyword">end</span>
6230557     xlabel(sprintf(<span class="string">'%s(%d,%d), raw values [%s]'</span>, QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
6240558     ylabel(sprintf(<span class="string">'BPM LS Fit [%s]'</span>, BPMUnitsString));
6250559     grid on
6260560     axis tight
6270561 <span class="keyword">end</span>
6280562
6290563
6300564
6310565 <span class="keyword">if</span> ~isempty(iStatus)
6320566     <span class="keyword">if</span> ~isempty(find(iStatus==BPMelem1))
6330567         fprintf(<span class="string">'   WARNING: BPM(%d,%d) has a bad status\n'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6340568         WarningString = sprintf(<span class="string">'BPM(%d,%d) has a bad status'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6350569     <span class="keyword">end</span>
6360570 <span class="keyword">end</span>
6370571 <span class="keyword">if</span> ~isempty(iBPMOutlier)
6380572     <span class="keyword">if</span> ~isempty(find(iBPMOutlier==BPMelem1))
6390573         fprintf(<span class="string">'   WARNING: BPM(%d,%d) removed due to outlier (based on std(BPM))\n'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6400574         WarningString = sprintf(<span class="string">'BPM(%d,%d) removed due to outlier (based on std(BPM))'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6410575     <span class="keyword">end</span>
6420576 <span class="keyword">end</span>
6430577 <span class="keyword">if</span> ~isempty(iSlope)
6440578     <span class="keyword">if</span> ~isempty(find(iSlope==BPMelem1))
6450579         fprintf(<span class="string">'   WARNING: BPM(%d,%d) slope is too small\n'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6460580         WarningString = sprintf(<span class="string">'BPM(%d,%d) slope is too small'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6470581     <span class="keyword">end</span>
6480582 <span class="keyword">end</span>
6490583 <span class="keyword">if</span> ~isempty(iCenterOutlier)
6500584     <span class="keyword">if</span> ~isempty(find(iCenterOutlier==BPMelem1))
6510585         fprintf(<span class="string">'   WARNING: BPM(%d,%d) removed due to outlier (based on all the centers)\n'</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6520586         WarningString = sprintf(<span class="string">'BPM(%d,%d) '</span>, QMS.BPMDev(1), QMS.BPMDev(2));
6530587     <span class="keyword">end</span>
6540588 <span class="keyword">end</span>
6550589
6560590
6570591 <span class="comment">% % Axes #3</span>
6580592 <span class="comment">% if all(FigureHandle ~= 0)</span>
6590593 <span class="comment">%     axes(AxesHandles(3));</span>
6600594 <span class="comment">%     iii=1:NN;</span>
6610595 <span class="comment">%     iii(iBad)=[];</span>
6620596 <span class="comment">%     for j = 1:length(iii)</span>
6630597 <span class="comment">%         i = iii(j);</span>
6640598 <span class="comment">%</span>
6650599 <span class="comment">%         if 1</span>
6660600 <span class="comment">%             % Plot fit</span>
6670601 <span class="comment">%             if QuadraticFit</span>
6680602 <span class="comment">%                 y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);</span>
6690603 <span class="comment">%             else</span>
6700604 <span class="comment">%                 y = bhat(i,2)*xx + bhat(i,1);</span>
6710605 <span class="comment">%             end</span>
6720606 <span class="comment">%             if all(FigureHandle ~= 0)</span>
6730607 <span class="comment">%                 plot(xx, y,'b'); hold on</span>
6740608 <span class="comment">%             end</span>
6750609 <span class="comment">%         else</span>
6760610 <span class="comment">%             % Plot error in fit</span>
6770611 <span class="comment">%             if QuadraticFit</span>
6780612 <span class="comment">%                 y = bhat(i,3)*x.^2 + bhat(i,2)*x + bhat(i,1);</span>
6790613 <span class="comment">%             else</span>
6800614 <span class="comment">%                 y = bhat(i,2)*x + bhat(i,1);</span>
6810615 <span class="comment">%             end</span>
6820616 <span class="comment">%             plot(x, y - merit(i,:)','b'); hold on</span>
6830617 <span class="comment">%         end</span>
6840618 <span class="comment">%     end</span>
6850619 <span class="comment">%     hold off</span>
6860620 <span class="comment">%     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));</span>
6870621 <span class="comment">%     ylabel(sprintf('Final %s Merit Fun [%s]', QMS.BPMFamily, BPMUnitsString));</span>
6880622 <span class="comment">%     grid on</span>
6890623 <span class="comment">%     axis tight</span>
6900624 <span class="comment">%     orient tall</span>
6910625 <span class="comment">% end</span>
6920626
6930627 <span class="keyword">if</span> ~isempty(QMS.OldCenter)
6940628     fprintf(<span class="string">'   Starting Offset %s(%d,%d) = %+f [%s]\n'</span>, QMS.BPMFamily, QMS.BPMDev, QMS.OldCenter, BPMUnitsString);
6950629 <span class="keyword">end</span>
6960630 fprintf(<span class="string">'   New Offset      %s(%d,%d) = %+f [%s]\n'</span>, QMS.BPMFamily, QMS.BPMDev, m2, BPMUnitsString);
6970631
6980632 QMS.Center = m2;
6990633 QMS.CenterSTD = s2;
7000634
7010635
7020636 <span class="keyword">if</span> all(FigureHandle ~= 0)
7030637     addlabel(1, 0, datestr(QMS.TimeStamp));
7040638     addlabel(0, 0, sprintf(<span class="string">'Offset %s(%d,%d)=%f, {\\Delta}%s(%d,%d)=%f, {\\Delta}%s(%d,%d)=%f'</span>, QMS.BPMFamily, QMS.BPMDev, QMS.Center, QMS.QuadFamily, QMS.QuadDev, QMS.QuadDelta, QMS.CorrFamily, QMS.CorrDevList(1,:), QMS.CorrDelta(1)));
7050639 <span class="keyword">end</span>
7060640
7070641 fprintf(<span class="string">'\n'</span>);
7080642
7090643
7100644 <span class="comment">% Get and plot errors</span>
7110645 <span class="keyword">if</span> all(FigureHandle ~= 0)
7120646     QMS = <a href="quaderrors.html" class="code" title="function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag)">quaderrors</a>(QMS, gcf+1);
7130647 <span class="keyword">else</span>
7140648     <span class="comment">% This is a cluge for now</span>
7150649     QMS = <a href="quaderrors.html" class="code" title="function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag)">quaderrors</a>(QMS, 0);
7160650 <span class="keyword">end</span>
7170651
7180652
7190653 <span class="comment">%QMS = orderfields(QMS);</span></pre></div>
720<hr><address>Generated on Mon 21-May-2007 15:29:18 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
721</body>
722</html>
Note: See TracBrowser for help on using the repository browser.