source: MML/trunk/mml/doc_html/mml/orbitcorrectionmethods.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: 55.3 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 orbitcorrectionmethods</title>
6  <meta name="keywords" content="orbitcorrectionmethods">
7  <meta name="description" content="ORBITCORRECTIONMETHODS - Some the orbit correction methods used on light sources">
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; orbitcorrectionmethods.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>orbitcorrectionmethods
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>ORBITCORRECTIONMETHODS - Some the orbit correction methods used on light sources</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 [OCS, SmatNoWeights, S, U, V] = orbitcorrectionmethods(OCS, Smat, S, U, V) </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">ORBITCORRECTIONMETHODS - Some the orbit correction methods used on light sources
31
32  [OCS, Smat, S, U, V] = orbitcorrectionmethods(OCS)                 % Get response matrix &amp; compute SVD correction
33  [OCS, Smat, S, U, V] = orbitcorrectionmethods(OCS, Smat)           % Compute SVD &amp; correction
34  [OCS, Smat, S, U, V] = orbitcorrectionmethods(OCS, Smat, S, U, V)  % Compute correction
35
36  OCS.FitRF - Determines which RF freqency method to use
37  See setorbit for information on all the different flags.
38 
39  NOTES
40  1. OCS.CM.Data is not changed by this function
41     OCS.CM.Delta is the delta correction
42  2. Both row and column weighting of the response matrix can be done.
43     The default is unity row weights and column weight by the std of each column.
44
45  See also <a href="setorbit.html" class="code" title="function [OCS, OCS0, V, S, ErrorFlag] = setorbit(varargin)">setorbit</a>, <a href="setorbitbump.html" class="code" title="function [OCS, OCS0, V, S, ErrorFlag] = setorbitbump(varargin)">setorbitbump</a>, <a href="rmdisp.html" class="code" title="function [DeltaRF, BPM, c, DispOrbit] = rmdisp(varargin)">rmdisp</a>, <a href="plotcm.html" class="code" title="function [DeltaRF, HCMEnergyChangeTotal, DeltaL] = plotcm(varargin)">plotcm</a></pre></div>
46
47<!-- crossreference -->
48<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
49This function calls:
50<ul style="list-style-image:url(../matlabicon.gif)">
51<li><a href="getdisp.html" class="code" title="function [Data, FileName] = getdisp(varargin)">getdisp</a>       GETDISP - Returns the dispersion for a family (from file)</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="gethbpmfamily.html" class="code" title="function Family = gethbpmfamily">gethbpmfamily</a>   GETHBPMFAMILY - Return the default horizontal BPM family</li><li><a href="gethcmfamily.html" class="code" title="function Family = gethcmfamily">gethcmfamily</a>       GETHCMFAMILY - Returns the default horizontal corrector family</li><li><a href="getmcf.html" class="code" title="function Alpha = getmcf(ModelString)">getmcf</a>       GETMCF - Returns the momentum compaction factor (MCF) stored in the AD or the model</li><li><a href="getrespmat.html" class="code" title="function [S, FileName] = getrespmat(varargin)">getrespmat</a> GETRESPMAT - Get response matrix data from a file</li><li><a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>        GETSP - Gets setpoint channels</li><li><a href="getvbpmfamily.html" class="code" title="function Family = getvbpmfamily">getvbpmfamily</a>      GETVBPMFAMILY - Return the default vertical BPM family</li><li><a href="getvcmfamily.html" class="code" title="function Family = getvcmfamily">getvcmfamily</a> GETVCMFAMILY - Returns the default vertical corrector family</li><li><a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>    HW2PHYSICS - Converts from 'Hardware' units to 'Physics' units</li><li><a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>    ISMEMBEROF - Returns turn if the membership information of a family (cell of strings)</li><li><a href="measbpmresp.html" class="code" title="function [Rmat, OutputFileName] = measbpmresp(varargin)">measbpmresp</a>   MEASBPMRESP - Measures the BPM response matrix in the horizontal and vertical planes</li><li><a href="measdisp.html" class="code" title="function [Dx, Dy, FileName] = measdisp(varargin);">measdisp</a>        MEASDISP - Measures the dispersion function</li><li><a href="physics2hw.html" class="code" title="function S = physics2hw(Family, Field, value, DeviceList, Energy)">physics2hw</a>     PHYSICS2HW - Converts from 'Physics' units to 'Hardware' units</li></ul>
52This function is called by:
53<ul style="list-style-image:url(../matlabicon.gif)">
54<li><a href="setorbit.html" class="code" title="function [OCS, OCS0, V, S, ErrorFlag] = setorbit(varargin)">setorbit</a>        SETORBIT - Orbit correction function</li></ul>
55<!-- crossreference -->
56
57
58<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
59<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [OCS, SmatNoWeights, S, U, V] = orbitcorrectionmethods(OCS, Smat, S, U, V)</a>
600002 <span class="comment">%ORBITCORRECTIONMETHODS - Some the orbit correction methods used on light sources</span>
610003 <span class="comment">%</span>
620004 <span class="comment">%  [OCS, Smat, S, U, V] = orbitcorrectionmethods(OCS)                 % Get response matrix &amp; compute SVD correction</span>
630005 <span class="comment">%  [OCS, Smat, S, U, V] = orbitcorrectionmethods(OCS, Smat)           % Compute SVD &amp; correction</span>
640006 <span class="comment">%  [OCS, Smat, S, U, V] = orbitcorrectionmethods(OCS, Smat, S, U, V)  % Compute correction</span>
650007 <span class="comment">%</span>
660008 <span class="comment">%  OCS.FitRF - Determines which RF freqency method to use</span>
670009 <span class="comment">%  See setorbit for information on all the different flags.</span>
680010 <span class="comment">%</span>
690011 <span class="comment">%  NOTES</span>
700012 <span class="comment">%  1. OCS.CM.Data is not changed by this function</span>
710013 <span class="comment">%     OCS.CM.Delta is the delta correction</span>
720014 <span class="comment">%  2. Both row and column weighting of the response matrix can be done.</span>
730015 <span class="comment">%     The default is unity row weights and column weight by the std of each column.</span>
740016 <span class="comment">%</span>
750017 <span class="comment">%  See also setorbit, setorbitbump, rmdisp, plotcm</span>
760018
770019 <span class="comment">%</span>
780020 <span class="comment">%  Written by Gregory J. Portmann</span>
790021 <span class="comment">%  Addition for Soleil by Laurent S. Nadolski</span>
800022 <span class="comment">%    Constraints on min and max RF frequency variations</span>
810023
820024 <span class="comment">%  T = V(:,OCS.SVDIndex) * diag(S(OCS.SVDIndex).^(-1)) * U(:,OCS.SVDIndex)';</span>
830025 <span class="comment">%  Delcm = T * (BPMWeight .* (GoalOrbitVec - StartOrbitVec));</span>
840026 <span class="comment">%</span>
850027
860028
870029 <span class="comment">% Initialize</span>
880030 FitRFDefault = 0;
890031
900032 <span class="comment">% Machine dependent</span>
910033 <span class="keyword">if</span> strcmpi(<a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'SubMachine'</span>),<span class="string">'Soleil'</span>)
920034     DeltaRFmin = 1e-7;   <span class="comment">% MHz minimum variation allowed by Master oscillator</span>
930035     DeltaRFmax = 100e-6; <span class="comment">% MHz maximum variation allowed in one step</span>
940036 <span class="keyword">end</span>
950037
960038
970039 <span class="comment">% Input checking (should add more)</span>
980040 <span class="keyword">if</span> nargin &lt; 2
990041     Smat = [];
1000042 <span class="keyword">end</span>
1010043 <span class="keyword">if</span> nargin &lt; 3
1020044     S = [];
1030045     U = [];
1040046     V = [];
1050047 <span class="keyword">end</span>
1060048
1070049
1080050 <span class="comment">% This function expects .BPM and .CM to be cell arrays.</span>
1090051 <span class="comment">% This gets put back at the end of the function.</span>
1100052 <span class="keyword">if</span> ~iscell(OCS.BPM)
1110053     NoCellBPMFlag = 1;
1120054     OCS.BPM = {OCS.BPM};
1130055     OCS.GoalOrbit = {OCS.GoalOrbit};
1140056     <span class="keyword">if</span> isfield(OCS, <span class="string">'BPMWeight'</span>)
1150057         OCS.BPMWeight = {OCS.BPMWeight};
1160058     <span class="keyword">end</span>
1170059 <span class="keyword">else</span>
1180060     NoCellBPMFlag = 0;
1190061 <span class="keyword">end</span>
1200062 <span class="keyword">if</span> ~iscell(OCS.CM)
1210063     NoCellCMFlag = 1;
1220064     OCS.CM = {OCS.CM};
1230065     <span class="keyword">if</span> isfield(OCS, <span class="string">'CMWeight'</span>)
1240066         <span class="keyword">if</span> ~iscell(OCS.CMWeight)
1250067             OCS.CMWeight = {OCS.CMWeight};
1260068         <span class="keyword">end</span>
1270069     <span class="keyword">end</span>
1280070 <span class="keyword">else</span>
1290071     NoCellCMFlag = 0;
1300072 <span class="keyword">end</span>
1310073
1320074
1330075 <span class="comment">% RF frequency methods in orbit correction</span>
1340076 <span class="keyword">if</span> ~isfield(OCS, <span class="string">'FitRF'</span>)
1350077     OCS.FitRF = FitRFDefault;
1360078 <span class="keyword">end</span>
1370079 <span class="keyword">if</span> isempty(OCS.FitRF)
1380080     OCS.FitRF = FitRFDefault;
1390081 <span class="keyword">end</span>
1400082
1410083
1420084 <span class="comment">%%%%%%%%%%%%%%%%</span>
1430085 <span class="comment">% The SVD Part %</span>
1440086 <span class="comment">%%%%%%%%%%%%%%%%</span>
1450087 <span class="keyword">if</span> isempty(U)
1460088     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1470089     <span class="comment">% Get the response matrix %</span>
1480090     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1490091     <span class="keyword">if</span> isempty(Smat)
1500092         <span class="keyword">for</span> i = 1:length(OCS.BPM)
1510093             Srow = [];
1520094             <span class="keyword">for</span> j = 1:length(OCS.CM)
1530095                 <span class="comment">%if ModelRespFlag == 1</span>
1540096                 <span class="keyword">if</span> any(strcmpi(OCS.Flags,<span class="string">'ModelResp'</span>))
1550097                     <span class="comment">% When using measbpmresp('Model')</span>
1560098                     <span class="comment">% Rmat(1,1) is always horizontal</span>
1570099                     <span class="comment">% Rmat(2,2) is always vertical</span>
1580100                     Rmat = <a href="measbpmresp.html" class="code" title="function [Rmat, OutputFileName] = measbpmresp(varargin)">measbpmresp</a>(<span class="string">'model'</span>,OCS.BPM{i},OCS.BPM{i}, OCS.CM{j},OCS.CM{j}, OCS.BPM{i}.Units);
1590101                     <span class="keyword">if</span> <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.BPM{i}.FamilyName, <a href="gethbpmfamily.html" class="code" title="function Family = gethbpmfamily">gethbpmfamily</a>) &amp;&amp; <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.CM{j}.FamilyName,<a href="gethcmfamily.html" class="code" title="function Family = gethcmfamily">gethcmfamily</a>)
1600102                         S = Rmat(1,1).Data;
1610103                     <span class="keyword">elseif</span> <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.BPM{i}.FamilyName, <a href="getvbpmfamily.html" class="code" title="function Family = getvbpmfamily">getvbpmfamily</a>) &amp;&amp; <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.CM{j}.FamilyName,<a href="gethcmfamily.html" class="code" title="function Family = gethcmfamily">gethcmfamily</a>)
1620104                         S = Rmat(2,1).Data;               
1630105                     <span class="keyword">elseif</span> <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.BPM{i}.FamilyName, <a href="gethbpmfamily.html" class="code" title="function Family = gethbpmfamily">gethbpmfamily</a>) &amp;&amp; <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.CM{j}.FamilyName,<a href="getvcmfamily.html" class="code" title="function Family = getvcmfamily">getvcmfamily</a>)
1640106                         S = Rmat(1,2).Data;
1650107                     <span class="keyword">elseif</span> <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.BPM{i}.FamilyName, <a href="getvbpmfamily.html" class="code" title="function Family = getvbpmfamily">getvbpmfamily</a>) &amp;&amp; <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.CM{j}.FamilyName,<a href="getvcmfamily.html" class="code" title="function Family = getvcmfamily">getvcmfamily</a>)
1660108                         S = Rmat(2,2).Data;
1670109                     <span class="keyword">else</span>
1680110                         error(<span class="string">'Problem getting the model response matrix'</span>);
1690111                     <span class="keyword">end</span>
1700112                 <span class="keyword">else</span>
1710113                     <span class="comment">% Get the golden BPM response matrix</span>
1720114                     [S, FileName] = <a href="getrespmat.html" class="code" title="function [S, FileName] = getrespmat(varargin)">getrespmat</a>(OCS.BPM{i}, OCS.CM{j}, <span class="string">'Numeric'</span>);
1730115                     <span class="keyword">if</span> any(find(isnan(S)))
1740116                         error(<span class="string">'Not all BPMs and correctors exist in the response matrix.'</span>);
1750117                     <span class="keyword">end</span>
1760118                 <span class="keyword">end</span>
1770119                 <span class="keyword">if</span> length(OCS.GoalOrbit{i}) ~= size(S,1)
1780120                     error(<span class="string">'The goal orbit vector length does not equal the response matrix row length'</span>);
1790121                 <span class="keyword">end</span>
1800122
1810123                 Srow = [Srow S];
1820124             <span class="keyword">end</span>
1830125             Smat = [Smat; Srow];
1840126         <span class="keyword">end</span>
1850127         <span class="comment">%SmatCheck = measbpmresp('Model','Numeric');</span>
1860128
1870129
1880130         <span class="comment">%%%%%%%%%%%%%%%%%%</span>
1890131         <span class="comment">% Get Dispersion %</span>
1900132         <span class="comment">%%%%%%%%%%%%%%%%%%</span>
1910133         Eta = [];
1920134         <span class="keyword">if</span> OCS.FitRF
1930135             <span class="keyword">for</span> i = 1:length(OCS.BPM)
1940136                 <span class="keyword">if</span> any(strcmpi(OCS.Flags,<span class="string">'ModelDisp'</span>))
1950137                     <span class="comment">%if strcmpi(DispFlag,'ModelDisp')</span>
1960138                     DispOrbit = <a href="measdisp.html" class="code" title="function [Dx, Dy, FileName] = measdisp(varargin);">measdisp</a>(OCS.BPM{i}, <span class="string">'Hardware'</span>, <span class="string">'Numeric'</span>, <span class="string">'Model'</span>);
1970139                 <span class="keyword">elseif</span> any(strcmpi(OCS.Flags,<span class="string">'MeasDisp'</span>))
1980140                     <span class="comment">%elseif strcmpi(DispFlag,'MeasDisp')</span>
1990141                     DispOrbit = <a href="measdisp.html" class="code" title="function [Dx, Dy, FileName] = measdisp(varargin);">measdisp</a>(OCS.BPM{i}, <span class="string">'Hardware'</span>, <span class="string">'Numeric'</span>);
2000142                 <span class="keyword">elseif</span> any(strcmpi(OCS.Flags,<span class="string">'GoldenDisp'</span>))
2010143                     <span class="comment">%elseif strcmpi(DispFlag,'GoldenDisp')</span>
2020144                     DispOrbit = <a href="getdisp.html" class="code" title="function [Data, FileName] = getdisp(varargin)">getdisp</a>(OCS.BPM{i}, <span class="string">'Hardware'</span>, <span class="string">'Numeric'</span>);
2030145                 <span class="keyword">end</span>
2040146                 
2050147                 <span class="keyword">if</span> strcmpi(OCS.BPM{i}.Units, <span class="string">'Physics'</span>)
2060148                     <span class="comment">% Put the RF change in physics units, not energy change</span>
2070149                     RFhw = <a href="physics2hw.html" class="code" title="function S = physics2hw(Family, Field, value, DeviceList, Energy)">physics2hw</a>(OCS.RF);
2080150                     DispOrbit = DispOrbit / (OCS.RF.Data / RFhw.Data);
2090151                 <span class="keyword">end</span>
2100152                 
2110153                 <span class="keyword">if</span> isempty(DispOrbit)
2120154                     error(<span class="string">'Did not find or generate dispersion orbit properly'</span>);
2130155                 <span class="keyword">end</span>
2140156                 <span class="keyword">if</span> any(isnan(DispOrbit))
2150157                     error(<span class="string">'The dispersion data has at least one NaN.'</span>);
2160158                 <span class="keyword">end</span>
2170159                 
2180160                 Eta = [Eta; DispOrbit];
2190161             <span class="keyword">end</span>
2200162         <span class="keyword">end</span>
2210163         OCS.Eta = Eta;
2220164     <span class="keyword">end</span>
2230165 <span class="keyword">end</span>
2240166
2250167
2260168 <span class="comment">% Save a weight free response matrix</span>
2270169 SmatNoWeights = Smat;
2280170
2290171
2300172 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
2310173 <span class="comment">% Build the orbit vectors %</span>
2320174 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
2330175 StartOrbitVec = [];
2340176 GoalOrbitVec = [];
2350177 BPMWeight = [];
2360178 <span class="keyword">for</span> i = 1:length(OCS.BPM)
2370179     StartOrbitVec = [StartOrbitVec; OCS.BPM{i}.Data(:)];
2380180     GoalOrbitVec  = [GoalOrbitVec; OCS.GoalOrbit{i}(:)];
2390181
2400182     <span class="keyword">if</span> ~isfield(OCS, <span class="string">'BPMWeight'</span>) || isempty(OCS.BPMWeight{i})
2410183         BPMWeight = [BPMWeight; ones(length(OCS.BPM{i}.Data),1)];
2420184     <span class="keyword">elseif</span> isscalar(OCS.BPMWeight{i})
2430185         BPMWeight = [BPMWeight; ones(length(OCS.BPM{i}.Data),1) * OCS.BPMWeight{i}];
2440186     <span class="keyword">else</span>
2450187         BPMWeight = [BPMWeight; OCS.BPMWeight{i}(:)];
2460188     <span class="keyword">end</span>
2470189 <span class="keyword">end</span>
2480190
2490191
2500192 <span class="comment">% Build column weight vector</span>
2510193 CMWeight = [];
2520194 <span class="keyword">for</span> i = 1:length(OCS.CM)
2530195     <span class="keyword">if</span> ~isfield(OCS, <span class="string">'CMWeight'</span>) || isempty(OCS.CMWeight{i})
2540196         <span class="comment">% When using all singular values this does not do anything</span>
2550197         <span class="comment">% The amp to radian conversion is probably a better normalizer</span>
2560198         <span class="comment">%CMWeight = 1 ./ sqrt(sum(Smat.^2));</span>
2570199         CMWeight = 1./std(Smat)';
2580200
2590201         <span class="comment">% No weight</span>
2600202         <span class="comment">%CMWeight = ones(size(Smat,2),1);</span>
2610203
2620204         <span class="keyword">break</span>;
2630205         <span class="comment">%CMWeight = [CMWeight; ones(length(OCS.CM{i}.Data),1)];</span>
2640206
2650207     <span class="keyword">elseif</span> isscalar(OCS.CMWeight{i})
2660208         CMWeight = [CMWeight; ones(length(OCS.CM{i}.Data),1) * OCS.CMWeight{i}];
2670209
2680210     <span class="keyword">else</span>
2690211         CMWeight = [CMWeight; OCS.CMWeight{i}(:)];
2700212     <span class="keyword">end</span>
2710213 <span class="keyword">end</span>
2720214
2730215
2740216 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
2750217 <span class="comment">% Add a column weight to the response matrix %</span>
2760218 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
2770219 <span class="comment">%if isempty(U)</span>
2780220     <span class="keyword">for</span> j = 1:length(CMWeight)
2790221         Smat(:,j) = Smat(:,j) * CMWeight(j);
2800222     <span class="keyword">end</span>
2810223 <span class="comment">%end</span>
2820224
2830225
2840226 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
2850227 <span class="comment">% Add a weighted disperion as a column of the response matrix %</span>
2860228 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
2870229 <span class="keyword">if</span> OCS.FitRF == 0
2880230     <span class="comment">% Don't fit the RF frequency</span>
2890231     OCS.DeltaRF = 0;
2900232 <span class="keyword">elseif</span> any(OCS.FitRF == [1 2 3 4 5])
2910233     <span class="comment">% Column weighting (changing the units or sensitivity) seems to be a good thing.</span>
2920234     <span class="comment">% (At least at the ALS working in hardware units.)</span>
2930235     <span class="comment">% Weight by more than the average std of Smat also seems to give good results but</span>
2940236     <span class="comment">% I'm not sure it's necessary or should be done.</span>
2950237     RFWeight = 10 * mean(std(Smat)) / std(OCS.Eta);
2960238     Smat = [Smat RFWeight*OCS.Eta];
2970239 <span class="keyword">else</span>
2980240     error(<span class="string">'Unknown RF correction method.'</span>);
2990241 <span class="keyword">end</span>
3000242
3010243
3020244
3030245 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
3040246 <span class="comment">% Add a row weight to the response matrix %</span>
3050247 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
3060248 <span class="comment">% Weighted LS</span>
3070249 <span class="comment">% Weighted LS was developed for heteroscedastic errors but it's used here for other reasons</span>
3080250 <span class="comment">% W*(Mmeas - Mmodel) = W*A*b + W*e,  where the e's are gausian, zero mean, independent, but not constant variance</span>
3090251 <span class="comment">% W is chosen to make E(Wee'W') = sigma * I  (ie, W * e has a constant variance)</span>
3100252 <span class="comment">%</span>
3110253 <span class="comment">% The new LS equations are,</span>
3120254 <span class="comment">% b = inv(Amod'*W' * W*Amod) * Amod'*W' * W*(Xgoal - Xmeasured);      % Parameter fit</span>
3130255 <span class="comment">% b = V(:,Ivec) * b;</span>
3140256 <span class="comment">% bvar = inv(Amod'*W'*W*Amod);</span>
3150257 <span class="comment">%</span>
3160258 <span class="comment">% Since the weight matrix, W, is only going to be diagonal, it is less memory</span>
3170259 <span class="comment">% to scaled Amod and (Mmeas-Model) by the diagonal term and not create the matrix W</span>
3180260 <span class="keyword">if</span> ~all(BPMWeight == 1)
3190261     <span class="keyword">for</span> i = 1:size(Smat,1)
3200262         Smat(i,:) = BPMWeight(i) * Smat(i,:);
3210263     <span class="keyword">end</span>
3220264 <span class="keyword">end</span>
3230265
3240266
3250267 <span class="comment">%%%%%%%%%%%%%%%%</span>
3260268 <span class="comment">% The SVD Part %</span>
3270269 <span class="comment">%%%%%%%%%%%%%%%%</span>
3280270 <span class="keyword">if</span> isempty(U)
3290271     <span class="keyword">if</span> any(any(isnan(Smat)))
3300272         error(<span class="string">'The response matrix has at least one NaN.'</span>);
3310273     <span class="keyword">end</span>
3320274
3330275     <span class="comment">% Do the SVD</span>
3340276     [U, S, V] = svd(Smat, 0);  <span class="comment">%'econ');</span>
3350277     S = diag(S);
3360278 <span class="keyword">end</span>
3370279
3380280
3390281 <span class="comment">% Determine the singular vector and error check</span>
3400282 <span class="keyword">if</span> isnumeric(OCS.SVDIndex)
3410283     <span class="keyword">if</span> length(OCS.SVDIndex) == 1
3420284         <span class="keyword">if</span> rem(OCS.SVDIndex,1) == 0
3430285             OCS.SVDIndex = 1:OCS.SVDIndex;
3440286         <span class="keyword">else</span>
3450287             <span class="comment">% Use threshold</span>
3460288             Ivec = find(S &gt; max(S)*OCS.SVDIndex);
3470289             OCS.SVDIndex = Ivec;
3480290         <span class="keyword">end</span>
3490291     <span class="keyword">end</span>
3500292 <span class="keyword">elseif</span> ischar(OCS.SVDIndex)
3510293     <span class="keyword">if</span> strcmpi(OCS.SVDIndex, <span class="string">'All'</span>)
3520294         OCS.SVDIndex = 1:length(S);
3530295     <span class="keyword">else</span>
3540296         error(<span class="string">'OCS.SVDIndex unknown'</span>);
3550297     <span class="keyword">end</span>
3560298 <span class="keyword">else</span>
3570299     error(<span class="string">'OCS.SVDIndex unknown'</span>);
3580300 <span class="keyword">end</span>
3590301
3600302 <span class="comment">% Nothing greater then the total number is singular values</span>
3610303 i = find(OCS.SVDIndex &gt; length(S));
3620304 <span class="keyword">if</span> ~isempty(i)
3630305     warning(<span class="string">'Requested number of singular values is greater than the total.  Removing %d values.'</span>,length(i));
3640306     OCS.SVDIndex(i) = length(S);
3650307 <span class="keyword">end</span>
3660308
3670309 <span class="comment">% No non-integers</span>
3680310 OCS.SVDIndex = round(OCS.SVDIndex);
3690311
3700312 <span class="comment">% Nothing less than zero</span>
3710313 i = find(OCS.SVDIndex &lt;= 0);
3720314 OCS.SVDIndex(i) = length(S);
3730315
3740316 <span class="comment">% No repeats</span>
3750317 OCS.SVDIndex = sort(OCS.SVDIndex);
3760318 i = find(diff(OCS.SVDIndex)==0);
3770319 OCS.SVDIndex(i) = [];
3780320
3790321 <span class="keyword">if</span> isempty(OCS.SVDIndex)
3800322     error(<span class="string">'Singular values vector is empty'</span>);
3810323 <span class="keyword">end</span>
3820324
3830325
3840326 <span class="keyword">if</span> OCS.FitRF == 2 || OCS.FitRF == 3
3850327     <span class="comment">% Constrained L.S. with SVD</span>
3860328     <span class="comment">% Constraint: Sum of the energy change to zero</span>
3870329     <span class="comment">% Note: Model required</span>
3880330     
3890331     R = [];
3900332     HCMEnergyChangeTotal = 0;
3910333     <span class="keyword">for</span> i = 1:length(OCS.CM)
3920334         <span class="keyword">if</span> <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.CM{i}.FamilyName, <span class="string">'HCM'</span>)
3930335             L = <a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Circumference'</span>);
3940336
3950337             <span class="keyword">if</span> ~isfield(OCS.CM{i}, <span class="string">'Dispersion'</span>)
3960338                 <span class="comment">% Dispersion at the correctors in physics units</span>
3970339                 [OCS.CM{i}.Dispersion, DyHCM] = modeldisp([], OCS.CM{i}.FamilyName, OCS.CM{i}.DeviceList, <span class="string">'Numeric'</span>, <span class="string">'Physics'</span>);
3980340             <span class="keyword">end</span>
3990341             
4000342             <span class="comment">% Present corrector setpoint in physics units</span>
4010343             HCM = <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i});
4020344
4030345             
4040346             <span class="comment">% Move ALS specific code to setorbitdefault ???</span>
4050347             <span class="keyword">if</span> strcmpi(<a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Machine'</span>),<span class="string">'ALS'</span>)
4060348                 <span class="comment">% For the ALS, either the HCMCHICANE families need to be included or</span>
4070349                 <span class="comment">% the chicane &quot;part&quot; of the HCMs needs to be removed.  This will remove the chicane.</span>
4080350                 Energy = <a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Energy'</span>);
4090351
4100352                 <span class="comment">% Sector 6</span>
4110353                 <span class="comment">%                   Off    1.9 GeV   1.5 GeV</span>
4120354                 <span class="comment">% HCMCHICANEM(6,1)  80.0    18.0       ?</span>
4130355                 <span class="comment">% HCMCHICANEM(6,1)  80.0    20.0       ?</span>
4140356                 <span class="comment">% HCM(6,1)           0.0    18.8       ?</span>
4150357                 ihcm = findrowindex([6 1], OCS.CM{i}.DeviceList);
4160358                 <span class="keyword">if</span> length(ihcm) == 1
4170359                     <span class="keyword">try</span>
4180360                         <span class="keyword">if</span> <a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>(<span class="string">'HCMCHICANEM'</span>,[6 1]) &lt; 70
4190361                             <span class="comment">% Assume sector 6 chicane is on</span>
4200362                             <span class="keyword">if</span> Energy == 1.9
4210363                                 HCM.Data(ihcm) = HCM.Data(ihcm) + <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}.FamilyName, OCS.CM{i}.Field, -18.8, [6 1]);
4220364                             <span class="keyword">else</span>
4230365                                 HCM.Data(ihcm) = HCM.Data(ihcm) + <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}.FamilyName, OCS.CM{i}.Field,  -18.8*1.5/1.9, [6 1]);
4240366                             <span class="keyword">end</span>
4250367                         <span class="keyword">end</span>
4260368                     <span class="keyword">catch</span>
4270369                         fprintf(<span class="string">'%s\n'</span>, lasterr);
4280370                         fprintf(<span class="string">'Problem reading HCMCHICANEM(6,1).  The chicane &quot;offset&quot; on HCM(6,1) will not be removed.\n\n'</span>);
4290371                     <span class="keyword">end</span>
4300372                 <span class="keyword">end</span>
4310373
4320374                 <span class="comment">% Sector 11</span>
4330375                 <span class="comment">%                    Off    1.9 GeV   1.5 GeV</span>
4340376                 <span class="comment">% HCMCHICANEM(11,1)  80.0    40.5      52.0</span>
4350377                 <span class="comment">% HCMCHICANEM(11,1)  80.0    40.5      52.0</span>
4360378                 <span class="comment">% HCM(10,8)           0.0   -17.0     -14.0</span>
4370379                 <span class="comment">% HCM(11,1)           0.0   -17.0     -14.0</span>
4380380                 ihcm = findrowindex([10 8], OCS.CM{i}.DeviceList);
4390381                 <span class="keyword">if</span> length(ihcm) == 1
4400382                     <span class="keyword">try</span>
4410383                         <span class="keyword">if</span> <a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>(<span class="string">'HCMCHICANEM'</span>,[11 1]) &lt; 60
4420384                             <span class="comment">% Assume sector 11 chicane is on</span>
4430385                             <span class="keyword">if</span> Energy == 1.9
4440386                                 HCM.Data(ihcm) = HCM.Data(ihcm) + <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}.FamilyName, OCS.CM{i}.Field, [17; 17], [10 8;11 1]);
4450387                             <span class="keyword">else</span>
4460388                                 HCM.Data(ihcm) = HCM.Data(ihcm) + <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}.FamilyName, OCS.CM{i}.Field, [14; 14], [10 8;11 1]);
4470389                             <span class="keyword">end</span>
4480390                         <span class="keyword">end</span>
4490391                     <span class="keyword">catch</span>
4500392                         fprintf(<span class="string">'%s\n'</span>, lasterr);
4510393                         fprintf(<span class="string">'Due to an error, the chicane &quot;offset&quot; on HCM(10,8) will not be removed.n\n'</span>);
4520394                     <span class="keyword">end</span>
4530395                 <span class="keyword">end</span>
4540396                 ihcm = findrowindex([11 1], OCS.CM{i}.DeviceList);
4550397                 <span class="keyword">if</span> length(ihcm) == 1
4560398                     <span class="keyword">try</span>
4570399                         <span class="keyword">if</span> <a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>(<span class="string">'HCMCHICANEM'</span>,[11 1]) &lt; 60
4580400                             <span class="comment">% Assume sector 11 chicane is on</span>
4590401                             <span class="keyword">if</span> Energy == 1.9
4600402                                 HCM.Data(ihcm) = HCM.Data(ihcm) + <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}.FamilyName, OCS.CM{i}.Field, [17; 17], [10 8;11 1]);
4610403                             <span class="keyword">else</span>
4620404                                 HCM.Data(ihcm) = HCM.Data(ihcm) + <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}.FamilyName, OCS.CM{i}.Field, [14; 14], [10 8;11 1]);
4630405                             <span class="keyword">end</span>
4640406                         <span class="keyword">end</span>
4650407                     <span class="keyword">catch</span>
4660408                         fprintf(<span class="string">'%s\n'</span>, lasterr);
4670409                         fprintf(<span class="string">'Due to an error, the chicane &quot;offset&quot; on HCM(11,1) will not be removed.n\n'</span>);
4680410                     <span class="keyword">end</span>
4690411                 <span class="keyword">end</span>
4700412             <span class="keyword">end</span>
4710413             <span class="comment">% END ALS only</span>
4720414             
4730415             
4740416             <span class="comment">% Energy change to remove</span>
4750417             HCMEnergyChangeTotal = HCMEnergyChangeTotal + sum(-1 * HCM.Data .* OCS.CM{i}.Dispersion / <a href="getmcf.html" class="code" title="function Alpha = getmcf(ModelString)">getmcf</a> / L);
4760418             
4770419             <span class="comment">% Energy scaling must be in physics units, hence the (HCM.Data./OCS.CM{i}.Data)</span>
4780420             <span class="keyword">if</span> strcmpi(OCS.CM{i}.Units, <span class="string">'Hardware'</span>)
4790421                 <span class="comment">% PhysicsUnitsScaling = (HCM.Data./OCS.CM{i}.Data);  %This does not work for zero setpoints</span>
4800422                 PhysicsUnitsScaling = <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(OCS.CM{i}. FamilyName,OCS.CM{i}.Field, 1, OCS.CM{i}.DeviceList);
4810423             <span class="keyword">else</span>
4820424                 PhysicsUnitsScaling = 1;
4830425             <span class="keyword">end</span>
4840426             HCMEnergyChangeScalar = -1 * PhysicsUnitsScaling .* OCS.CM{i}.Dispersion / <a href="getmcf.html" class="code" title="function Alpha = getmcf(ModelString)">getmcf</a> / L;
4850427             R = [R HCMEnergyChangeScalar(:)'];
4860428             
4870429             <span class="comment">% Sum of correctors to zero</span>
4880430             <span class="comment">%R = [R ones(size(HCMEnergyChange(:)'))];</span>
4890431         <span class="keyword">else</span>
4900432             R = [R zeros(1,size(OCS.CM{i}.DeviceList,1))];
4910433         <span class="keyword">end</span>
4920434     <span class="keyword">end</span>
4930435     
4940436     <span class="keyword">if</span> OCS.FitRF == 2
4950437         <span class="comment">% Also remove the energy change due to the present corrector setpoints</span>
4960438         r = -HCMEnergyChangeTotal/2;    <span class="comment">% Not sure about the divide by 2</span>
4970439     <span class="keyword">elseif</span> OCS.FitRF == 3
4980440         <span class="comment">% Only remove energy change due to the incremental change in the correctors</span>
4990441         r = 0;
5000442     <span class="keyword">end</span>
5010443     
5020444     <span class="comment">% Add a zero for no constaint on the RF</span>
5030445     R = [R 0];
5040446     
5050447
5060448     <span class="comment">% Projection onto the columns of Smat*V (or U)</span>
5070449     <span class="comment">% Since the correctors are V*b, the restriction in corrector strength is R*V</span>
5080450     X = Smat * V(:,OCS.SVDIndex);
5090451     R = R * V(:,OCS.SVDIndex);
5100452     b = inv(X'*X)*X' * (BPMWeight .* (GoalOrbitVec - StartOrbitVec));
5110453     br = b - inv(X'*X)*R'*inv(R*inv(X'*X)*R')*(R*b-r);
5120454     Delcm = V(:,OCS.SVDIndex) * br;
5130455
5140456 <span class="keyword">elseif</span> OCS.FitRF == 4
5150457     <span class="comment">% Constrained L.S. with SVD</span>
5160458     <span class="comment">% Constraint: dot(DispersionX,   Smat * dHCM) = 0</span>
5170459     <span class="comment">%         or  dot(DispersionX' * Smat,  dHCM) = 0 or total dot(DisperionX, Smat * HCM)</span>
5180460     R = [];
5190461     DispDotSmatHCM = 0;
5200462     <span class="keyword">for</span> i = 1:length(OCS.BPM)
5210463         <span class="keyword">if</span> strcmpi(OCS.BPM{i}.FamilyName, <span class="string">'BPMx'</span>) || <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.BPM{i}.FamilyName, <span class="string">'BPMx'</span>)
5220464             <span class="keyword">if</span> ~isfield(OCS.BPM{i}, <span class="string">'Dispersion'</span>)
5230465                 <span class="comment">% Dispersion at the correctors in physics units</span>
5240466                 [OCS.BPM{i}.Dispersion, Dy] = <a href="getdisp.html" class="code" title="function [Data, FileName] = getdisp(varargin)">getdisp</a>(OCS.BPM{i}.FamilyName, OCS.BPM{i}.DeviceList, <span class="string">'Numeric'</span>,  OCS.BPM{i}.Units);
5250467             <span class="keyword">end</span>
5260468             R = [R OCS.BPM{i}.Dispersion(:)'];
5270469             
5280470             HCM = <a href="physics2hw.html" class="code" title="function S = physics2hw(Family, Field, value, DeviceList, Energy)">physics2hw</a>(OCS.CM{i});
5290471
5300472             
5310473             <span class="comment">% Move ALS specific code to setorbitdefault ???</span>
5320474             <span class="keyword">if</span> strcmpi(<a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Machine'</span>),<span class="string">'ALS'</span>)
5330475                 <span class="comment">% For the ALS, either the HCMCHICANE families need to be included or</span>
5340476                 <span class="comment">% the chicane &quot;part&quot; of the HCMs needs to be removed.  This will remove the chicane.</span>
5350477                 Energy = <a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'Energy'</span>);
5360478
5370479                 <span class="comment">% Sector 6</span>
5380480                 <span class="comment">%                   Off    1.9 GeV   1.5 GeV</span>
5390481                 <span class="comment">% HCMCHICANEM(6,1)  80.0    18.0       ?</span>
5400482                 <span class="comment">% HCMCHICANEM(6,1)  80.0    20.0       ?</span>
5410483                 <span class="comment">% HCM(6,1)           0.0    18.8       ?</span>
5420484                 ihcm = findrowindex([6 1], OCS.CM{i}.DeviceList);
5430485                 <span class="keyword">if</span> length(ihcm) == 1
5440486                     <span class="keyword">try</span>
5450487                         <span class="keyword">if</span> <a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>(<span class="string">'HCMCHICANEM'</span>,[6 1]) &lt; 70
5460488                             <span class="comment">% Assume sector 6 chicane is on</span>
5470489                             <span class="keyword">if</span> Energy == 1.9
5480490                                 HCM.Data(ihcm) = HCM.Data(ihcm) - 18.8;
5490491                             <span class="keyword">else</span>
5500492                                 HCM.Data(ihcm) = HCM.Data(ihcm) - 18.8*1.5/1.9;
5510493                             <span class="keyword">end</span>
5520494                         <span class="keyword">end</span>
5530495                     <span class="keyword">catch</span>
5540496                         fprintf(<span class="string">'%s\n'</span>, lasterr);
5550497                         fprintf(<span class="string">'Due to an error, the chicane &quot;offset&quot; on HCM(6,1) will not be removed.\n\n'</span>);
5560498                     <span class="keyword">end</span>
5570499                 <span class="keyword">end</span>
5580500
5590501                 <span class="comment">% Sector 11</span>
5600502                 <span class="comment">%                    Off    1.9 GeV   1.5 GeV</span>
5610503                 <span class="comment">% HCMCHICANEM(11,1)  80.0    40.5      52.0</span>
5620504                 <span class="comment">% HCMCHICANEM(11,1)  80.0    40.5      52.0</span>
5630505                 <span class="comment">% HCM(10,8)           0.0   -17.0     -14.0</span>
5640506                 <span class="comment">% HCM(11,1)           0.0   -17.0     -14.0</span>
5650507                 ihcm = findrowindex([10 8], OCS.CM{i}.DeviceList);
5660508                 <span class="keyword">if</span> length(ihcm) == 1
5670509                     <span class="keyword">try</span>
5680510                         <span class="keyword">if</span> <a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>(<span class="string">'HCMCHICANEM'</span>,[11 1]) &lt; 60
5690511                             <span class="comment">% Assume sector 11 chicane is on</span>
5700512                             <span class="keyword">if</span> Energy == 1.9
5710513                                 HCM.Data(ihcm) = HCM.Data(ihcm) + 17;
5720514                             <span class="keyword">else</span>
5730515                                 HCM.Data(ihcm) = HCM.Data(ihcm) + 14;
5740516                             <span class="keyword">end</span>
5750517                         <span class="keyword">end</span>
5760518                     <span class="keyword">catch</span>
5770519                         fprintf(<span class="string">'%s\n'</span>, lasterr);
5780520                         fprintf(<span class="string">'Problem reading HCMCHICANEM(11,1).  The chicane &quot;offset&quot; on HCM(10,8) will not be removed.n\n'</span>);
5790521                     <span class="keyword">end</span>
5800522                 <span class="keyword">end</span>
5810523                 ihcm = findrowindex([11 1], OCS.CM{i}.DeviceList);
5820524                 <span class="keyword">if</span> length(ihcm) == 1
5830525                     <span class="keyword">try</span>
5840526                         <span class="keyword">if</span> <a href="getsp.html" class="code" title="function [SP, tout, DataTime, ErrorFlag] = getsp(Family, varargin)">getsp</a>(<span class="string">'HCMCHICANEM'</span>,[11 1]) &lt; 60
5850527                             <span class="comment">% Assume sector 11 chicane is on</span>
5860528                             <span class="keyword">if</span> Energy == 1.9
5870529                                 HCM.Data(ihcm) = HCM.Data(ihcm) + 17;
5880530                             <span class="keyword">else</span>
5890531                                 HCM.Data(ihcm) = HCM.Data(ihcm) + 14;
5900532                             <span class="keyword">end</span>
5910533                         <span class="keyword">end</span>
5920534                     <span class="keyword">catch</span>
5930535                         fprintf(<span class="string">'%s\n'</span>, lasterr);
5940536                         fprintf(<span class="string">'Due to an error, the chicane &quot;offset&quot; on HCM(11,1) will not be removed.n\n'</span>);
5950537                     <span class="keyword">end</span>
5960538                 <span class="keyword">end</span>
5970539             <span class="keyword">end</span>
5980540             <span class="comment">% END ALS only</span>
5990541
6000542
6010543             <span class="comment">% Assumes 1 HCM family and it's first???</span>
6020544             DispDotSmatHCM = DispDotSmatHCM + OCS.BPM{i}.Dispersion'*Smat(1:size(OCS.BPM{i}.DeviceList,1),1:1:size(OCS.CM{i}.DeviceList,1))*HCM.Data;
6030545         <span class="keyword">else</span>
6040546             R = [R zeros(1,size(OCS.BPM{i}.DeviceList,1))];
6050547         <span class="keyword">end</span>
6060548     <span class="keyword">end</span>
6070549     <span class="comment">% [OCS.BPM{1}.Dispersion, Dy] = getdisp(OCS.BPM{1}.FamilyName, OCS.BPM{1}.DeviceList, 'Numeric', OCS.BPM{1}.Units);</span>
6080550     <span class="comment">% RR = [OCS.BPM{1}.Dispersion(:); zeros(size(OCS.BPM{2}.DeviceList,1),1)]';</span>
6090551
6100552     <span class="comment">%r = 0;</span>
6110553     r = -1*DispDotSmatHCM;  <span class="comment">% Not sure if this should be divided by 2</span>
6120554     R = [R*Smat(:,1:end-1) 0];
6130555
6140556
6150557     <span class="comment">% Projection onto the columns of Smat*V (or U)</span>
6160558     <span class="comment">% Since the correctors are V*b, the restriction in corrector strength is R*V</span>
6170559     X = Smat * V(:,OCS.SVDIndex);
6180560     R = R * V(:,OCS.SVDIndex);
6190561     b = inv(X'*X)*X' * (BPMWeight .* (GoalOrbitVec - StartOrbitVec));
6200562     br = b - inv(X'*X)*R'*inv(R*inv(X'*X)*R')*(R*b-r);
6210563     Delcm = V(:,OCS.SVDIndex) * br;
6220564     
6230565     <span class="comment">%DelcmNoR = V(:,OCS.SVDIndex) * b;</span>
6240566     <span class="comment">%fprintf('   dRF=%g  dRF(NoR)=%g\n',  RFWeight * Delcm(end),  RFWeight * DelcmNoR(end));</span>
6250567     
6260568 <span class="keyword">elseif</span> OCS.FitRF == 5
6270569     <span class="comment">% Constrained L.S. with SVD</span>
6280570     <span class="comment">% Constraint: dot(HCM(Dispersion), dHCM) = 0</span>
6290571     R = [];
6300572     <span class="keyword">for</span> i = 1:length(OCS.CM)
6310573         <span class="keyword">if</span> <a href="ismemberof.html" class="code" title="function  [IsTest, Index] = ismemberof(FamilyName, Field, MemberString)">ismemberof</a>(OCS.CM{i}.FamilyName, <span class="string">'HCM'</span>)
6320574             <span class="comment">% Assume CM{i} is the same plane as BPM{i}</span>
6330575             <span class="keyword">if</span> ~isfield(OCS.BPM{i}, <span class="string">'Dispersion'</span>)
6340576                 <span class="comment">% Dispersion at the horizontal BPMs</span>
6350577                 [OCS.BPM{i}.Dispersion, Dy] = <a href="getdisp.html" class="code" title="function [Data, FileName] = getdisp(varargin)">getdisp</a>(OCS.BPM{i}.FamilyName, OCS.BPM{i}.DeviceList, <span class="string">'Numeric'</span>,  OCS.BPM{i}.Units);
6360578
6370579                 <span class="comment">% Find the corrector pattern of the dispersion orbit</span>
6380580                 <span class="comment">%</span>
6390581                 <span class="comment">% OCS.SVDIndex ????</span>
6400582                 <span class="comment">% Smat???  make an array</span>
6410583                 <span class="comment">% Why is Rhcm so different from Rbpmx?</span>
6420584                 
6430585                 <span class="comment">% Assumes 1 HCM family and it's first???</span>
6440586                 SmatNoEta = Smat(1:size(OCS.BPM{i}.DeviceList,1),1:size(OCS.CM{i}.DeviceList,1));
6450587                 [UU, SS, VV] = svd(SmatNoEta, <span class="string">'econ'</span>);
6460588                 SS = diag(SS);
6470589
6480590                 X = SmatNoEta * VV(:,OCS.SVDIndex);
6490591                 b = inv(X'*X)*X' * OCS.BPM{i}.Dispersion;
6500592                 OCS.BPM{i}.DispersionCorrectors = VV(:,OCS.SVDIndex) * b;
6510593             <span class="keyword">end</span>
6520594             R = [R OCS.BPM{i}.DispersionCorrectors(:)'];
6530595         <span class="keyword">else</span>
6540596             R = [R zeros(1,size(OCS.CM{i}.DeviceList,1))];
6550597         <span class="keyword">end</span>
6560598     <span class="keyword">end</span>
6570599
6580600     r = 0;  <span class="comment">% Should it be an absolute conditional ???</span>
6590601     R = [R 0];
6600602
6610603
6620604     <span class="comment">% Projection onto the columns of Smat*V (or U)</span>
6630605     <span class="comment">% Since the correctors are V*b, the restriction in corrector strength is R*V</span>
6640606     X = Smat * V(:,OCS.SVDIndex);
6650607     R = R * V(:,OCS.SVDIndex);
6660608     b = inv(X'*X)*X' * (BPMWeight .* (GoalOrbitVec - StartOrbitVec));
6670609     br = b - inv(X'*X)*R'*inv(R*inv(X'*X)*R')*(R*b-r);
6680610     Delcm = V(:,OCS.SVDIndex) * br;
6690611     
6700612     <span class="comment">%DelcmNoR = V(:,OCS.SVDIndex) * b;</span>
6710613     <span class="comment">%fprintf('   dRF=%g  dRF(NoR)=%g\n',  RFWeight * Delcm(end),  RFWeight * DelcmNoR(end));</span>
6720614
6730615 <span class="keyword">else</span>
6740616
6750617     <span class="comment">% L.S. with SVD algorithm</span>
6760618     <span class="comment">% The V matrix columns are the singular vectors in the corrector magnet space</span>
6770619     <span class="comment">% The U matrix columns are the singular vectors in the BPM space</span>
6780620     <span class="comment">% Smat = U*S*V' where U'*U=I and V*V'=I</span>
6790621     <span class="comment">%</span>
6800622     <span class="comment">% b is the coef. of the projection onto the columns of Smat*V(:,Ivec) - the corrector space (same space as spanned by U)</span>
6810623     <span class="comment">% Sometimes it's interesting to look at the size of these coefficients with singular value number.</span>
6820624     b = U(:,OCS.SVDIndex)' * (BPMWeight .* (GoalOrbitVec - StartOrbitVec));
6830625     b = diag(S(OCS.SVDIndex).^(-1)) * b;
6840626
6850627     <span class="comment">% Convert the b vector back to coefficents of response matrix</span>
6860628     Delcm = V(:,OCS.SVDIndex) * b;
6870629 <span class="keyword">end</span>
6880630
6890631
6900632 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
6910633 <span class="comment">% Remove Weights for the Correction %</span>
6920634 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
6930635 <span class="comment">% Separate the RF from the correctors</span>
6940636 <span class="keyword">if</span> OCS.FitRF
6950637     <span class="comment">% RF frequency is the last actuator</span>
6960638     OCS.DeltaRF = RFWeight * Delcm(end);
6970639     <span class="keyword">if</span> strcmpi(<a href="getfamilydata.html" class="code" title="function [Data, ErrorFlag] = getfamilydata(Family, Field1, Field2, DeviceList)">getfamilydata</a>(<span class="string">'SubMachine'</span>), <span class="string">'Soleil'</span>)
6980640         <span class="comment">% Begin Add by Laurent % To be done with variables</span>
6990641         <span class="comment">% test for avoiding to move too often RF frequency</span>
7000642         <span class="comment">%fprintf('Prog %e \n',DelRF);</span>
7010643         <span class="keyword">if</span> abs(OCS.DeltaRF) &lt; DeltaRFmin <span class="comment">% minimum RF is .1 Hz</span>
7020644             OCS.DeltaRF = 0;
7030645         <span class="keyword">end</span>
7040646         <span class="keyword">if</span> abs(OCS.DeltaRF) &gt; DeltaRFmax <span class="comment">% maximum is 100 Hz</span>
7050647             error(<span class="string">'RF change too large DeltaRF = %f Hz (&gt; %f Hz) \nRF frequency shift has to be fixed first by hand.\nSee expert steprf command'</span>, <span class="keyword">...</span>
7060648                 <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(<span class="string">'RF'</span>,<span class="string">'Setpoint'</span>,OCS.DeltaRF), <a href="hw2physics.html" class="code" title="function S = hw2physics(Family, Field, value, DeviceList, Energy)">hw2physics</a>(<span class="string">'RF'</span>,<span class="string">'Setpoint'</span>,DeltaRFmax));
7070649         <span class="keyword">end</span>
7080650     <span class="keyword">end</span>
7090651     Delcm = Delcm(1:end-1);
7100652 <span class="keyword">end</span>
7110653
7120654 <span class="comment">% Remove column weights from Delcm</span>
7130655 Delcm = Delcm .* CMWeight;
7140656     
7150657     
7160658 <span class="comment">% Predict Orbit</span>
7170659 <span class="keyword">if</span> OCS.FitRF
7180660     OrbitDeltaTotal = [SmatNoWeights OCS.Eta] * [Delcm; OCS.DeltaRF];
7190661 <span class="keyword">else</span>
7200662     OrbitDeltaTotal = SmatNoWeights * Delcm;
7210663 <span class="keyword">end</span>
7220664 <span class="keyword">for</span> j = 1:length(OCS.BPM)
7230665     N = length(OCS.BPM{j}.Data);
7240666     OCS.BPM{j}.PredictedOrbitDelta = OrbitDeltaTotal(1:N);
7250667     OrbitDeltaTotal(1:N) = [];
7260668 <span class="keyword">end</span>
7270669
7280670
7290671 <span class="comment">% Don't put the Delcm into OCS.CM so this function can get</span>
7300672 <span class="comment">% called multiple times before actually correcting the orbit</span>
7310673 <span class="keyword">for</span> i = 1:length(OCS.CM)
7320674     N = length(OCS.CM{i}.Data);
7330675     <span class="comment">%OCS.CM{i}.Data = OCS.CM{i}.Data + Delcm(1:N);</span>
7340676     OCS.CM{i}.Delta = Delcm(1:N);
7350677     Delcm(1:N) = [];
7360678
7370679     OCS.CMWeight{i} = CMWeight(1:N);
7380680     CMWeight(1:N) = [];
7390681 <span class="keyword">end</span>
7400682
7410683
7420684 <span class="comment">% Remove the cell array if the length is one</span>
7430685 <span class="keyword">if</span> NoCellBPMFlag
7440686     OCS.BPM = OCS.BPM{1};
7450687     OCS.GoalOrbit = OCS.GoalOrbit{1};
7460688     OCS.BPMWeight = OCS.BPMWeight{1};
7470689 <span class="keyword">end</span>
7480690 <span class="keyword">if</span> NoCellCMFlag
7490691     OCS.CM = OCS.CM{1};
7500692     OCS.CMWeight = OCS.CMWeight{1};
7510693 <span class="keyword">end</span>
7520694
7530695
7540696 <span class="comment">% Output statistics???</span>
7550697 <span class="comment">% 1. Correctable orbit</span>
7560698 <span class="comment">% 2. Not correctable orbit</span>
7570699 <span class="comment">% 3. Error propagation</span>
7580700
7590701</pre></div>
760<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>
761</body>
762</html>
Note: See TracBrowser for help on using the repository browser.