clear all
close all

dpath = '134.89.12.189/data/';
dataDirs = dir([dpath 'R*']);

rho = 1025;  %kg/m^3
Area = 9.6;  %m^2
j = 0;
for i = 1:length(dataDirs)
  if(dataDirs(i).isdir == 1)
     ppath = [dpath dataDirs(i).name '/'];
     
     matfiles = dir([ppath 'D*.mat']);

     K = length(matfiles);
     for k = 1:K
       load([ppath matfiles(k).name]);
         j = j+1;  
         %MaxVolt(j) = max(RectifiedVoltage)
         FiltdownVelocity = filtfilt(ones(1,400)/400,1,downVelocity);
         downVelocity = downVelocity-FiltdownVelocity;
       
         PistonDistance = PistonDistance*.0254;  %Convert to meters
         ElongationRate = [diff(PistonDistance)./diff(time) 0];

         zdot = downVelocity+ElongationRate;
         zddot = [diff(zdot)./diff(time) 0];
       
         
         idx = find(ElongationRate > .1);
         
         N = length(idx);
         A = -ones(N,3);
         B = -4.45*LoadCellForce(idx)';
         A(:,1) = .5*rho*Area*zdot(idx)'.*abs(zdot(idx)');
         A(:,2) = zddot(idx)';
       
         x = A\B;
       
         disp(['Cd = ' num2str(x(1))]);
         disp(['virtual mass = ' num2str(x(2))]);
         disp(['Fmean = ' num2str(x(3))]);
       
       

         j = j+1;
         Cd(j) = x(1);
         vm(j) = x(2);
         Fmean(j) = x(3);
       
       
         close all
       
         figure
         ah1 = subplot(3,1,1);
         plot(time-time(1),downVelocity,time-time(1),ElongationRate,time-time(1),zdot)
         legend('downVel','ElongationRate','zdot');
         ah2 = subplot(3,1,2);
         %plot(time-time(1),-LoadCellForce,time-time(1),(-x(3)+x(1)*.5*rho*Area*zdot.*abs(zdot)+x(2)*zddot)/4.45);
         ah3 = subplot(3,1,3);
         plot(time(idx)-time(idx(1)),B,time(idx)-time(1),A*x);
         linkaxes([ah1 ah2 ah3],'x');
     
     end
  end
end