close all
clear all

MeasuredResultsPath = '.\MeasuredResults\';
WaveDataResultsPath = '..\CDIPWaveRiderData\';

PlotDirectory = '.\SeaStateComparisonPlots\';

files = dir([WaveDataResultsPath 'WaveData*.mat']);
  %Put in chronological order.
  for i = 1:length(files)
    load([WaveDataResultsPath files(i).name])
    SerialDate(i) = sdate;
    BL(i) = BudalLimit;
    Hs(i) = SigWaveHeight;
    Tp(i) = PeakPeriod;  
  
    datestr(sdate)
  end
  [Dates IX] = sort(SerialDate);
  files = files(IX);
  SerialDate = SerialDate(IX);
  BL = BL(IX);
  Hs = Hs(IX);
  Tp = Tp(IX);
  
  clear  Dates IX i;


  for file_num = 1:length(files)
    wavedata_files{file_num} = files(file_num).name;
  end


Meas_meanElecPow = NaN(1,length(files));
Meas_PercentofBudal = NaN(1,length(files));
Meas_StdPTOElongation = NaN(1,length(files));
Meas_MaxPTOElongation(i) = NaN(1,length(files));
Meas_MinPTOElongation(i) = NaN(1,length(files));
Meas_StdPTOElongation(i) = NaN(1,length(files));

Meas_StdPTOElongationRate = NaN(1,length(files));
Meas_MaxPTOElongationRate(i) = NaN(1,length(files));
Meas_MinPTOElongationRate(i) = NaN(1,length(files));
Meas_StdPTOElongationRate(i) = NaN(1,length(files));

plot_num = 0;
for i = 1:length(files);
       
  n = strfind(files(i).name,'WaveData');
  wavefilenm = files(i).name(n:end)
    
  if(exist([MeasuredResultsPath 'MeasuredResults_for' wavefilenm]))
     MeasResults = LoadMeasResults([MeasuredResultsPath 'MeasuredResults_for' wavefilenm]);
     Meas_meanElecPow(i) = MeasResults.meanElecPow;
     Meas_PercentOfBudal(i) = Meas_meanElecPow(i)/BL(i);
     
     Meas_MeanPTOElongation(i) = MeasResults.MeanPTOElongation;
     Meas_MaxPTOElongation(i) = MeasResults.MaxPTOElongation;
     Meas_MinPTOElongation(i) = MeasResults.MinPTOElongation;
     Meas_StdPTOElongation(i) = MeasResults.StdPTOElongation;
              
     Meas_StdPTOElongationRate(i) = MeasResults.StdPTOElongationRate; 
     Meas_MaxPTOElongationRate(i) = MeasResults.MaxPTOElongationRate;
     Meas_MinPTOElongationRate(i) = MeasResults.MinPTOElongationRate;
     Meas_StdPTOElongationRate(i) = MeasResults.StdPTOElongationRate;
  end
  
  
  
end

idx = find(~isnan(Meas_meanElecPow));


%Plot figures for mean power versus Hs
figure
[AX,H1,H2] = plotyy(SerialDate(idx),Meas_meanElecPow(idx),SerialDate(idx),Hs(idx))
%legend(['Mean Power = ' num2str(nanmean(Meas_meanElecPow(36:602)),3) ' W']);
set(gca,'FontSize',14)


set(AX(1),'Ylim',[0 1000]);
set(AX(2),'Ylim',[0 2]);

set(get(AX(1),'Ylabel'),'String','Rectified Power (W)') 
set(get(AX(1),'Ylabel'),'FontSize',14) 
set(AX(1),'FontSize',14)
set(get(AX(2),'Ylabel'),'String','Significant Wave Height (m)') 
set(get(AX(2),'Ylabel'),'FontSize',14) 
set(AX(2),'FontSize',14)
datetick(AX(1),'x','mmm-dd','keeplimits','keepticks') 
set(AX(2),'XTick',[])

print('-dpdf',[PlotDirectory 'RectifiedPowerAndHs']);
print('-dpng',[PlotDirectory 'RectifiedPowerAndHs']);

figure
plot(Hs(idx),Meas_meanElecPow(idx),'x')
set(gca,'FontSize',14)
ylabel('Rectified Electrical Power');
xlabel('Significant Wave Height (m)');

print('-dpdf',[PlotDirectory 'RectifiedPowerVsHs']);
print('-dpng',[PlotDirectory 'RectifiedPowerVsHs']);


%Plot figures for Budal percentage versus Hs
figure
[AX,H1,H2] = plotyy(SerialDate(idx),100*Meas_PercentOfBudal(idx),SerialDate(idx),Hs(idx))
set(gca,'FontSize',14)

set(AX(1),'Ylim',[0 14]);
set(AX(2),'Ylim',[0 2]);

set(get(AX(1),'Ylabel'),'String','% of Budal Limit Captured') 
set(get(AX(1),'Ylabel'),'FontSize',14) 
set(AX(1),'FontSize',14)
set(get(AX(2),'Ylabel'),'String','Significant Wave Height (m)') 
set(get(AX(2),'Ylabel'),'FontSize',14) 
set(AX(2),'FontSize',14)
datetick(AX(1),'x','mmm-dd','keeplimits','keepticks') 
set(AX(2),'XTick',[])

print('-dpdf',[PlotDirectory 'BudalPercentageAndHs']);
print('-dpng',[PlotDirectory 'BudalPercentageAndHs']);

figure
plot(Hs(idx),100*Meas_PercentOfBudal(idx),'x')
set(gca,'FontSize',14)
ylabel('% of Budal Limit Captured');
xlabel('Significant Wave Height (m)');

print('-dpdf',[PlotDirectory 'BudalPercentagePowerVsHs']);
print('-dpng',[PlotDirectory 'BudalPercentagePowerVsHs']);



%Plot figures for Ram behavior
figure
[AX,H1,H2] = plotyy(SerialDate(idx),[],SerialDate(idx),Hs(idx))
set(gca,'FontSize',14)

set(AX(1),'Ylim',[0 14]);
set(AX(2),'Ylim',[0 2]);

set(get(AX(1),'Ylabel'),'String','% of Budal Limit Captured') 
set(get(AX(1),'Ylabel'),'FontSize',14) 
set(AX(1),'FontSize',14)
set(get(AX(2),'Ylabel'),'String','Significant Wave Height (m)') 
set(get(AX(2),'Ylabel'),'FontSize',14) 
set(AX(2),'FontSize',14)
datetick(AX(1),'x','mmm-dd','keeplimits','keepticks') 
set(AX(2),'XTick',[])

print('-dpdf',[PlotDirectory 'BudalPercentageAndHs']);
print('-dpng',[PlotDirectory 'BudalPercentageAndHs']);

figure
plot(Hs(idx),100*Meas_PercentOfBudal(idx),'x')
set(gca,'FontSize',14)
ylabel('% of Budal Limit Captured');
xlabel('Significant Wave Height (m)');

print('-dpdf',[PlotDirectory 'BudalPercentagePowerVsHs']);
print('-dpng',[PlotDirectory 'BudalPercentagePowerVsHs']);