Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/eqmotion.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 % EQMOTION solves numerically a given linear equation of motion
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: EQMOTION solves numerically a given linear equation
5 % of motion:
6 % d^2 x dx
7 % F(t) = alpha2 ------- + alpha1 ------ + alpha0 (x-x0)
8 % dt^2 dt
9 %
10 % CALL: eqmotion(a)
11 % b = eqmotion(a,pl)
12 %
13 % INPUTS: a - analysis object(s) containing data as a function of
14 % time.
15 % pl - parameter list containing input parameters.
16 %
17 % OUTPUTS: b - analysis object(s) containing output data as a function
18 % of time.
19 %
20 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'eqmotion')">Parameters Description</a>
21 %
22 % NOTE: Derivative estimation is performed with the parabolic fit
23 % approximation by default [1, 2]. Try to change D#COEFF to use another
24 % method. D0COEFF is used to calculate a five point data smoother to be
25 % applied to the third term at the second member of the equation above. If
26 % you do not whant to smooth data (before the multiplication with alpha0)
27 % you have to input NaN for D0COEFF.
28 % See also help for ao/diff and utils.math.fpsder.
29 %
30 % REFERENCES:
31 % [1] L. Ferraioli, M. Hueller and S. Vitale, Discrete derivative
32 % estimation in LISA Pathfinder data reduction, Class. Quantum Grav.,
33 % 7th LISA Symposium special issue.
34 % [2] L. Ferraioli, M. Hueller and S. Vitale, Discrete derivative
35 % estimation in LISA Pathfinder data reduction
36 % http://arxiv.org/abs/0903.0324v1
37 %
38 % VERSION: $Id: eqmotion.m,v 1.13 2011/04/11 10:24:45 mauro Exp $
39 %
40 % SEE ALSO: ao/diff, utils.math.fpsder
41 %
42 %
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44
45 function varargout = eqmotion(varargin)
46
47 % Check if the method was called by another method
48 callerIsMethod = utils.helper.callerIsMethod;
49
50 %%% Check if this is a call for parameters
51 if utils.helper.isinfocall(varargin{:})
52 varargout{1} = getInfo(varargin{3});
53 return
54 end
55
56 %%% Collect input variable names
57 in_names = cell(size(varargin));
58 for ii = 1:nargin,in_names{ii} = inputname(ii);end
59
60 %%% Collect all AOs
61 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
62 pli = utils.helper.collect_objects(varargin(:), 'plist', in_names);
63
64 %%% Decide on a deep copy or a modify
65 %%% REMARK: If you create a new AO (call the constructor) then
66 %%% it is not necessay to copy the input-AOs !!!!!!!!!!!!!!!!!!!!!!!!!
67 bs = copy(as, nargout);
68
69 %%% Combine plists
70 pl = combine(pli, getDefaultPlist);
71
72 %%% Get Parameters
73 alpha0 = find(pl,'ALPHA0');
74 alpha1 = find(pl,'ALPHA1');
75 alpha2 = find(pl,'ALPHA2');
76 X0 = find(pl,'X0');
77 d0c = find(pl,'D0COEFF');
78 d1c = find(pl,'D1COEFF');
79 d2c = find(pl,'D2COEFF');
80 tunits = find(pl,'TARGETUNITS');
81
82 % check if the params are AOs
83 if ~isa(tunits,'unit')
84 tunits = unit(tunits);
85 end
86 if ~isa(alpha0,'ao')
87 alpha0 = cdata(alpha0);
88 alpha0.setYunits(tunits./unit(as.yunits));
89 alpha0 = ao(alpha0);
90 alpha0.simplifyYunits;
91 end
92 if ~isa(alpha1,'ao')
93 alpha1 = cdata(alpha1);
94 alpha1.setYunits(tunits.*unit('s')./unit(as.yunits));
95 alpha1 = ao(alpha1);
96 alpha1.simplifyYunits;
97 end
98 if ~isa(alpha2,'ao')
99 alpha2 = cdata(alpha2);
100 alpha2.setYunits(tunits.*(unit('s').^2)./unit(as.yunits));
101 alpha2 = ao(alpha2);
102 alpha2.simplifyYunits;
103 end
104 if ~isa(X0,'ao')
105 if isempty(X0)
106 X0 = cdata(0);
107 X0.setYunits(as.yunits);
108 X0 = ao(X0);
109 else
110 X0 = cdata(X0);
111 X0.setYunits(as.yunits);
112 X0 = ao(X0);
113 end
114 end
115 if isa(d0c,'ao')
116 d0c = d0c.data.y;
117 end
118 if isa(d1c,'ao')
119 d1c = d1c.data.y;
120 end
121 if isa(d2c,'ao')
122 d2c = d2c.data.y;
123 end
124
125 %%% go through analysis objects
126 for kk = 1:numel(bs)
127
128 %%% Calculate derivatives
129 if ~isnan(d0c) % do the smoothing
130 a0 = diff(bs(kk),plist('method', 'FPS', 'ORDER', 'ZERO', 'COEFF', d0c));
131 else
132 a0 = copy(bs(kk),1); % just use input data as they are
133 end
134 a1 = diff(bs(kk),plist('method', 'FPS', 'ORDER', 'FIRST', 'COEFF', d1c));
135 a2 = diff(bs(kk),plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', d2c));
136
137 %%% Calculate Force
138 b0 = (a0 - X0);
139 b0 = b0*alpha0;
140 b1 = a1*alpha1;
141 b2 = a2*alpha2;
142 bs(kk) = b2 + b1 + b0;
143 % simplify units
144 bs(kk).simplifyYunits(plist('prefixes', false));
145
146 %%% Set Name
147 bs(kk).name = sprintf('eqmotion(%s)', ao_invars{kk});
148
149 if ~callerIsMethod
150 %%% Set Name
151 bs(kk).name = sprintf('eqmotion(%s)', ao_invars{kk});
152 %%% Add History
153 bs(kk).addHistory(getInfo('None'), pl, ao_invars(kk), [as.hist(kk)]);
154 end
155
156
157 end
158
159 %%% Output
160 if nargout == numel(bs)
161 % List of outputs
162 for ii = 1:numel(bs)
163 varargout{ii} = bs(ii);
164 end
165 else
166 % Single output
167 varargout{1} = bs;
168 end
169
170 end
171
172 %--------------------------------------------------------------------------
173 % Get Info Object
174 %--------------------------------------------------------------------------
175 function ii = getInfo(varargin)
176 if nargin == 1 && strcmpi(varargin{1}, 'None')
177 sets = {};
178 pl = [];
179 else
180 sets = {'Default'};
181 pl = getDefaultPlist;
182 end
183 % Build info object
184 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: eqmotion.m,v 1.13 2011/04/11 10:24:45 mauro Exp $', sets, pl);
185 end
186
187 %--------------------------------------------------------------------------
188 % Get Default Plist
189 %--------------------------------------------------------------------------
190 function plout = getDefaultPlist()
191 persistent pl;
192 if exist('pl', 'var')==0 || isempty(pl)
193 pl = buildplist();
194 end
195 plout = pl;
196 end
197
198 function pl = buildplist()
199 pl = plist();
200
201 % ALPHA0
202 p = param({'ALPHA0','Zero order coefficient. Input a cdata ao with the proper units or a number.'}, ...
203 {1, {0}, paramValue.OPTIONAL});
204 pl.append(p);
205
206 % ALPHA1
207 p = param({'ALPHA1','First order coefficient. Input a cdata ao with the proper units or a number.'},...
208 {1, {0}, paramValue.OPTIONAL});
209 pl.append(p);
210
211 % ALPHA2
212 p = param({'ALPHA2','Second order coefficient. Input a cdata ao with the proper units or a number.'}, ...
213 {1, {0}, paramValue.OPTIONAL});
214 pl.append(p);
215
216 % X0
217 p = param({'X0','Data offset. Input a cdata ao with the proper units or a number.'}, ...
218 {1, {0}, paramValue.OPTIONAL});
219 pl.append(p);
220
221 % D0COEFF
222 p = param({'D0COEFF','Data smoother coefficient.'}, ...
223 {1, {-3/35}, paramValue.OPTIONAL});
224 pl.append(p);
225
226 % D1COEFF
227 p = param({'D1COEFF','First derivative coefficient.'}, ...
228 {1, {-1/5}, paramValue.OPTIONAL});
229 pl.append(p);
230
231 % D2COEFF
232 p = param({'D2COEFF','Second derivative coefficient.'}, ...
233 {1, {2/7}, paramValue.OPTIONAL});
234 pl.append(p);
235
236 % Target units
237 p = param({'TARGETUNITS','Set this parameter if you input just numbers for the ALPHA# coefficients.'}, ...
238 {1, {'N'}, paramValue.OPTIONAL});
239 pl.append(p);
240
241 end
242 % END
243