Context Navigation

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

Last change on this file since 4 was 4, checked in by zhangj, 11 years ago
Initial import--MML version from SOLEIL@2013
File size: 38.3 KB

Line
1	function [AM, tout, DataTime, ErrorFlag] = getpvmodel(varargin)
2	%GETPVMODEL - Get the model value
3	% [Value, tout, DataTime, ErrorFlag] = getpvmodel(Family, Field, DeviceList, t)
4	% [Value, tout, DataTime, ErrorFlag] = getpvmodel(Family, DeviceList, t)
5	% [Value, tout, DataTime, ErrorFlag] = getpvmodel(DataStructure, t)
6	%
7	% INPUTS
8	% 1. Family - Family Name
9	% Data Structure
10	% Accelerator Object
11	% AT family name (or 'All' for every AT element)
12	% 'DCCT', 'TUNE', 'Energy'
13	% 2. Field - Accelerator Object field name ('Monitor', 'Setpoint', etc) {'Monitor'}
14	% AT field name
15	% or
16	% 'K','Quad','Quadrupole'
17	% 'K2','Sext','Sextupole'
18	% 'K3','Octupole'
19	% 'KS1','SkewQuad','Skew'
20	% 'Roll' or 'Tilt' - for a magnet rolls
21	% 'RollX' or 'TiltX' and 'RollY' or 'TiltY' - for corrector magnet rolls
22	% 'DX' or 'XShift' - for a magnet shift (see setshift)
23	% 'ClosedOrbit' - [x Px y Py dP dL] (dL is only calculated if the cavity is on)
24	%
25	% 'xTurns', 'PxTurns', 'yTurns', 'PyTurns', dpTurns', 'dLTurns' (single turn data) (BPM Number x N turns)
26	% 'FirstTurn' - 6-dimensional first turn or single pass
27	% 'Turns' - 3 dimensional matrix (BPM Number x N turns x 6)
28	% Note: Units can only be 'Physics' for these field types. 'xTurns' or 'yTurns' can
29	% be converted to 'Hardware' units by calling physics2hw as a separate call.
30	% Note: getpvmodel('TwissData','ClosedOrbit') - Gets the 6-dim start condition at the first AT element
31	% getpvmodel('TwissData', Field) - Gets the TwissData.(Field) twiss parameters at the first AT element
32	%
33	% 3. DeviceList ([Sector Device #] or [element #]) {Default: Entire family}
34	% 4. 'Physics' - Use physics units (optional override of units)
35	% 'Hardware' - Use hardware units (optional override of units) {Default: Hardware}
36	% 5. 'Struct' will return a data structure {Default for data structure inputs}
37	% 'Numeric' will return numeric outputs {Default for non-data structure inputs}
38	%
39	% OUTPUTS
40	% 1. Value - Model value
41	% 2. tout - Time it took to compute the output [seconds]
42	% 3. DataTime - Time (in seconds) since 00:00:00, Jan 1, 1970
43	% (seconds + milliseconds * i)
44	% 4. ErrorFlag
45	%
46	% NOTES
47	% 1. If Family is a cell array, then DeviceList and Field can also be a cell arrays
48
49	%
50	% Written by Gregory J. Portmann
51	% Revisions
52	% January 25, 2005 Laurent S. Nadolski
53	% ATparamgroup part modified to handle structure
54	% February 2005: BEND return always an energy (??? commented in 2007)
55	% 12/12/05 BPMx, BPMz, HCOR and VCOR getroll and getcrunch commented since too slow
56	% add in bpm and corrector other names (machine dependent?)
57
58	global THERING THERINGCELL
59
60	GCRFlag = 1; % Activate Gain Crunch Roll compulation
61
62	ErrorFlag = 0;
63
64	% Active NoiseFlag if BPM in Simulation mode (AO)
65	NoiseFlag = getfamilydata('BPMx','Simulated','NoiseActivated') \| ...
66	getfamilydata('BPMz','Simulated','NoiseActivated');
67	BPMsigma = 1e-6;
68
69	% Input parsing
70	StructOutputFlag = 0;
71	NumericOutputFlag = 0;
72	UnitsFlag = {};
73	for i = length(varargin):-1:1
74	if isstruct(varargin{i})
75	% Ignore structures
76	elseif iscell(varargin{i})
77	% Ignore cells
78	elseif strcmpi(varargin{i},'struct')
79	StructOutputFlag = 1;
80	varargin(i) = [];
81	elseif strcmpi(varargin{i},'numeric')
82	NumericOutputFlag = 1;
83	StructOutputFlag = 0;
84	varargin(i) = [];
85	elseif strcmpi(varargin{i},'Simulator') \|\| ...
86	strcmpi(varargin{i},'Model') \|\| ...
87	strcmpi(varargin{i},'Online') \|\| ...
88	strcmpi(varargin{i},'Manual')
89	% Remove and ignore
90	varargin(i) = [];
91	elseif strcmpi(varargin{i},'physics')
92	UnitsFlag = {'Physics'};
93	varargin(i) = [];
94	elseif strcmpi(varargin{i},'hardware')
95	UnitsFlag = {'Hardware'};
96	varargin(i) = [];
97	elseif strcmpi(varargin{i},'archive')
98	% Just remove
99	varargin(i) = [];
100	elseif strcmpi(varargin{i},'noarchive')
101	% Just remove
102	varargin(i) = [];
103	end
104	end
105
106	if isempty(varargin)
107	error('Must have at least a family name input');
108	else
109	Family = varargin{1};
110	if length(varargin) >= 2
111	Field = varargin{2};
112	end
113	if length(varargin) >= 3
114	DeviceList = varargin{3};
115	end
116	if length(varargin) >= 4
117	t = varargin{4};
118	else
119	t = 0;
120	end
121	end
122
123
124	%%%%%%%%%%%%%%%%%%%%%
125	% Cell Array Inputs %
126	%%%%%%%%%%%%%%%%%%%%%
127	if iscell(Family)
128	for i = 1:length(Family)
129	if length(varargin) < 2
130	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, UnitsFlag{:});
131	elseif length(varargin) < 3
132	if iscell(Field)
133	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, Field{i}, UnitsFlag{:});
134	else
135	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, Field, UnitsFlag{:});
136	end
137	else
138	if iscell(Field)
139	if iscell(DeviceList)
140	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, Field{i}, DeviceList{i}, UnitsFlag{:}, t);
141	else
142	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, Field{i}, DeviceList, UnitsFlag{:}, t);
143	end
144	else
145	if iscell(DeviceList)
146	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, Field, DeviceList{i}, UnitsFlag{:}, t);
147	else
148	[AM{i}, tout{i}, DataTime{i}, ErrorFlag{i}] = getpvmodel(Family{i}, Field, DeviceList, UnitsFlag{:}, t);
149	end
150	end
151	end
152	end
153	return
154	end
155
156
157	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
158	% Family or data structure inputs beyond this point %
159	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
160	if isstruct(Family)
161	% Only change StructOutputFlag if 'numeric' is not on the input line
162	if ~NumericOutputFlag
163	StructOutputFlag = 1;
164	end
165
166	% Data structure inputs
167	Field = '';
168	if length(varargin) >= 2
169	% Look for t as the second input
170	if isnumeric(varargin{2})
171	t = varargin{2};
172	else
173	Field = varargin{2};
174	end
175	end
176	if length(varargin) >= 3
177	DeviceList = varargin{3};
178	else
179	DeviceList = [];
180	end
181	if length(varargin) >= 4
182	t = varargin{4};
183	end
184	if isfield(Family,'FamilyName')
185	Family = Family.FamilyName;
186	else
187	error('For data structure inputs FamilyName field must exist')
188	end
189	if isempty(Field)
190	if isfield(Family,'Field')
191	Field = Family.Field;
192	end
193	end
194	if isempty(DeviceList)
195	if isfield(Family,'DeviceList')
196	DeviceList = Family.DeviceList;
197	end
198	end
199	if isempty(UnitsFlag)
200	if isfield(Family,'Units')
201	UnitsFlag{1} = Family.Units;
202	end
203	end
204	else
205	% Family string input
206	if length(varargin) < 2
207	Field = '';
208	end
209	if length(varargin) < 3
210	DeviceList = [];
211	end
212	if isnumeric(Field)
213	if length(varargin) >= 3
214	t = DeviceList;
215	end
216	DeviceList = Field;
217	Field = '';
218	end
219	end
220
221	if isempty(Field)
222	Field = 'Monitor';
223	end
224	if isempty(DeviceList)
225	if isfamily(Family)
226	DeviceList = family2dev(Family);
227	end
228	else
229	if isfamily(Family)
230	if (size(DeviceList,2) == 1)
231	DeviceList = elem2dev(Family, DeviceList);
232	end
233	end
234	end
235
236	if isempty(UnitsFlag)
237	UnitsFlag = 'Hardware';
238	else
239	UnitsFlag = UnitsFlag{1};
240	end
241
242
243	% Look to see it the AT model needs to be changed for this family
244	ATModelNumber = getfamilydata(Family, 'AT', 'ATModel');
245	if ~isempty(ATModelNumber)
246	% Change THERING
247	THERING = THERINGCELL{ATModelNumber};
248
249	% Set AD.Circumference
250	setfamilydata(findspos(THERING,length(THERING)+1), 'Circumference');
251
252	if isfield(THERING{1}, 'MachineType')
253	setfamilydata(THERING{1}.MachineType, 'MachineType');
254	end
255	if isfield(THERING{1}, 'Energy')
256	setfamilydata(THERING{1}.Energy, 'Energy');
257	end
258	if isfield(THERING{1}, 'InjectionEnergy')
259	setfamilydata(THERING{1}.InjectionEnergy, 'InjectionEnergy');
260	end
261	if isfield(THERING{1}, 'MCF')
262	setfamilydata(THERING{1}.MCF, 'MCF');
263	else
264	% Recompute the MCF if it's likely if a new accelerator
265	%try
266	% setfamilydata(getmcf('Model'), 'MCF');
267	%catch
268	%end
269	end
270	end
271
272
273	% Simulator (AT)
274	if isempty(THERING)
275	error('Simulator variable is not setup properly.');
276	end
277
278
279	t0 = gettime;
280
281	%%%%%%%%%%%%
282	% Get Data %
283	%%%%%%%%%%%%
284
285	% 20 hour lifetime, refill at midnight to 500 mamps
286	tau = 20;
287	BeamCurrent = 500;
288
289
290	% Families that do not require at AT field
291	if strcmp(Family, 'TUNE')
292	AM = modeltune;
293
294	% Add random errors
295	%AM = Tune1' + rand(2,1)*.0001;
296
297	elseif strcmp(Family, 'RF')
298
299	iCavity = findcells(THERING,'Frequency');
300	AM = THERING{iCavity(1)}.Frequency;
301
302	elseif strcmpi(Family, 'DCCT')
303
304	a = clock;
305	AM = BeamCurrent * exp(-(6060a(4)+60*a(5)+a(6))/tau/60/60);
306
307	elseif strcmpi(Family, 'LifeTime')
308
309	% 12 hour lifetime, refill at midnight to 500 mamps
310	AM = tau;
311	a = clock;
312
313	elseif any(strcmpi(Family,{'Energy','GeV'}))
314
315	AM = getenergymodel;
316	Family = 'GeV'; % Just so hw2physics works (ALS)
317
318	elseif strcmp(Family, 'TwissData')
319
320	if isfield(THERING{1}, 'TwissData')
321	TwissData = THERING{1}.TwissData;
322	else
323	TwissData = getfamilydata('TwissData');
324	end
325
326	if any(strcmpi(Field, {'xTurns','PxTurns','yTurns','PyTurns','dPTurns','dLTurns'}))
327	if ~isfield(TwissData, 'ClosedOrbit')
328	error('ClosedOrbit twiss parameter not found.');
329	end
330	if strcmpi(Field, 'xTurns')
331	AM = TwissData.ClosedOrbit(1);
332	elseif strcmpi(Field, 'PxTurns')
333	AM = TwissData.ClosedOrbit(2);
334	elseif strcmpi(Field, 'yTurns')
335	AM = TwissData.ClosedOrbit(3);
336	elseif strcmpi(Field, 'PyTurns')
337	AM = TwissData.ClosedOrbit(4);
338	elseif strcmpi(Field, 'dPTurns')
339	AM = TwissData.dP;
340	elseif strcmpi(Field, 'dLTurns')
341	AM = TwissData.dL;
342	end
343	else
344	if nargin == 1
345	AM = TwissData;
346	else
347	AM = TwissData.(Field);
348	end
349	end
350	return;
351
352	else
353	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354	% Families that require an AT field %
355	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356
357	% Find the index for where the desired data is in the total device list
358	if isfamily(Family)
359	DeviceListTotal = family2dev(Family, 0);
360	[DeviceIndex, iNotFound] = findrowindex(DeviceList, DeviceListTotal);
361	if length(iNotFound) > 0
362	% Device not found
363	for i = 1:length(iNotFound)
364	fprintf(' No devices to get for Family %s(%d,%d)\n', Family, DeviceList(i,1), DeviceList(i,2));
365	end
366	error(sprintf('%d Devices not found', length(iNotFound)));
367	end
368
369	% Find the AT structure if it exists
370	AT = getfamilydata(Family, Field, 'AT');
371	if isempty(AT)
372	if any(strcmpi(Field,{'Setpoint','Monitor','Sum','Set','Read','Readback'}))
373	AT = getfamilydata(Family, 'AT');
374	if isempty(AT)
375	% AT.ATType must be defined
376	%warning('Simulator not setup for %s.%s family (data filled with NaN)\n', Family, Field);
377	fprintf('WARNING: Simulator not setup for %s.%s family (data filled with NaN)\n', Family, Field);
378	AM = NaN * ones(length(DeviceIndex),1);
379	tout = gettime - t0;
380	ErrorFlag = 1;
381	DataTime = 0+0*sqrt(-1);
382	return
383	end
384	else
385	AT.ATType = Field;
386	end
387	end
388
389	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
390	% Special function for simulator %
391	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
392	if isfield(AT, 'SpecialFunctionGet')
393	AM = feval(AT.SpecialFunctionGet, Family, Field, DeviceList);
394
395	if strcmpi(UnitsFlag, 'Hardware')
396	if isfamily(Family, Field)
397	AM = physics2hw(Family, Field, AM, DeviceList, getenergymodel);
398	else
399	UnitsFlag = 'Physics';
400	end
401	else
402	UnitsFlag = 'Physics';
403	end
404
405	% Expand if there is a time vector input
406	if length(t) > 1
407	AM(:,1:length(t)) = AM * ones(1,length(t));
408	tout(1,1:length(t)) = t + gettime - t0;
409	DataTime(1:size(AM,1),1:length(t)) = fix(t0) + rem(t0,1)1e9sqrt(-1);
410	else
411	tout = t + gettime - t0;
412	DataTime = fix(t0) + rem(t0,1)1e9sqrt(-1);
413	end
414
415	if StructOutputFlag
416	% Structure output for channel name method
417	AM.Data = AM;
418	AM.FamilyName = Family;
419	AM.Field = Field;
420	AM.DeviceList = DeviceList;
421	AM.Mode = 'Simulator';
422	AM.t = t; % Matlab time at start of measurement
423	AM.tout = tout; % Matlab time at end of measurement
424	AM.DataTime = DataTime;
425	AM.TimeStamp = clock;
426	AM.Units = UnitsFlag;
427	if strcmpi(AM.Units,'Hardware')
428	AM.UnitsString = getfamilydata(Family, Field, 'HWUnits');
429	else
430	AM.UnitsString = getfamilydata(Family, Field, 'PhysicsUnits');
431	end
432	AM.DataDescriptor = [];
433	AM.CreatedBy = 'getpvmodel';
434	end
435
436	return
437	end
438	else
439	% Look for an AT family
440	if strcmpi(Family, 'All')
441	AT.ATIndex = (1:length(THERING))';
442	else
443	AT.ATIndex = findcells(THERING, 'FamName', Family);
444	end
445	DeviceIndex = 1:length(AT.ATIndex);
446	if isempty(DeviceList)
447	DeviceList = [ones(length(DeviceIndex),1) DeviceIndex(:)];
448	end
449	AT.ATType = Field;
450	end
451
452
453	% Get methods
454	if isfield(AT, 'ATParameterGroup')
455	ATIndexList = AT.ATIndex(DeviceIndex,1); % Only use the first column
456	for i = 1:size(DeviceList,1)
457	if isstruct(AT.ATParameterGroup)
458	% Introduced for SOLEIL
459	PGField = AT.ATParameterGroup(1).FieldName;
460	AM(i,1) = getfield(THERING{ATIndexList(1)}, PGField, ...
461	AT.ATParameterGroup(1).FieldIndex);
462
463	elseif iscell(AT.ATParameterGroup)
464	% If cell, get the parameter group
465	% Note: this only works if the first FieldName sets the group parameter value
466	PGField = AT.ATParameterGroup{DeviceIndex(i)}(1).FieldName;
467	AM(i,:) = getfield(THERING{ATIndexList(i)}, PGField, AT.ATParameterGroup{DeviceIndex(i)}(1).FieldIndex);
468
469	% If getfield disappears use:
470	%DataField = THERING{ATIndexList(i)}.(PGField);
471	%Index1 = AT.ATParameterGroup{DeviceIndex(i)}(1).FieldIndex{1};
472	%Index2 = AT.ATParameterGroup{DeviceIndex(i)}(1).FieldIndex{2};
473	%AM(i,1) = DataField(Index1, Index2);
474	elseif ischar(AT.ATParameterGroup)
475	% If string, get the field
476	PGField = AT.ATParameterGroup;
477	AM(i,:) = THERING{ATIndexList(i)}.PGField;
478	end
479	end
480
481	else
482
483	% Make sure AT.Index exists
484	if ~isfield(AT, 'ATIndex')
485	AT.ATIndex = family2atindex(Family, DeviceListTotal);
486	end
487
488	% For split magnet reduce the get to the first magnet in the split
489	% If K*Leff (not K) is being varied then use Nsplits as a multiplier
490	if size(AT.ATIndex,2) > 1
491	Nsplits = ones(size(AT.ATIndex));
492	Nsplits(find(isnan(AT.ATIndex)))=0;
493	Nsplits = sum(Nsplits')';
494	ATIndexList = AT.ATIndex(:,1);
495	else
496	ATIndexList = AT.ATIndex;
497	end
498	ATIndexList = ATIndexList(DeviceIndex);
499
500
501	% Check the DeviceIndex
502	if isempty(ATIndexList)
503	AM = [];
504	DataTime = 0+0*sqrt(-1);
505	return;
506	end
507
508
509	% Switch on simulation method
510	if any(strcmpi(AT.ATType, {'x','BPMx','Horizontal','y','BPMy', 'z','BPMz','Vertical','ClosedOrbit'}))
511	% Need to consider AT 'x' at say Family='Quad' -> rolls etc?
512	[ATIndexList, isort] = sort(ATIndexList);
513
514	[CavityState, PassMethod, iCavity] = getcavity;
515	if isempty(CavityState)
516	% No cavity in AT model
517	if isradiationon
518	fprintf(' Turning radiation off since there is no RF cavity.\n');
519	setradiation off;
520	end
521	Orbit = findsyncorbit(THERING, 0, ATIndexList);
522	else
523	if strcmpi(deblank(CavityState(1,:)), 'On') \|\| isradiationon
524	% findorbit6 recommended when cavity or radiation is on
525	if strcmpi(deblank(CavityState(1,:)), 'Off')
526	% Turn off cavity in AT model
527	fprintf(' Turning the RF cavity on, since radiation is on.\n');
528	setcavity('On');
529	end
530	Orbit = findorbit6(THERING, ATIndexList);
531	else
532	% Cavity is off in AT model
533	%setcavity('Off');
534
535	% This is a way to simulate the effect of the RF without a cavity
536	C = 2.99792458e8;
537	CavityFrequency = THERING{iCavity(1)}.Frequency;
538	CavityHarmNumber = THERING{iCavity(1)}.HarmNumber;
539	L = findspos(THERING,length(THERING)+1); % getfamilydata('Circumference')
540	f0 = C * CavityHarmNumber / L;
541	DeltaRF = CavityFrequency - f0; % Hz
542	Orbit = findsyncorbit(THERING, -CDeltaRFCavityHarmNumber/CavityFrequency^2, ATIndexList);
543	%Orbit = findsyncorbit(THERING, 0, ATIndexList);
544
545	% Reset cavity PassMethod
546	%setcavity(PassMethod);
547	end
548	end
549
550	if any(strcmpi(AT.ATType, {'x','BPMx','Horizontal'}))
551
552	% Roll and Crunch are corrected here (BPM gain errors are corrected in real2raw/raw2real)
553	if isfamily(Family, Field)
554	% Only corrector for gain/roll errors is using a ML family
555	% Roll and Crunch are corrected here (BPM gain errors are corrected in real2raw/raw2real)
556	NDevices = length(ATIndexList);
557	for i = 1:NDevices,
558	if isfield(THERING{ATIndexList(i)}, 'GCR')
559	GCR(i,:) = THERING{ATIndexList(i)}.GCR;
560	else
561	if GCRFlag
562	if i == 1
563	Crunch = getcrunch(Family, Field, DeviceList);
564	Roll = getroll(Family, Field, DeviceList);
565	end
566	m = gcr2loco(1, 1, Crunch(i), Roll(i));
567	GCR(i,:) = [1 1 Crunch(i) Roll(i)];
568	else % Laurent
569	GCR(i,:) = [1 1 0 0];
570	end
571	end
572	end
573
574	x = Orbit(1,:)'; %meters
575	y = Orbit(3,:)'; %meters
576
577
578	if NoiseFlag
579	x = x + randncut(size(x))*BPMsigma;
580	y = y + randncut(size(y))*BPMsigma;
581	end
582
583	AM = ((x - GCR(:,3) .* y) .* cos(GCR(:,4)) + (y + GCR(:,3) .* x) .* sin(GCR(:,4))) ./ sqrt(1-GCR(:,3).^2);
584
585	% Same as:
586	%Crunch = GCR(:,3);
587	%Roll = GCR(:,4);
588	%a = ( Crunch .* sin(Roll) + cos(Roll)) ./ sqrt(1 - Crunch.^2);
589	%b = (-Crunch .* cos(Roll) + sin(Roll)) ./ sqrt(1 - Crunch.^2);
590	%AM = a .* x + b .* y;
591	else
592	AM = Orbit(1,:)';
593	end
594
595	elseif any(strcmpi(AT.ATType, {'y','BPMy', 'z', 'BPMz', 'Vertical'}))
596
597	if isfamily(Family, Field)
598	% Only corrector for gain/roll errors is using a ML family
599	% Roll and Crunch are corrected here (BPM gain errors are corrected in real2raw/raw2real)
600	for i = 1:length(ATIndexList)
601	if isfield(THERING{ATIndexList(i)}, 'GCR')
602	GCR(i,:) = THERING{ATIndexList(i)}.GCR;
603	else
604	if GCRFlag
605	GCR(i,:) = [1 1 0 0]; % [Gx Gy Crunch Roll]
606	if i == 1
607	Crunch = getcrunch(Family, Field, DeviceList); %Laurent
608	Roll = getroll(Family, Field, DeviceList); % Laurent
609	end
610	m = gcr2loco(1, 1, Crunch(i), Roll(i));
611	GCR(i,:) = [1 1 Crunch(i) Roll(i)];
612	else % Laurent
613	GCR(i,:) = [1 1 0 0];
614	end
615	end
616	end
617
618	x = Orbit(1,:)';
619	y = Orbit(3,:)';
620
621	if NoiseFlag
622	x = x + randn(size(x))*1e-7;
623	y = y + randn(size(y))*1e-7;
624	end
625
626
627	AM = ((y - GCR(:,3) .* x) .* cos(GCR(:,4)) - (x + GCR(:,3) .* y) .* sin(GCR(:,4))) ./ sqrt(1-GCR(:,3).^2);
628
629	% Same as:
630	%Crunch = GCR(:,3);
631	%Roll = GCR(:,4);
632	%c = (-Crunch .* cos(Roll) - sin(Roll)) ./ sqrt(1 - Crunch.^2);
633	%d = (-Crunch .* sin(Roll) + cos(Roll)) ./ sqrt(1 - Crunch.^2);
634	%AM = c .* x + d .* y;
635	else
636	AM = Orbit(3,:)';
637	end
638	elseif strcmpi(AT.ATType, 'ClosedOrbit')
639	AM = Orbit';
640	end
641
642	AM(isort,:) = AM;
643
644	% Add noise to the BPMs
645	%AM = AM + .1e-6*randn(size(AM));
646
647	elseif any(strcmpi(AT.ATType, {'HCM', 'HCOR', 'FHCOR', 'CH'}))
648
649	if any(strcmpi(Field,{'Setpoint','Monitor','Sum','Set','Read','Readback','Kick'}))
650	% Bending Angle (Radians)
651	for i=1:length(ATIndexList)
652	% The CM gain errors are corrected in real2raw/raw2real
653	% Coupling: Magnet roll is part of the AT model
654	% The gain is part of hw2physics/physics2hw
655	if isfield(THERING{ATIndexList(i)}, 'Roll')
656	Roll = THERING{ATIndexList(i)}.Roll;
657	else
658	% Knowing the cross-plane family name can be a problem
659	% If the HCM family has the same AT index, then use it.
660	if GCRFlag
661	try
662	%HCMFamily = gethcmfamily;
663	HCMFamily = 'HCOR';
664	Roll = [0 getroll(Family, Field, DeviceList(i,:))];
665	HCMDevList = family2dev(HCMFamily);
666	iHCM = findrowindex(DeviceList(i,:), HCMDevList);
667	if ~isempty(iHCM)
668	ATIndexHCM = family2atindex(HCMFamily, DeviceList(i,:));
669	if ATIndexHCM == ATIndexList(i)
670	Roll = [getroll(HCMFamily, Field, DeviceList(i,:)) Roll(2)];
671	else
672	Roll = [0 0];
673	end
674	end
675	catch
676	Roll = [0 0];
677	end
678	else
679	Roll = [0 0];
680	end
681	end
682	AM(i,1) = [cos(Roll(2)) sin(Roll(2))] * THERING{ATIndexList(i)}.KickAngle(:) / (cos(Roll(1)-Roll(2)));
683
684	if size(AT.ATIndex,2) > 1
685	AM(i,1) = Nsplits(i) * AM(i,1);
686	end
687	end
688
689	% Add noise (microradian)
690	%AM = AM + 1e-6*randn(length(AM),1);
691	else
692	if isfield(THERING{ATIndexList(1)}, Field)
693	for i=1:length(ATIndexList)
694	AM(i,:) = THERING{ATIndexList(i)}.(Field);
695	end
696	else
697	AM = zeros(length(ATIndexList),1); % or error???
698	end
699	end
700
701
702	elseif any(strcmpi(AT.ATType, {'PZ'}))
703	AM = scriptSOFB_XBPM('Model', DeviceList);
704
705	elseif any(strcmpi(AT.ATType, {'VCM', 'VCOR', 'FVCOR', 'CV'}))
706
707	if any(strcmpi(Field,{'Setpoint','Monitor','Sum','Set','Read','Readback','Kick'}))
708	% Bending Angle (Radians)
709	for i=1:length(ATIndexList)
710	% The CM gain errors are corrected in real2raw/raw2real
711	% Coupling: Magnet roll is part of the AT model
712	% The gain is part of hw2physics/physics2hw
713	if isfield(THERING{ATIndexList(i)}, 'Roll')
714	Roll = THERING{ATIndexList(i)}.Roll;
715	else
716	% Knowing the cross-plane family name can be a problem
717	% If the VCM family has the same AT index, then use it.
718	try
719	if GCRFlag
720	%VCMFamily = getvcmfamily;
721	VCMFamily = 'VCOR';
722	Roll = [getroll(Family, Field, DeviceList(i,:)) 0];
723	VCMDevList = family2dev(VCMFamily);
724	iVCM = findrowindex(DeviceList(i,:), VCMDevList);
725	if ~isempty(iVCM)
726	ATIndexVCM = family2atindex(VCMFamily, DeviceList(i,:));
727	if ATIndexVCM == ATIndexList(i)
728	Roll = [Roll(1) getroll(VCMFamily, Field, DeviceList(i,:))];
729	else
730	Roll = [0 0];
731	end
732	end
733	else
734	Roll = [0 0];
735	end
736	catch err
737	Roll = [0 0];
738	end
739	end
740
741	AM(i,1) = [-sin(Roll(1)) cos(Roll(1))] * THERING{ATIndexList(i)}.KickAngle(:) / (cos(Roll(1)-Roll(2)));
742
743	if size(AT.ATIndex,2) > 1
744	AM(i,1) = Nsplits(i) * AM(i,1);
745	end
746	end
747	% Add noise (microradian)
748	%AM = AM + 1e-6*randn(length(AM),1);
749	else
750	if isfield(THERING{ATIndexList(1)}, Field)
751	for i=1:length(ATIndexList)
752	AM(i,:) = THERING{ATIndexList(i)}.(Field);
753	end
754	else
755	AM = zeros(length(ATIndexList),1);
756	end
757	end
758
759	elseif any(strcmpi(AT.ATType,{'K','Quad','Quadrupole'}))
760	% Quadrupole
761	if isfield(THERING{ATIndexList(1)}, 'K')
762	for i = 1:length(ATIndexList)
763	AM(i,1) = THERING{ATIndexList(i)}.K;
764	end
765	else
766	for i = 1:length(ATIndexList)
767	AM(i,1) = THERING{ATIndexList(i)}.PolynomB(2);
768	end
769	end
770	% Add noise
771	%AM = AM + 1e-3*randn(length(AM),1);
772
773	elseif any(strcmpi(AT.ATType,{'K2','Sext','Sextupole'}))
774	% Sextupole
775	for i = 1:length(ATIndexList)
776	AM(i,1) = THERING{ATIndexList(i)}.PolynomB(3);
777	end
778	% Add noise
779	%AM = AM + 1e-3*randn(length(AM),1);
780
781	elseif any(strcmpi(AT.ATType,{'K3','Octupole'}))
782	% Octupole
783	for i = 1:length(ATIndexList)
784	AM(i,1) = THERING{ATIndexList(i)}.PolynomB(4);
785	end
786	% Add noise
787	%AM = AM + 1e-3*randn(length(AM),1);
788
789	elseif any(strcmpi(AT.ATType,{'KS','KS1','SkewQ','SkewQuad', 'QT'}))
790	% SkewQuad
791	for i = 1:length(ATIndexList)
792	AM(i,1) = THERING{ATIndexList(i)}.PolynomA(2);
793	end
794	% Add noise
795	%AM = AM + 1e-3*randn(length(AM),1);
796
797	elseif strcmpi(AT.ATType, 'BEND')
798	% BEND
799	for i = 1:length(ATIndexList)
800	AM(i,1) = THERING{ATIndexList(i)}.BendingAngle;
801	% Modif: Laurent to be consistent, has to return an energy
802	% AM(i,1) = getenergymodel;
803	end
804
805	elseif strcmpi(AT.ATType, 'Photon BPM')
806	% Spear only
807	if strcmpi(Family, 'BLOpen')
808	AM = ones(length(DeviceIndex),1);
809	elseif strcmpi(Family, 'BLSum')
810	AM = ones(length(DeviceIndex),1);
811	elseif strcmpi(Family, 'BLErr')
812	% Calculate photon beam location
813	% for ii=1:length(DeviceIndex)
814	% if strcmp(AO.BLType,'dipole')
815	% AM(ii)=CalcDipolePBPM(orbit(3,BPMATIndex),orbit(4,UpstreamBPMIndex),AO.BLLength);
816	% elseif strcmp(AO.BLType,'id')
817	% AM(ii)=CalcIDAngle(DeviceIndex,'simulator');
818	% AM(ii)=CalcIDError(DeviceIndex,'simulator');
819	% end
820	% end
821	% Add noise
822	%AM = AM + 1e-3*randn(length(AM),1);
823	end
824
825	elseif strcmpi(AT.ATType, 'Kicker') \|\| strcmpi(AT.ATType, 'kickeramp')
826
827	ATIndexList = AT.ATIndex(DeviceIndex);
828	if any(strcmpi(Field,{'Setpoint','Monitor'}))
829	% KickAngle
830	for i=1:length(ATIndexList)
831	AM(i,1) = THERING{ATIndexList(i)}.KickAngle(1); % This only allows for a horizontal kick
832	end
833	else
834	if isfield(THERING{ATIndexList(1)}, Field)
835	for i=1:length(ATIndexList)
836	AM(i,:) = THERING{ATIndexList(i)}.(Field);
837	end
838	else
839	AM = NaN * ones(length(ATIndexList),1);
840	end
841	end
842
843	elseif strcmpi(AT.ATType, 'Roll') \|\| strcmpi(AT.ATType, 'Tilt')
844	% Roll or Tilt of a magnet
845	for i = 1:length(ATIndexList)
846	if isfield(THERING{ATIndexList(i)}, 'R1') && isfield(THERING{ATIndexList(i)}, 'R2')
847	R1 = THERING{ATIndexList(i)}.R1;
848	AM(i,1) = acos(R1(1,1));
849	else
850	error(sprintf('%s(%d,%d) must have a R1 & R2 field in the model to be rolled.', Family, DeviceList(i,:)));
851	end
852	end
853
854	elseif strcmpi(AT.ATType, 'XShift') \|\| strcmpi(AT.ATType, 'DX')
855	% Roll or Tilt of a magnet
856	for i = 1:length(ATIndexList)
857	if isfield(THERING{ATIndexList(i)}, 'T1')
858	T1 = THERING{ATIndexList(i)}.T1;
859	AM(i,1) = T1(1);
860	else
861	error(sprintf('%s(%d,%d) must have a T1 field in the model to be shifted.', Family, DeviceList(i,:)));
862	end
863	end
864
865	elseif strcmpi(AT.ATType, 'YShift') \|\| strcmpi(AT.ATType, 'DY')
866	% Roll or Tilt of a magnet
867	for i = 1:length(ATIndexList)
868	if isfield(THERING{ATIndexList(i)}, 'T1')
869	T1 = THERING{ATIndexList(i)}.T1;
870	AM(i,1) = T1(3);
871	else
872	error(sprintf('%s(%d,%d) must have a T1 field in the model to be shifted.', Family, DeviceList(i,:)));
873	end
874	end
875
876	elseif strcmpi(AT.ATType, 'RollX') \|\| strcmpi(AT.ATType, 'RollY')
877	% Roll or Tilt
878	for i=1:length(ATIndexList)
879	% Corrector magnet roll
880	% The .Roll field is just a middle layer way to store the roll.
881	% Otherwise, one can not tell the difference between x/y kicks and coupling
882	if isfield(THERING{ATIndexList(i)}, 'KickAngle')
883	if isfield(THERING{ATIndexList(i)}, 'Roll')
884	Roll = THERING{ATIndexList(i)}.Roll;
885	else
886	Roll = [0 0];
887	end
888	if strcmpi(AT.ATType, 'RollX')
889	AM(i,1) = Roll(1);
890	else
891	AM(i,1) = Roll(2);
892	end
893	else
894	error(sprintf('%s(%d,%d) must be a KickAngle field in the model to be rolled.', Family, DeviceList(i,:)));
895	end
896	end
897
898	elseif any(strcmpi(AT.ATType,{'FirstTurn','linepass','Turns', 'xTurns','PxTurns','yTurns','PyTurns','dPTurns','dLTurns'}))
899	% Turn-by-turn data
900
901	Nturns = round(max(abs(t))); % getfamilydata('NTURNS');
902	t = Nturns; % Just incase it changed
903	if isempty(Nturns) \|\| Nturns < 1
904	Nturns = 1;
905	end
906
907	% Initial launch condition
908	if isfield(THERING{1}, 'TwissData')
909	TwissData = THERING{1}.TwissData;
910	else
911	TwissData = getfamilydata('TwissData');
912	end
913	if ~isfield(TwissData, 'dP')
914	TwissData.dP = 0;
915	end
916	if ~isfield(TwissData, 'L')
917	TwissData.dL = 0;
918	end
919	if isempty(TwissData.ClosedOrbit)
920	fprintf(' 6-dim initial condition set to zeros. Use setpvmodel(''TwissData'',''ClosedOrbit'') to change it.\n');
921	TwissData.ClosedOrbit = [0 0 0 0]';
922	%x0 = [.001 0, 0.001, 0]';
923	TwissData.dP = 0;
924	TwissData.dL = 0;
925	end
926	x0 = [TwissData.ClosedOrbit(:); TwissData.dP; TwissData.dL];
927
928	[AM, ATIndexList, LostBeam] = getturns(x0, Nturns, Family, DeviceList);
929
930	if LostBeam
931	error('Tracking failed due to beam loss.');
932	end
933
934	if any(strcmpi(AT.ATType,{'FirstTurn', 'linepass', 'Turns'}))
935	%AM = squeeze(AM(:,:,1:6));
936	elseif strcmpi(AT.ATType,'xTurns')
937	AM = squeeze(AM(:,:,1));
938	elseif strcmpi(AT.ATType,'PxTurns')
939	AM = squeeze(AM(:,:,2));
940	elseif strcmpi(AT.ATType,'yTurns')
941	AM = squeeze(AM(:,:,3));
942	elseif strcmpi(AT.ATType,'PyTurns')
943	AM = squeeze(AM(:,:,4));
944	elseif strcmpi(AT.ATType,'dPTurns')
945	AM = squeeze(AM(:,:,5));
946	elseif strcmpi(AT.ATType,'dLTurns')
947	AM = squeeze(AM(:,:,6));
948	else
949	error('Input of unknown type');
950	end
951
952	%if size(AM, 1) ~= length(ATIndexList)
953	if size(AM, 2) ~= Nturns
954	AM = AM';
955	end
956
957	% Gain/Roll coordinate change???
958
959	% Add noise
960	%AM = AM + 1e-3*randn(length(AM),1);
961
962	elseif strcmpi(AT.ATType, 'Septum')
963
964	AM = 0;
965
966	elseif strcmpi(AT.ATType, 'null')
967	% JR default do-nothing behaviour
968	AM = NaN;
969
970	else
971	if isfield(THERING{ATIndexList(1)}, Field)
972	for i = 1:length(ATIndexList)
973	AM(i,:) = THERING{ATIndexList(i)}.(Field);
974	end
975	else
976	%error(sprintf('ATType unknown for Family %s', Family));
977	AM = NaN * ones(length(ATIndexList),1);
978	end
979	% Add noise
980	%AM = AM + 1e-3*randn(length(AM),1);
981	end
982	end
983	end
984
985
986	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
987	% Change to hardware units if requested %
988	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
989	if strcmpi(UnitsFlag, 'Hardware')
990	if isfamily(Family, Field)
991	AM = physics2hw(Family, Field, AM, DeviceList, getenergymodel);
992	else
993	persistent WarningFlag
994	if isempty(WarningFlag)
995	fprintf('\n You are asking for hardware units from a nonstandard MML family (%s.%s).\n', Family, Field);
996	fprintf(' Whenever conversion factors are not known for model data, physics units are returned.\n\n');
997	WarningFlag = 1;
998	end
999
1000	UnitsFlag = 'Physics';
1001	end
1002	else
1003	UnitsFlag = 'Physics';
1004	end
1005
1006
1007	% Expand if there is a time vector input
1008	if length(t) > 1
1009	AM(:,1:length(t)) = AM * ones(1,length(t));
1010	tout(1,1:length(t)) = t + gettime - t0;
1011	else
1012	if strcmpi(Field,'Turns')
1013	tout = gettime - t0;
1014	else
1015	tout = t + gettime - t0;
1016	end
1017	%DataTime = fix(t0) + rem(t0,1)1e9sqrt(-1);
1018	end
1019	DataTime(1:size(AM,1),1:length(t)) = fix(t0) + rem(t0,1)1e9sqrt(-1);
1020
1021
1022	%%%%%%%%%%%%%%%%%%%%%
1023	% Structure Outputs %
1024	%%%%%%%%%%%%%%%%%%%%%
1025	if StructOutputFlag
1026	% Structure output for channel name method
1027	AM.Data = AM;
1028	AM.FamilyName = Family;
1029	AM.Field = Field;
1030	AM.DeviceList = DeviceList;
1031	AM.Mode = 'Simulator';
1032	AM.t = t; % Matlab time at start of measurement
1033	AM.tout = tout; % Matlab time at end of measurement
1034	AM.DataTime = DataTime;
1035	AM.TimeStamp = clock;
1036	AM.Units = UnitsFlag;
1037	if strcmpi(AM.Units,'Hardware')
1038	AM.UnitsString = getfamilydata(Family, Field, 'HWUnits');
1039	else
1040	AM.UnitsString = getfamilydata(Family, Field, 'PhysicsUnits');
1041	end
1042	AM.DataDescriptor = [];
1043	AM.CreatedBy = 'getpvmodel';
1044	end
1045

Note: See TracBrowser for help on using the repository browser.

Download in other formats: