close all
clear all


data = load('wPlates-CS.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.
idx2 = find(vel < -25);
vel(idx) = 0;
vel(idx2) = 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 = [2 520 1325 2375 2993 3586 4193 4621 5174 5561 6526 7014 8225 8814 9587 10630 12120 13060 14750];  
stop  = [118 872 1565 2588 3160 3719 4345 4756 5294 5721 6662 7153 8345 8936 9695 10750 12230 13180 14870];

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;

V = 2250;
m = 2000/2.2; %Mass in kg
g = 9.81;

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

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

x = A\B;
Cd = x(1);


Force = m*g + Cd*.5*rho*Area*vel(idx).*abs(vel(idx)) - rho*V*g;
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: Cd = ' num2str(Cd,2)],'Location','NorthWest');

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