close all
clear all


data = load('wPlates-A36-sin.csv');  %%"Inch","lbf","PSI","PSI","Volt","degC","degC","cm"
dt = .04;

t = data(:,1);
pos = data(:,2);
force = data(:,3);


clear data;

vel = [0; diff(pos)/dt];
idx = find(vel > 25); %Remove outliers due to starting and stopping logging.
vel(idx) = 0;
%Filter numerical derivatives to smooth.
b = ones(1,8)/8;
a = [1];
vel = filtfilt(b,a,vel);

accel = [0; diff(vel)/dt];
idx = find(accel > 30); %Remove outliers due to starting and stopping logging.
accel(idx) = 0;
accel = filtfilt(b,a,accel);


force = filtfilt(b,a,force);
force = force-180;  %Correct for weight of ram.




fh = figure('Position',[6 38 1267 690]);
ah1 =subplot(3,1,1)
plot(pos)
ylabel('pos (in)');
ah2 = subplot(3,1,2)
plot(vel);
ylabel('vel (in/s)');
ah3 = subplot(3,1,3)
plot(force)
ylabel('force(lbs)');
xlabel('time (s)');
linkaxes([ah1 ah2 ah3],'x');



      
start = [243 4100 9204 12480];  
stop  = [2867 6005 10160 13960];   

idx = [];
for i = 1:length(start)
    idx = [idx start(i):stop(i)];
end

%idx = start(j):stop(j);

pos = pos*.0254;  %Convert to meters.
vel = vel*.0254;
accel = accel*.0254;

W = mean(force(idx))*4.45;  %Wet weight in newtons
m = 2000/2.2;  %Mass in kg

Area = 14.34;  %plate area in m^2
rho = 1025;  %water density.

F = force(idx)*4.45;  %newtons
A = [accel(idx) .5*rho*Area*vel(idx).*abs(vel(idx))];
B = W+m*accel(idx);

x = A\(F-B)
mu = x(1);
Cd = x(2);


Force = W+(m+mu)*accel(idx)+Cd*.5*rho*Area*vel(idx).*abs(vel(idx));
fh = figure('Position',[6 38 1267 690]);

subplot(2,1,1)
[ah1,h11,h12] = plotyy(t,pos/.0254,t,vel/.0254);
%set(ah1(1),'YLim',[-1 1]);
%set(ah1(1),'YTick',-1:2/8:1);
set(ah1(2),'YLim',[-25 25]);
set(ah1(2),'YTick',-25:50/10:25);
set(get(ah1(1),'Ylabel'),'String','Position (in)');
set(get(ah1(2),'Ylabel'),'String','Vel (in/s)');

ah2 = subplot(2,1,2)
plot(t(idx),force(idx),t(idx),Force/4.45);
xlabel('seconds');
ylabel('force (lbs)');
legend('Measured Force',['Estimated Force: mu = ' num2str(mu,4) 'kg, Cd = ' num2str(Cd,2)],'Location','NorthWest');

linkaxes([ah1(1) ah1(2) ah2],'x'); 