close all
clear all

WaveDataPath = './CDIPWaveRiderData/';


WaveDataFiles = dir([WaveDataPath '*.mat']);

n_WaveDataFiles = length(WaveDataFiles);

Hs = NaN(n_WaveDataFiles,1);
Tp = NaN(n_WaveDataFiles,1);
StartTime = NaN(n_WaveDataFiles,1);
EndTime = NaN(n_WaveDataFiles,1);

RectPower = NaN(n_WaveDataFiles,1);

f = cell(n_WaveDataFiles,1);
Pxx = cell(n_WaveDataFiles,1);

%Read in all wave-data
for i = 1:n_WaveDataFiles
    disp(num2str(i));
    WaveData = load([WaveDataPath WaveDataFiles(i).name]);
    Hs(i) = WaveData.SigWaveHeight;
    Tp(i) = WaveData.PeakPeriod;
    StartTime(i) =  WaveData.sdate;
    EndTime(i) = StartTime(i)+1/24/2;
    f{i} = WaveData.f;
    Pxx{i} = WaveData.PxxHeave;    
end



%Now work through all measured data and make statistics of data
%coorelating to timeperiods of each wave-data spectra integration time.
%This is done bit by bit to avoid loading the entire dataset into memory.



%Assembled all filenames of results in chronological order
results_path = './DecodedResults/'
dataDirs = dir([results_path 'R*']);
i = 0;
for k = 1:1:length(dataDirs)
  if(dataDirs(k).isdir == 1)
     dpath = [results_path dataDirs(k).name '/'];
     files = dir([dpath 'D*']);
     for j = 1:length(files)
       i = i + 1;
       data_files{i,1} = [results_path  dataDirs(k).name '/' files(j).name]
     end   
  end
end


i = 1;
meas_data = load(data_files{i});
meas_date = meas_data.serialdate;
meas_power = meas_data.RectifiedVoltage.*meas_data.RectifiedCurrent;

time_epsilon = .05/(24*60);  %3 seconds

j = 1;
while(j < n_WaveDataFiles)

idx = find((meas_date < EndTime(j)) & (meas_date > StartTime(j)));

if(~isempty(idx))  %Got hits.
  disp(['j = ' num2str(j) '   length idx = ' num2str(length(idx))]);


  StartDiff = StartTime(j)-meas_date(idx(1));
  EndDiff = EndTime(j)-meas_date(idx(end));

  disp(['StartTime = ' datestr(StartTime(j)) ' meas_date(idx(1)) = ' datestr( meas_date(idx(1)))]);
  disp(['EndTime = ' datestr(EndTime(j)) ' meas_date(idx(end)) = ' datestr( meas_date(idx(end)))]);

  disp(['Start time diff = ' num2str(StartDiff*24*60) 'min   End time diff = ' num2str(EndDiff*24*60) 'min ']);

  %if(length(idx) >= 15000)
  
else  %No hits in this data so fake the start/end times so the next file gets loaded
   StartDiff = time_epsilon+1;    
end

if((abs(StartDiff) < time_epsilon) & (abs(EndDiff) < time_epsilon))
%Have a complete half hour of data (or close anyway), compute statistics for this time-window
  %and remove data since it won't be needed again.
  RectPower(j) = mean(meas_power(idx));
  j = j+1;  %Got this one so increase to next wave-data date range
  meas_date = meas_date(idx(end):end);
  meas_power = meas_power(idx(end):end);  

else
  %Don't have a complete half-hour of data, load the next data file and concatenate
  i = i+1;
  data_files{i}
  meas_data = load(data_files{i});
  meas_date = [meas_date meas_data.serialdate];
  meas_power = [meas_power meas_data.RectifiedVoltage.*meas_data.RectifiedCurrent];
  
  %If the first date of the new file is after the needed date to complete
  %the current wave-data time window, this means the system stopped and that time doesn't exist,
  % so just increase j and that wave-data window will remain empty.
  if(meas_data.serialdate(1) > EndTime(j))
      j = j + 1;
      disp('Increasing j without finding complete data window');
     % pause
  end
end
  
    
end

   

%Remove wave-data periods in which power data wasn't collected

idx = find(~isnan(RectPower));

RectPower = RectPower(idx);
Hs = Hs(idx);
Tp = Tp(idx);
StartTime = StartTime(idx);
EndTime = EndTime(idx);

f = f(idx);
Sy = Pxx(idx);

figure
plot(Hs,RectPower);


save('SpectrasAndPowers.mat','RectPower','Hs','Tp','StartTime','EndTime','f','Sy');








