close all
clear all

V = (5000/2.2)/1025; 

filename = '156Data';

g = 9.81;
ff = 0:.002:.58;
ff = ff';
alpha = 8.1e-3;
beta = .74;
w = 2*pi*ff;
df = ff(2)-ff(1);
 Hs = 1:2;
 clear S_Hs;
 for i = 1:length(Hs)
   U19 = sqrt(Hs(i)*g/0.21);
   w0 = g/U19;
   S_Hs(i,:) = 2*pi*(alpha*g^2)*exp(-beta*(w0./w).^4)./(w.^5);
   S_Hs(i,1) = S_Hs(i,2);
 
 
   Pb_PM(i) = BudalLimit(S_Hs(i,:)',ff,V);
 end

 


cmd = (['wc ' filename]);

[status result] = system(cmd);

if(status > 0) 
  error([cmd ' returns error']);
end

n_lines = str2num(strtok(result));
n_LinesPerRecord = 74;
n_Records = n_lines/n_LinesPerRecord;

fid = fopen(filename);

%tline = cell{n_LinesPerRecord,1};

for i = 1:n_Records
    
  for j = 1:n_LinesPerRecord
    tline{j,1} = fgetl(fid); 
  end
  disp(['Processing ' tline{1}]);
  [serialdate(i) freq{i} bw{i} energy{i} Hs(i) Tp(i)] = ParseDataWellSpectrum(tline);
  
  xi = freq{i}(1):.005:freq{i}(end);
  yi_energy = interp1(freq{i},energy{i},xi);
  
  Pb_meas(i) = BudalLimit(yi_energy,xi,V);
  
end

fclose(fid);

load('../History.mat');


for i = 1:4:n_Records-4
  figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
  for j = 0:3
    if ((i+j) <= n_Records)
      subplot(2,2,j+1)
      plot(ff,S_Hs,freq{i+j},energy{i+j});
      xlabel('freq (Hz)');
      ylabel('Energy (m^2/Hz)');
      title(['Wave Conditions: ' datestr(serialdate(i+j))]);
      legend(['Hs = 1m  (P_b = ' num2str(Pb_PM(1),4) ' W)'],['Hs = 2m  (P_b = ' num2str(Pb_PM(2),4) ' W)'],['Measured (P_b  = ' num2str(Pb_meas(i+j),4) ' W']);   
      PxxHeave = energy{i+j};
      f = freq{i+j};
      save(['WaveSpectrum_'  datestr(serialdate(i+j),'mmmm_dd_yyyy_HH-MM') '.mat'],'f','PxxHeave','serialdate',');
    
    end  
  end
  print(['WaveConditions_' datestr(serialdate(i),'mmmm_dd_yyyy_HH-MM')],'-dpdf');
  close all

end






xi = MeanSerialDate(1):(MeanSerialDate(end)-MeanSerialDate(1))/1000:MeanSerialDate(end);
yi_max = interp1(serialdate,Pb_meas,xi);
yi_meas = interp1(MeanSerialDate,MeanRectifiedPower,xi);
yi_Hs = interp1(serialdate,Hs,xi);
yi_Tp = interp1(serialdate,Tp,xi);




figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
set(gca,'FontSize',20);
plot(yi_Hs,yi_meas,'x');
title('June 26, 2012 through July 10, 2012');
xlabel('Significant Wave Height (m)');
ylabel('Rectified Power (W)');
print('RectifiedPowerVsWaveHeight','-dpdf');


figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
set(gca,'FontSize',20);
plot(yi_Tp,yi_meas,'x');
title('June 26, 2012 through July 10, 2012');
xlabel('Peak Period (Tp)');
ylabel('Rectified Power (W)');
print('RectifiedPowerVsPeakPeriod','-dpdf');



figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');

[AX,H1,H2] = plotyy(xi,yi_meas,xi,yi_Hs);
set(AX(1),'FontSize',20);
set(AX(2),'FontSize',20);
title('Electrical Power and Sig. Wave Height for June 26, 2012 through July 10, 2012');
legend('Rectified Power','Sig. Wave Height');

xlabel('Date');
datetick('x','keepticks');


set(AX(1),'YLim',[0 1000]);
set(AX(1),'YTick',0:200:1000);

set(AX(2),'YLim',[0 2]);
set(AX(2),'YTick',0:.2:2);

set(get(AX(1),'Ylabel'),'String','Rectified Power (W)'); 
set(get(AX(2),'Ylabel'),'String','Sig. Wave Height (m)');
set(get(AX(1),'Ylabel'),'FontSize',20);
set(get(AX(2),'Ylabel'),'FontSize',20);
set(AX(2),'XTick',[]);

print('ElectricalPowerVsWaveHeight','-dpdf');


figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
set(gca,'FontSize',20);
plot(xi,yi_meas,xi,yi_max);
title('Electrical Power for June 26, 2012 through July 10, 2012');
legend('Rectified Power','Max Absorbable Power');
xlabel('Date');
ylabel('Watts');
datetick('x','keepticks');
print('ElectricalPowerComparison','-dpdf');


figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
[AX,H1,H2] = plotyy(xi,100*yi_meas./yi_max,xi,yi_Hs);
set(H1,'LineWidth',2)
set(AX(1),'FontSize',20);
set(AX(2),'FontSize',20);
title('Performance for June 26, 2012 through July 10, 2012');

xlabel('Date');
datetick('x','keepticks');

set(AX(1),'YLim',[0 20]);
set(AX(1),'YTick',0:2:20);

set(AX(2),'YLim',[0 2]);
set(AX(2),'YTick',0:.2:2);

set(get(AX(1),'Ylabel'),'String','Percentage of Theoretical Maximum'); 
set(get(AX(2),'Ylabel'),'String','Sig. Wave Height (m)');
set(get(AX(1),'Ylabel'),'FontSize',20);
set(get(AX(2),'Ylabel'),'FontSize',20);
set(AX(2),'XTick',[]);

print('PercentageOfTheoreticalMax','-dpdf');



figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
set(gca,'FontSize',20);
plot(yi_Hs,100*yi_meas./yi_max,'x');
title('June 26, 2012 through July 10, 2012');
xlabel('Significant Wave Height (m)');
ylabel('Percentage of Theoretical Maximum');
print('PercentageOfMaxVsWaveHeight','-dpdf');


figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
set(gca,'FontSize',20);
plot(yi_Tp,100*yi_meas./yi_max,'x');
title('June 26, 2012 through July 10, 2012');
xlabel('Peak Period (s)');
ylabel('Percentage of Theoretical Maximum');
print('PercentageOfMaxVsPeakPeriod','-dpdf');


figure('Position',[19 45 1538 1093],'PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5],'PaperOrientation','landscape');
set(gca,'FontSize',20);
scatter3(yi_Hs(10:920),yi_Tp(10:920),100*yi_meas(10:920)./yi_max(10:920));
title('June 26, 2012 through July 10, 2012');
xlabel('Significant Wave Height (m)');
ylabel('Peak Period (s)');
zlabel('Percentage of Theoretical Maximum');

print('PercentageOfMaxVsPeakPeriodAndWaveHeight','-dpdf');