close all
clear all


DecodedResultsPath = '..\DecodedResults\';

directories = dir(DecodedResultsPath);

PistonVel = [];
RectifiedVoltage = [];
RectifiedCurrent = [];

max_speed = 2*39.4;
max_power = 12000;
n_speeds = 100;
n_powers = 100;
n_hist = zeros(n_speeds,n_powers);
NN = 0;

AllSpeeds = [];
AllPowers = [];

energy_total = 0;

[X,Y] = meshgrid(0:max_speed/(n_speeds-1):max_speed,0:max_power/(n_powers-1):max_power);

edges{1} = X(1,:);
edges{2} = Y(:,1);
kk = 0;

dt = 0.1;

for i = 3:length(directories)
    data_dir = [DecodedResultsPath directories(i).name];
    
    if(isdir(data_dir))
       files = dir(data_dir);       
       for j = 3:length(files)
          data = load([data_dir '\' files(j).name]);
          files(j).name
          if((data.serialdate(1) > datenum('03-Sep-2013 21:34:15')) && ...
             (data.serialdate(1) < datenum('08-Oct-2013 07:49:30')))
            datestr(data.serialdate(1))
            PistonVel = [0 diff(data.PistonDistance)./dt];
            PistonSpd = abs(PistonVel);
            RectifiedVoltage = data.RectifiedVoltage;
            RectifiedCurrent = data.RectifiedCurrent;
            RectifiedPower = RectifiedVoltage.*RectifiedCurrent;
            LoadCellForce = data.LoadCellForce;
            BuoyDepthPressure = data.BuoyDepthPressure;
            dat = [PistonSpd' RectifiedPower'];
            kk=kk+1;
            PistonTravel(kk) = sum(PistonSpd)*dt;
            mean_LoadCell_log(kk) = mean(LoadCellForce)
            max_LoadCell_log(kk) = max(LoadCellForce);
            mean_speed_log(kk) = mean(dat(:,1))
            max_speed_log(kk) = max(dat(:,1));
            max_power_log(kk) = max(dat(:,2));
            SpeedBins = 0:.01:2/0.0254;
            SpeedHist{kk} = histc(dat(:,1),SpeedBins);
            ForceBins = 0:10000/100:10000;
            ForceHist{kk} = histc(LoadCellForce,ForceBins);
            energy_total = sum([energy_total ; dat(:,2)]);
            n_hist = n_hist+hist3(dat,'Edges',edges);
            NN = NN+length(PistonSpd);
            
            AllSpeeds = [AllSpeeds PistonSpd];
            AllPowers = [AllPowers RectifiedPower];
            
            
          end
          clear data;
       end
       
    end
end
    


figure
h = bar3(n_hist/NN);
colormap('Hot');
for n=1:numel(h)
    cdata=get(h(n),'zdata');
    set(h(n),'cdata',cdata,'facecolor','interp')
end
colorbar
colormap(flipud(colormap))

figure
contour(X',Y',n_hist,100);
title('Contour plot of percentage of time spent at each power level');
xlabel('Speed (in/s)');
ylabel('Power (W)');
colorbar



e_hist = (1/energy_total)*n_hist.*Y';
figure
contour(X',Y',e_hist,20);
title('Contour plot of percentage of energy captured at each power level');
xlabel('Speed (in/s)');
ylabel('Power (W)');
colorbar

Ap = 0.25*pi*(2^2-1.375^2);  %in^2
RPM = 0.95*X*Ap*2.54^3*60/4.9;  %4.9cm^3/rev  Factor of 0.95 is a wag at the volumetric efficiency of hydraulic motor
max_rpm_log = 0.95*max_speed_log*Ap*2.54^3*60/4.9; 
figure
[C,ch] = contour(RPM',Y',e_hist,10);
title('Contour plot of percentage of energy captured at each power level');
xlabel('Speed (RPM)');
ylabel('Power (W)');
xlim([0 max(max_rpm_log)]);
ylim([0 max(max_power_log)]);
colorbar
%set(ch,'ShowText','on','TextStep',get(ch,'LevelStep')*2)




figure('Position',[520   173   766   925],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
plot(1:kk,mean_LoadCell_log);
xlabel('hour');
ylabel('Mean Load Cell Force (lbs)');
subplot(2,1,2)
plot(1:kk,max_LoadCell_log);
xlabel('hour');
ylabel('Max Load Cell Force (lbs)');


figure('Position',[520   173   766   925],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
plot(1:kk,mean_speed_log*.0254,1:kk,PistonTravel*.0254/3600);
xlabel('hour');
ylabel('Mean Piston Speed (m/s)');
subplot(2,1,2)
plot(1:kk,max_speed_log*.0254);
xlabel('hour');
ylabel('Mean Piston Speed (m/s)');

TotalForceHist = zeros(1,length(ForceHist{1}));
for k = 1:length(ForceHist)
  TotalForceHist = TotalForceHist + ForceHist{k};
end

TotalSpeedHist = zeros(length(SpeedHist{1}),1);
for k = 1:length(ForceHist)
  TotalSpeedHist = TotalSpeedHist + SpeedHist{k};
end


figure('Position',[520   173   766   925],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
bar(ForceBins,TotalForceHist)
subplot(2,1,2)
bar(SpeedBins,TotalSpeedHist)
xlabel('speed (in/s)')
ylabel('Occurences');

delta = 0.3;  %in^3/s
Ap = 1.37;  %in^2
figure('Position',[520   173   766   925],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
bar(SpeedBins*Ap*60/delta,TotalSpeedHist)
xlabel('speed (RPM)')
ylabel('Occurences');

AllRPMs = AllSpeeds*60*Ap/delta;


figure('Position',[520   173   766   925],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
plot(AllRPMs,AllPowers,'x')
xlabel('RPM')
ylabel('Watts')

RPMBins = 0:250:10000;
[bincounts,ind] = histc(AllRPMs,RPMBins);

for i = 1:length(RPMBins)
    idx = find(ind == i);
    
    Joules(i) = sum(AllPowers(idx))*.1;
end

Total_Energy = sum(Joules);

figure('Position',[520   173   766   925],'PaperOrientation','portrait','PaperPosition',[0 0 8.5 11]);
subplot(2,1,1)
bar(RPMBins,100*Joules/Total_Energy);
xlabel('RPM');
ylabel('Percentage of Total Energy');
title('Energy Distribution during 1 month deployment');
subplot(2,1,2)
plot(RPMBins,cumsum(100*Joules/Total_Energy))
xlabel('RPM');
ylabel('Cumulative percentage');
grid on





