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 = [8 525 1328 2389 2993 3588 4195 4632 5170 5597 6545 8248 9606 12150];  
stop  = [114 861 1563 2586 3158 3719 4341 4752 5292 5690 6648 8329 9680 12210];

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

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

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

idx = [];
meanspeed = [];
meanforce = [];
meanaccel = [];

for i = 1:length(start)
    idx = [idx start(i):stop(i)];
    meanspeed = [meanspeed mean(vel(start(i):stop(i)))];
    meanforce = [meanforce mean(force(start(i):stop(i)))];
    meanaccel = [meanaccel mean(accel(start(i):stop(i)))];
    if meanspeed(i)>0
        Cd(i) = (meanforce(i) + rho*V*g - m*g)/(.5*rho*Area*(meanspeed(i)^2));
    else
        Cd(i) = (m*g - meanforce(i) - rho*V*g)/(.5*rho*Area*(meanspeed(i)^2));
    end
    
end