Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_ltp_ifo2acc.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % MDC1 Conversion to acceleration
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Test script for the MDC1 conversion to acceleration of IFO data
4 %
5 % $Id: test_ao_ltp_ifo2acc.m,v 1.6 2009/12/10 15:35:38 luigi Exp $
6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 %% Clear
8
9 clear all
10
11 %% Loading data
12
13 currPath = cd;
14 cd 'C:\Users\Luigi\Dati\mock_data\interferometer\';
15
16 j=1;
17 fs = 10;
18
19 disp(sprintf('*** Reading data %d', j));
20 in = load(sprintf('mockdata_16_48_17_11_2007_%d.dat', j));
21 d1 = in(:,2);
22 d12 = in(:,3);
23
24 % Build O1
25 eval('ts = tsdata(d1, fs);');
26 eval(sprintf('o1_%d = ao(ts);', j));
27 eval(sprintf('o1_%d = setName(o1_%d, ''o1'');', j, j));
28 eval(sprintf('o1_%d = setXunits(o1_%d, ''s'');', j, j));
29 eval(sprintf('o1_%d = setYunits(o1_%d, ''m'');', j, j));
30
31 % Build O12
32 eval('ts = tsdata(d12, fs);');
33 eval(sprintf('o12_%d = ao(ts);', j));
34 eval(sprintf('o12_%d = setName(o12_%d, ''o12'');', j, j));
35 eval(sprintf('o12_%d = setXunits(o12_%d, ''s'');', j, j));
36 eval(sprintf('o12_%d = setYunits(o12_%d, ''m'');', j, j));
37
38 cd(currPath)
39
40 eval(sprintf('ch1 = o1_%d;',j));
41 eval(sprintf('ch2 = o12_%d;',j));
42
43 %% Save
44
45 % save(ch1, 'C:\Dati\mock_data\interferometer\mdc1_ch1.mat')
46 % save(ch2, 'C:\Dati\mock_data\interferometer\mdc1_ch12.mat')
47
48 %% Check data spectra
49
50 ch1xx = ch1.psd;
51 ch2xx = ch2.psd;
52 iplot(ch1xx,ch2xx)
53
54 %% Convert to acceleration
55
56 % default settings
57 [a1,ad] = ltp_ifo2acc(ch1,ch2);
58
59 %% check spectra
60
61 a1xx = a1.psd;
62 adxx = ad.psd;
63 iplot(a1xx,adxx)
64
65
66 %% ************************************************************************
67
68 %% Set useful params
69 % corresponding to default values
70
71 fs = ch1.fs; % Hz - sampling frequency
72
73 TAU_th = 0.1; % s - time constant of thrusters actuation
74 DELTAT_Th = 0.315; % s - delay introduced to command the force to thrusters system
75
76 TAU_el = 0.01; % s - time constant of electrostatic actuation
77 DELTAT_el = 0.305; % s - delay introduced to command the force to electrostatic actuators
78
79 w1 = 0.00114i; % TM1 stiffness
80 w2 = 0.00141i; % TM2 stiffness
81
82 S11 = 1; % Calibration of IFO CH1
83 S1D = 0; % Cross-talk diff channel to channel 1
84 SD1 = -1e-4; % Cross-talk channel 1 to diff channel
85 SDD = 1; % Calibration of IFO diff channel
86
87 m1 = 1.96; % TM1 mass
88 m2 = 1.96; % TM2 mass
89 msc = 422.7; % SC mass
90 Gamma = 4.92e-9; % Gravitational gradient
91
92
93 %% Do step-by-step Checks
94
95 % get mdc1 ifo noise psd models
96 currPath = cd;
97 cd 'C:\Users\Luigi\ltp_data_analysis\MDCs\MDC1_UTN\';
98
99
100 So = matrix([cd '\CSD_mod.mat']);
101
102 cd(currPath)
103
104 % frequency range for TFs calculation
105 f = So.objs(1,1).x;
106
107 % set units for So
108 for ii=1:numel(So.objs)
109 So.objs(ii).setYunits('m^2 Hz^-1');
110 end
111
112
113 %% Get intermediate steps
114
115 % get just Free Dynamics
116 [a1fd,adfd] = ltp_ifo2acc(ch1,ch2,plist('Hdf',0,'Hsus',0,'Adf',0,'Asus',0));
117
118 % get actuated force per unit of mass
119 [a1c,adc] = ltp_ifo2acc(ch1,ch2,plist('FreeDyn_1',0,'FreeDyn_Delta',0));
120
121
122 %% IFO TFs ***
123
124 % Getting models
125 IFO11 = smodel('S11');
126 IFO11.setParams({'S11'},{S11});
127 IFO11.setXvar('f');
128 IFO11.setXvals(f);
129 IFO11.setXunits('Hz');
130 IFO11.setName;
131
132 IFO1D = smodel('S1D');
133 IFO1D.setParams({'S1D'},{S1D});
134 IFO1D.setXvar('f');
135 IFO1D.setXvals(f);
136 IFO1D.setXunits('Hz');
137 IFO1D.setName;
138
139 IFOD1 = smodel('SD1');
140 IFOD1.setParams({'SD1'},{SD1});
141 IFOD1.setXvar('f');
142 IFOD1.setXvals(f);
143 IFOD1.setXunits('Hz');
144 IFOD1.setName;
145
146 IFODD = smodel('SDD');
147 IFODD.setParams({'SDD'},{SDD});
148 IFODD.setXvar('f');
149 IFODD.setXvals(f);
150 IFODD.setXunits('Hz');
151 IFODD.setName;
152
153 % Build the matrix
154 IFO = matrix([IFO11 IFO1D;IFOD1 IFODD]);
155
156
157 %% Free Dynamics TFs ***
158
159 % Getting models
160 FD11 = smodel('(2.*pi.*1i.*f).^2 + (1+m1/msc)*(w1^2) + (m2/msc)*(w2^2)');
161 FD11.setParams({'m1','m2','msc','w1','w2'},{m1,m2,msc,w1,w2});
162 FD11.setXvar('f');
163 FD11.setXvals(f);
164 FD11.setYunits('s^-2');
165 FD11.setXunits('Hz');
166 FD11.setName;
167
168 FD12 = smodel('(m2/msc)*(w2^2) + Gamma');
169 FD12.setParams({'m2','msc','w2','Gamma'},{m2,msc,w2,Gamma});
170 FD12.setXvar('f');
171 FD12.setXvals(f);
172 FD12.setYunits('s^-2');
173 FD12.setXunits('Hz');
174 FD12.setName;
175
176 FD21 = smodel('(w2^2) - (w1^2)');
177 FD21.setParams({'w1','w2'},{w1,w2});
178 FD21.setXvar('f');
179 FD21.setXvals(f);
180 FD21.setYunits('s^-2');
181 FD21.setXunits('Hz');
182 FD21.setName;
183
184 FD22 = smodel('(2.*pi.*1i.*f).^2 + (w2^2) - 2*Gamma');
185 FD22.setParams({'w2','Gamma'},{w2,Gamma});
186 FD22.setXvar('f');
187 FD22.setXvals(f);
188 FD22.setYunits('s^-2');
189 FD22.setXunits('Hz');
190 FD22.setName;
191
192 % Build the matrix
193 FD = matrix([FD11 FD12;FD21 FD22]);
194
195 % Transfer function from IFO output --> total acceleration
196 FDO = FD*inv(IFO);
197
198 FDO11 = eval(FDO.objs(1,1));
199 FDO12 = eval(FDO.objs(1,2));
200 FDO21 = eval(FDO.objs(2,1));
201 FDO22 = eval(FDO.objs(2,2));
202
203 FDOe = matrix([FDO11 FDO12;FDO21 FDO22]);
204
205 % expected PSD for total acceleration
206 Sfd = FDOe*So*conj(transpose(FDOe));
207
208 % eval results
209 Sfd11 = (2/fs).*abs(Sfd.objs(1,1));
210 Sfd12 = (2/fs).*Sfd.objs(1,2);
211 Sfd22 = (2/fs).*abs(Sfd.objs(2,2));
212
213 % get PSD from data
214 plpsd = plist('nfft',1e5);
215 Pfd11 = psd(a1fd,plpsd);
216 Pfd22 = psd(adfd,plpsd);
217 % plcsd = plist('nfft',5e3);
218 % Pfd12 = csd(a1fd,adfd,plcsd);
219
220 %% plot
221 plpl1 = plist('Legends', {'S_{a1a1} data','S_{a1a1}','S_{aDaD} data','S_{aDaD}'},...
222 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
223 iplot(sqrt(Pfd11),sqrt(Sfd11),sqrt(Pfd22),sqrt(Sfd22),plpl1)
224 % iplot(Pfd11,Sfd11,Pfd22,Sfd22,plpl1)
225
226 % plpl2 = plist('Legends', {'S_{a1aD} data','S_{a1aD}'},...
227 % 'LineStyles', {'', '--'},'LineWidths',{'All',3});
228 % iplot(Pfd12,Sfd12,plpl2)
229
230
231 %% Controllers TFs ***
232
233 % Get expected transfer functions
234 Cdf = miir(plist('built-in','DFACS_x_Hdf'));
235 Cdfr = resp(Cdf,plist('f',f,'bank','parallel'));
236
237 Csus = miir(plist('built-in','DFACS_M3_x_Hsus'));
238 Csusr = resp(Csus,plist('f',f,'bank','parallel'));
239
240 % Control matrix
241 H = matrix([Cdfr ao(0);ao(0) Csusr]);
242
243 %% Actuation TFs ***
244 % Get expected transfer functions
245 % ThActr = ao(plist('fsfcn', '10.*(1./(1 + (2.*pi.*1i.*f)*0.1)).*(1./(2.*pi.*1i.*f)).*(exp(-(2.*pi.*1i.*f).*0.315)-exp(-(2.*pi.*1i.*f).*(0.315+0.1)))', 'f', f));
246 % ElActr = ao(plist('fsfcn', '10.*(1./(1 + (2.*pi.*1i.*f)*0.01)).*(1./(2.*pi.*1i.*f)).*(exp(-(2.*pi.*1i.*f).*0.305)-exp(-(2.*pi.*1i.*f).*(0.305+0.1)))', 'f', f));
247
248 % Get expected transfer functions
249 pl = plist('built-in','mdc1_Adf');
250 Adf = miir(pl);
251 Adfr = resp(Adf,plist('f',f));
252
253 pl = plist('built-in','mdc1_Asus');
254 Asus = miir(pl);
255 Asusr = resp(Asus,plist('f',f));
256
257 % Actuation matrix
258 A = matrix([Adfr ao(0);ao(0) Asusr]);
259
260 %% Get expected spectra and measured
261
262 C = matrix([ao(1) ao(-m2/msc);ao(0) ao(-1)]);
263 Hnorm = matrix([ao(1/msc) ao(0);ao(0) ao(1/m2)]);
264 tfc = C*A*Hnorm*H;
265 for ii =1:numel(tfc.objs)
266 tfc.objs(ii).setYunits(tfc.objs(ii).yunits/unit('kg'));
267 tfc.objs(ii).simplifyYunits;
268 end
269 Sc = tfc*So*conj(transpose(tfc));
270
271 % get PSD from data
272 plpsd = plist('nfft',1e5);
273 Pc11 = psd(a1c,plpsd);
274 Pc22 = psd(adc,plpsd);
275 % plcsd = plist('nfft',5e3);
276 % Pc12 = csd(a1c,adc,plcsd);
277
278 % plot
279 plpl1 = plist('Legends', {'S_{a1a1} data','S_{a1a1}','S_{aDaD} data','S_{aDaD}'},...
280 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
281 % iplot(Pc11,(2/fs).*abs(Sc.objs(1,1)),Pc22,(2/fs).*abs(Sc.objs(2,2)),plpl1)
282 iplot(sqrt(Pc11),sqrt((2/fs).*abs(Sc.objs(1,1))),sqrt(Pc22),sqrt((2/fs).*abs(Sc.objs(2,2))),plpl1)
283
284 % plpl2 = plist('Legends', {'S_{a1aD} data','S_{a1aD}'},...
285 % 'LineStyles', {'', '--'},'LineWidths',{'All',3});
286 % iplot(Pc12,Sc.objs(1,2),plpl2)
287
288
289 %% Full model of acceleration noise
290
291 TF = FDOe + tfc;
292 Stot = TF*So*conj(transpose(TF));
293
294 % plot
295 plpl = plist('Legends', {'Ch1 data','Ch1','Ch12 data','Ch12'},...
296 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
297 iplot(sqrt(a1xx),sqrt((2/fs).*abs(Stot.objs(1,1))),sqrt(adxx),sqrt((2/fs).*abs(Stot.objs(2,2))),plpl)
298 % iplot(a1xx,(2/fs).*abs(Stot.objs(1,1)),adxx,(2/fs).*abs(Stot.objs(2,2)),plpl)
299
300
301
302 %% Test with zero input
303
304 nsecs = ch1.nsecs;
305
306 ch1_0 = ao(plist('tsfcn', 'zeros(size(t))', 'fs', 10, 'nsecs', nsecs, 'yunits','m'));
307 ch2_0 = ao(plist('tsfcn', 'zeros(size(t))', 'fs', 10, 'nsecs', nsecs, 'yunits','m'));
308
309 %% %%% %%% %%%
310
311 [a1_0,ad_0] = ltp_ifo2acc(ch1_0,ch2_0);
312
313 % output must be zero
314 plpsd = plist('nfft',1e5);
315 a1_0xx = psd(a1_0,plpsd);
316 ad_0xx = psd(ad_0,plpsd);
317
318 % plot
319 plpl = plist('Legends', {'Ch1 data','Ch12 data'},...
320 'LineStyles', {'', ''},'LineWidths',{'All',3});
321 iplot(sqrt(a1_0xx),sqrt(ad_0xx),plpl)
322
323 %% %%% %%% %%%
324
325 [a1_1,ad_1] = ltp_ifo2acc(ch1,ch2_0);
326
327 TF = FDOe + tfc;
328 So1 = copy(So,1);
329 for ii = 2:numel(So1.objs)
330 So1.objs(ii).setY(zeros(size(So1.objs(ii).x)));
331 end
332 Stot1 = TF*So1*conj(transpose(TF));
333
334 plpsd = plist('nfft',1e5);
335 a1_1xx = psd(a1_1,plpsd);
336 ad_1xx = psd(ad_1,plpsd);
337
338 % plot
339 plpl = plist('Legends', {'Ch1 data','Ch1','Ch12 data','Ch12'},...
340 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
341 iplot(sqrt(a1_1xx),sqrt((2/fs).*abs(Stot1.objs(1,1))),sqrt(ad_1xx),sqrt((2/fs).*abs(Stot1.objs(2,2))),plpl)
342
343 %% %%% %%% %%%
344
345 [a1_2,ad_2] = ltp_ifo2acc(ch1_0,ch2);
346
347 TF = FDOe + tfc;
348 So2 = copy(So,1);
349 for ii = 1:numel(So2.objs)-1
350 So2.objs(ii).setY(zeros(size(So2.objs(ii).x)));
351 end
352 Stot2 = TF*So2*conj(transpose(TF));
353
354 plpsd = plist('nfft',1e5);
355 a1_2xx = psd(a1_2,plpsd);
356 ad_2xx = psd(ad_2,plpsd);
357
358 % plot
359 plpl = plist('Legends', {'Ch1 data','Ch1','Ch12 data','Ch12'},...
360 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
361 iplot(sqrt(a1_2xx),sqrt((2/fs).*abs(Stot2.objs(1,1))),sqrt(ad_2xx),sqrt((2/fs).*abs(Stot2.objs(2,2))),plpl)
362
363
364
365
366
367
368
369
370
371