Mercurial > hg > ltpda

comparison m-toolbox/test/MDC1/test_mdc1_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_mdc1_ifo2acc.m,v 1.2 2009/09/24 15:17:36 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 %% Check data spectra
44
45 ch1xx = ch1.psd;
46 ch2xx = ch2.psd;
47 iplot(ch1xx,ch2xx)
48
49 %% Set useful params
50
51 fs = ch1.fs; % Hz - sampling frequency
52
53 TAU_th = 0.1; % s - time constant of thrusters actuation
54 DELTAT_Th = 0.315; % s - delay introduced to command the force to thrusters system
55
56 TAU_el = 0.01; % s - time constant of electrostatic actuation
57 DELTAT_el = 0.305; % s - delay introduced to command the force to electrostatic actuators
58
59 w1 = 0.00114i; % TM1 stiffness
60 w2 = 0.00141i; % TM2 stiffness
61
62 S11 = 1; % Calibration of IFO CH1
63 S1D = 0; % Cross-talk diff channel to channel 1
64 SD1 = -1e-4; % Cross-talk channel 1 to diff channel
65 SDD = 1; % Calibration of IFO diff channel
66
67 %% Load built in model for Drag Free
68
69 pl = plist('built-in','mdc1_Cdf','gain',1);
70 Cdf = rational(pl);
71
72 %% Load Built in model for Thrusters
73
74 pl = plist('built-in','mdc1_actuation','fs',fs,'tau',TAU_th,'DeltaT',DELTAT_Th);
75 ThAct = miir(pl);
76
77 %% Load built in model for Suspension controller
78
79 pl = plist('built-in','mdc1_Csus','gain',1);
80 Csus = rational(pl);
81
82 %% Load Built in model for Electrostatic actuators
83
84 pl = plist('built-in','mdc1_actuation','fs',fs,'tau',TAU_el,'DeltaT',DELTAT_el);
85 ElAct = miir(pl);
86
87 %% Convert to acceleration
88
89 pl = plist('w1',w1,'w2',w2,'S11',S11,'S1D',S1D,'SD1',SD1,'SDD',SDD,...
90 'Cdf',Cdf,'Csus',Csus,'ThAct',ThAct,'ElAct',ElAct,'method','parfit');
91 [a1,ad] = ltp_ifo2acc(ch1,ch2,pl);
92
93 %% plotting output psd
94
95 a1xx = a1.psd;
96 adxx = ad.psd;
97 iplot(a1xx,adxx)
98
99 %% plotting output lpsd
100
101 a1xx = a1.lpsd;
102 adxx = ad.lpsd;
103 iplot(a1xx,adxx,plist('LineWidths',{'All',3}))
104
105 %% Do step-by-step Checks
106
107 % frequency range for TFs calculation
108 f = [logspace(-6,log10(5),300)]';
109
110
111 %% Get intermediate steps
112
113 %%% IFO2POS conversion
114 ifod = S11*SDD - S1D*SD1;
115 x1 = (ch1.*SDD - ch2.*S1D)./ifod;
116 xd = (ch2.*S11 - ch1.*SD1)./ifod;
117
118
119 %%% Free Dynamics
120 sw1 = w1^2;
121 sw2 = w2^2;
122
123 D0COEFF = -3/35;
124 D1COEFF = -1/5;
125 D2COEFF = 2/7;
126
127 pldyn = plist('ALPHA0',sw1,'ALPHA1',0,'ALPHA2',1,...
128 'D0COEFF',D0COEFF,'D1COEFF',D1COEFF,'D2COEFF',D2COEFF,...
129 'X0',0,'TARGETUNITS',x1.yunits./(unit('s').^2));
130 a1fd = eqmotion(x1,pldyn);
131
132 % Total Differential force per unit of mass
133 pldyn = plist('ALPHA0',sw2,'ALPHA1',0,'ALPHA2',1,...
134 'D0COEFF',D0COEFF,'D1COEFF',D1COEFF,'D2COEFF',D2COEFF,...
135 'X0',0,'TARGETUNITS',x1.yunits./(unit('s').^2));
136 adfd = eqmotion(xd,pldyn);
137
138 % calculate dynamical cross talk
139 adct = x1.*(sw2 - sw1);
140 adct.setYunits(adfd.yunits);
141
142 % Total differential force
143 adfd = adfd + adct;
144
145 %%% DF Controller
146
147 parfrCdf = parfrac(Cdf);
148 filtCdf = miir(parfrCdf,plist('FS',fs));
149 ac1 = filter(ch1, filtCdf);
150
151 %%% Thrusters actuation
152 acomm1 = filter(ac1, ThAct);
153
154 %%% Sus controller
155 parfrCsus = parfrac(Csus);
156 filtCsus = miir(parfrCsus,plist('FS',fs));
157 acd = filter(ch2, filtCsus);
158
159 %%% Electrostatic actuation
160 acommd = filter(acd, ElAct);
161
162
163 %% Free Dynamics TFs ***
164
165 % get transfer functions from data
166 tffd1 = ltfe(ch1,a1fd);
167 tffd2 = ltfe(ch2,adfd);
168
169 % Get expected transfer functions
170 FD1 = ao(plist('fsfcn', '(2.*pi.*1i.*f).^2 + -1.3e-6', 'f', f));
171 FD2 = ao(plist('fsfcn', '-(-1e-4).*(2.*pi.*1i.*f).^2 + (-2e-6) - (-1.3e-6) - (-1e-4)*(-2e-6) + (2.*pi.*1i.*f).^2 + (-2e-6)', 'f', f));
172
173 % plot
174 plpl = plist('Legends', {'Dyn ch1','Dyn ch1 data','Dyn ch12','Dyn ch12 data'},...
175 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
176 iplot(abs(FD1),abs(tffd1),abs(FD2),abs(tffd2),plpl)
177
178 %% Controllers TFs ***
179
180 % get transfer functions from data
181 tffd1 = ltfe(ch1,ac1);
182 tffd2 = ltfe(ch2,acd);
183
184 % Get expected transfer functions
185 pl = plist('built-in','mdc1_Cdf','gain',1);
186 Cdf = rational(pl);
187 Cdfr = resp(Cdf,plist('f',f));
188
189 pl = plist('built-in','mdc1_Csus','gain',1);
190 Csus = rational(pl);
191 Csusr = resp(Csus,plist('f',f));
192
193 % plot
194 plpl = plist('Legends', {'Cdf','Cdf data','Csus','Csus data'},...
195 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
196 iplot(Cdfr,tffd1,Csusr,tffd2,plpl)
197
198 %% Actuation TFs
199
200 % get transfer functions from data
201 tffd1 = ltfe(ac1,acomm1);
202 tffd1.simplifyYunits;
203 tffd2 = ltfe(acd,acommd);
204 tffd2.simplifyYunits;
205
206 % TAU_th = 0.1; % s - time constant of thrusters actuation
207 % DELTAT_Th = 0.315; % s - delay introduced to command the force to thrusters system
208 %
209 % TAU_el = 0.01; % s - time constant of electrostatic actuation
210 % DELTAT_el = 0.305; % s - delay introduced to command the force to electrostatic actuators
211
212 % Get expected transfer functions
213 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));
214 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));
215
216 % plot
217 plpl = plist('Legends', {'Thrusters','Thrusters data','Electrostatic','Electrostatic data'},...
218 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
219 iplot(ThActr,tffd1,ElActr,tffd2,plpl)
220
221
222 %% Full model of acceleration noise
223
224 % get mdc1 ifo noise psd models
225 CSD = get_2D_test_obj_tf_psd(f);
226
227 % free dynamics
228 FD11 = ao(plist('fsfcn', '(2.*pi.*1i.*f).^2 + -1.3e-6', 'f', f));
229 FD21 = ao(plist('fsfcn', '-(-1e-4).*(2.*pi.*1i.*f).^2 + (-2e-6) - (-1.3e-6) - (-1e-4)*(-2e-6)', 'f', f));
230 FD22 = ao(plist('fsfcn', '-(-1e-4).*(2.*pi.*1i.*f).^2 + (-2e-6) - (-1.3e-6) - (-1e-4)*(-2e-6) + (2.*pi.*1i.*f).^2 + (-2e-6)', 'f', f));
231
232 Cdfr = resp(Cdf,plist('f',f));
233 Csusr = resp(Csus,plist('f',f));
234 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));
235 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));
236
237
238 TF11 = FD11 + Cdfr.*ThActr;
239
240 TF22 = FD22 + Csusr.*ElActr;
241
242 Acc1 = (2/fs).*TF11.*CSD(1,1).*conj(TF11);
243 Accd = (2/fs).*(FD21.*CSD(2,1).*conj(FD21) + TF22.*CSD(2,2).*conj(TF22));
244
245 % plot
246 plpl = plist('Legends', {'Ch1 data','Ch1','Ch12 data','Ch12'},...
247 'LineStyles', {'', '--', '', '--'},'LineWidths',{'All',3});
248 iplot(a1xx,abs(Acc1),adxx,abs(Accd),plpl)
249
250
251
252
253
254
255
256
257
258
259
260
261
262