close all
clear all


data = load('PlateinTestTank_2012-8-3-CSwoPlates.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 = [1617 2648 3557 4350 5483 5941 7425 8250 8932 9486 10090 10670 12520 12930 13440 14290 17680 18110 18740 19260];  
stop  = [1992 3048 3691 4522 5594 6049 7508 8345 9020 9565 10140 10740 12850 12990 13500 14350 17810 18260 19000 19530];

g=9.81;
m = 2100/2.2;  %Mass in kg

Area = 9.63;  %plate area in m^2
rho = 1025;  %water density.
V=.204; %volume in m^3

idx = [];
for i = 1:length(start)
    idx = [idx start(i):stop(i)];
    meanspeed(i)=mean(vel(idx));
    meanforce(i)=mean(force(idx));
    Cd(i)=(meanforce(i) +rho*V*g - m*g)/(.5*rho*Area*(meanspeed(i))^2);
end
%idx = start(j):stop(j);
meanspeed
meanforce
Cd
pos = pos*.0254;  %Convert to meters.
vel = meanspeed*.0254;
accel = accel*.0254;

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'); 