close all
clear all
j = 7;

data = load('Plate Test MBARI Tank 01.DAT');  %%"Inch","lbf","PSI","PSI","Volt","degC","degC","cm"
dt = .1;


pos = data(:,1);
tension = data(:,2);


t = 0:dt:dt*(length(pos)-1);

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


tension = filtfilt(b,a,tension);



fh = figure('Position',[6 38 1267 690]);
subplot(2,1,1)
[ah1,h11,h12] = plotyy(t,pos,t,vel);
%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)');





Amp   = [36    36    36    36    36      36       36      18  18   18    18      18   18     18    18    18  18];
Per   = [50    25    17    14    12.5    11       10      50  25   17    12.5    10    8     7.14  6.25  5.5  5];
start = [1750  4000  5550  7250  9000   10400   11500   12800  15700  17300  18450  19200  20050  20850  21400  22700  23300  23850];
stop  = [3750  5200  6850  8350  9950   11250   12250   15600  16800  18050  19000  19650  20450  21250  21850  23050  23650  24150];

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

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

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

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

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

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

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


Force = W+(m+mu)*accel(idx_36)+Cd*.5*rho*Area*vel(idx_36).*abs(vel(idx_36));

subplot(2,1,2)
plot(t(idx_36),tension(idx_36),t(idx_36),Force/4.45);
xlabel('seconds');
ylabel('tension (lbs)');
legend('Measured Force','Estimated