% Test ao/smallvectorfit functionality.
%
% A Monsky, M Hueller  29-03-10
%
% $Id: test_smallvectorfit.m,v 1.5 2010/04/01 08:29:21 anneke Exp $
%

% Parameters
datafiles = '2008-07-11_12-26';
test_data = sprintf('%s/example_file.dat', datafiles);
fs = 800000.0/24640.0;
%% Load data
% unfiltered Data
disp('*** Loading phase data....')
dpl  = plist('filename', test_data, 'fs', fs, 'columns', 1, 'yunits', '');
dpl2 = plist('filename', test_data, 'fs', fs, 'columns', 2, 'yunits', '');
dpl3 = plist('filename', test_data, 'fs', fs, 'columns', 3, 'yunits', '');

phir  = ao(dpl);
phi1  = ao(dpl2);
phi12 = ao(dpl3);
iplot(phir, phi12)

%% User input 
% starttime = input('enter start time: '); 
% endtime = input('enter end time: ');
%%
starttime = 40;
endtime   = 50;
%% Split and filter
pl = plist('type', 'lowpass', 'order', 1, 'gain', 1.0, 'fs', fs, 'fc', 1);
f = mfir(pl);
spl = plist('split_type', 'times', 'times',[starttime endtime]);
phi1_split = split(filter(phi1,f), spl);
phi12_split = split(filter(phi12,f), spl);
phir_split = split(filter(phir,f), spl);


%% Compute psi
psi1_split  = phi1_split - phir_split;
psi12_split = phi12_split - phir_split;
%% XYdata
xy_psi1  = ao(plist('xvals', phir_split.data.y, 'yvals', psi1_split.data.y, ...
  'xunits', phir_split.yunits, 'yunits', psi1_split.yunits));  % This could be replaced by ao/scatterData
xy_psi12 = ao(plist('xvals', phir_split.data.y, 'yvals', psi12_split.data.y, ...
  'xunits', phir_split.yunits, 'yunits', psi1_split.yunits)); % This could be replaced by ao/scatterData
%% correct for Doppler effect
phir_doppler  = dopplercorr(phir_split, plist('bin',50));
phi1_doppler  = dopplercorr(phi1_split, plist('bin',50));
phi12_doppler = dopplercorr(phi12_split, plist('bin',50));
%% Psi
psi1_doppler = phi1_doppler - phir_doppler;
psi12_doppler = phi12_doppler - phir_doppler;
psi12_doppler.data.growT;
%% XYdata
xy_psi1_doppler = ao(plist('xvals', phir_doppler.data.y, 'yvals', psi1_doppler.data.y, ...
  'xunits', phir_doppler.yunits, 'yunits', psi1_doppler.yunits));  % This could be replaced by ao/scatterData
xy_psi12_doppler = ao(plist('xvals', phir_doppler.data.y, 'yvals', psi12_doppler.data.y, ...
  'xunits', phir_doppler.yunits, 'yunits', psi12_doppler.yunits));  % This could be replaced by ao/scatterData
iplot(xy_psi1, xy_psi1_doppler)
%% Fit
[fit_psi1 fit_psi12] = smallvectorfit(phir_doppler, phi1_doppler, psi1_doppler, phi12_doppler, psi12_doppler);
%% Linear combination
[svec_psi1 svec_psi12] = smallvector_lincom(...
  phir_doppler, phi1_doppler, psi1_doppler, phi12_doppler, psi12_doppler, fit_psi1, fit_psi12);

%% create an xy data object from the linear combination
xy_klvek_psi1 = scatterData(phir_doppler,svec_psi1);
xy_klvek_psi12 = scatterData(phir_doppler,svec_psi12);
%% Plot data and fit together
iplot(xy_psi1_doppler, xy_klvek_psi1)
iplot(xy_psi12_doppler, xy_klvek_psi12)
%% read out small vector coefficiients 
fit_klvek_psi1 = find(fit_psi1.procinfo,'coef_psi1');
fit_klvek_psi12 = find(fit_psi12.procinfo,'coef_psi12');