%% compare_exter_inter_pos.m
clear; close all
foldername = '/Volumes/Persistence Lab/Fault Prognostication/data_logs/Integration-test-2017-03-07';
f_num = 1;

%% prep external position data
type = 'external';
fnames = getTestbedDatafiles(foldername, type);
filename = [foldername filesep fnames{f_num}]; % last file on list
data = readinTestbedData(filename);

x = data;
x_pos = x.position.val; % x_pos = -(x_pos - mean(x_pos(1:100)));  % remove background (pos @ 0mm) and revers sign
x_vol = x.voltage.val;
x_cur = x.current.val;

x_p_time = datenum(datetime(x.position.time_lv,'ConvertFrom','posixtime'));
x_v_time = datenum(datetime(x.voltage.time_lv,'ConvertFrom','posixtime'));
x_c_time = datenum(datetime(x.current.time_lv,'ConvertFrom','posixtime'));

%% prep intrinsic data struct 
type = 'intrinsic';
fnames = getTestbedDatafiles(foldername, type);
filename = [foldername filesep fnames{f_num}]; % last file on list
data = readinTestbedData(filename);

y = data;

% put intrinsic data on the same time-grid
in_posi = interp1(y.position.time_lv, y.position.val, x.position.time_lv);
in_voli = interp1(y.voltage.time_lv, y.voltage.val, x.voltage.time_lv);
in_curi = interp1(y.current.time_lv, y.current.val, x.position.time_lv);

%% plot
figure
set(gcf, 'Position', [0 0 1 1])

ax = axes;
ax.Position = [0.1 0.7 0.8 0.25];
plot(ax, x_p_time, in_posi); hold on
plot(ax, x_p_time, x_pos)
title(fnames{f_num})
ylabel('Position (mm)')
legend('intrinsic','external')

ax1 = axes;
ax1.Position = [0.1 0.4 0.8 0.25];
plot(ax1, x_v_time, in_voli); hold on
plot(ax1, x_v_time, x_vol)
ylabel('Voltage (v)')
legend('intrinsic','external')

ax2 = axes;
ax2.Position = [0.1 0.1 0.8 0.25];
plot(ax2, x_c_time, in_curi); hold on
plot(ax2, x_c_time, x_cur)
ylabel('Current (mAmp)')
xlabel('Time (UTC)')
legend('intrinsic','external')

linkaxes([ax ax1 ax2],'x')
dynamicDateTicks([ax ax1 ax2], 'linked')
set([ax ax1], 'XTickLabel',[])

%% determine position signal slope deviation

[xb_p,xb_S,xb_mu] = polyfit(x_p_time, x_pos, 1);
[inb_p,inb_S,inb_mu] = polyfit(y.position.time_lv, y.position.val, 1);

div = 100 * abs(abs(inb_p(1)) - abs(xb_p(1))) / abs(inb_p(1));
disp(['Position signal slope deviation = ' num2str(div) '%.'])

