close all
clear all

rho = 1025;
g = 9.81;

c0 = rho*pi*g/4;
V = 2.2;

path = 'WaveRiderData\01-05-2012- PowerBuoy Deployment\';
wave_samplefreq = 1.28;

newdata = importdata([path 'DisplacementData.txt'],',',4);

wavedata = newdata.data(round(47*60/wave_samplefreq):end,2:4);  %Strip first 47 minutes to start data at 11:00 Pacifict Time.

first_samp_datenum = datenum(2012,1,5,19,0,0)-8/24;

time = 0:length(wavedata)-1;
time = time/wave_samplefreq;

dt = time(2)-time(1);

Pb_history = zeros(1,(length(wavedata)/4608));

for nn = 0:(length(wavedata)/4608)-1
    nn

if (nn > 0)
  idx = (nn-1)*4608+1:nn*4608;
  n_window = 64;
else
  idx = 1:length(wavedata);
  n_window = 256; 
end

heave = wavedata(idx,1)/100;  %heave displacement in meters
north = wavedata(idx,2)/100;  
west = wavedata(idx,3)/100;

%clear wavedata;
t = time(idx)-time(idx(1));
t_datenum = first_samp_datenum+time/(24*3600);

[f, PxxHeave] = tseries2spectrum(t,heave,n_window);

save([path 'WaveConditions_' datestr(t_datenum(idx(1)),'mmmm-dd-yyyy HH_MM') '-'  datestr(t_datenum(idx(end)),'mmmm-dd-yyyy-HH_MM') ' (PDT)']);

%Compute Pierson Moskowitz spectrum 
alpha = 8.1e-3;
beta = .74;
g = 9.81;
w = 2*pi*f;

 Hs = 1:3;
 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);
 end

%Compute energy flux for each spectrum (PM and measured)

Integrand = PxxHeave./f;
Integrand(1) = 0;

P_measured = 0.5*rho*g*trapz(f,Integrand);

df = f(2)-f(1);
%Compute upper bound for measured spectrum 
  %Pow_meas = 2*df*PxxHeave.*(f.^2);  % Compute upper bound of power extractable from each component
  %Pb_meas = 2*V*c0*sqrt(sum(Pow_meas))
  Pb_meas = BudalLimit(PxxHeave,f,V);


save([path 'WaveConditions_' datestr(t_datenum(idx(1)),'mmmm-dd-yyyy HH_MM') '-'  datestr(t_datenum(idx(end)),'mmmm-dd-yyyy-HH_MM') ' (PDT)'],'t','t_datenum','f','PxxHeave','heave','north','west');

for i = 1:length(Hs)
  Integrand = S_Hs(i,:)'./f;
  Integrand(1) = 0; 
  P_Hs(i) = 0.5*rho*g*trapz(f,Integrand);
end


if (nn > 0)
%Reconstruct wave profile for sanity check.
tt = 0:.1:t(end);
[zt] = spectrum2tseries(f,PxxHeave,tt);

%Compute upper bound for each PM spectrum 
for i = 1:length(Hs)  % Compute upper bound of power extractable from each component
Pow(i,:) = 2*df*S_Hs(i,:).*(f'.^2);
Pb_Hs(i) = 2*V*c0*sqrt(sum(Pow(i,:)));
end

figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
subplot(2,1,1)
plot(t,heave,tt,zt);
set(gca,'FontSize',14);
xlim([0 240]);
xlabel('time (seconds)');
ylabel('Elevation (m)');
legend('Measured Wave Elevations','Reconstructed Wave Elevations');
title(['Wave conditions between ' datestr(t_datenum(idx(1)),'mmmm-dd-yyyy HH:MM:SS') ' and '  datestr(t_datenum(idx(end)),'mmmm-dd-yyyy HH:MM:SS') ' (PDT)']);


subplot(2,1,2)
plot(f,PxxHeave,f,S_Hs(1,:),f,S_Hs(2,:),f,S_Hs(3,:));
set(gca,'FontSize',14);
xlabel('Frequency (Hz)');
ylabel('m^2/Hz');
foo = axis; foo(2) = .4; foo(4) = 20; axis(foo); clear foo;

%legend(['Measured Wave Spectrum (P = ' num2str(P_measured,5) ' W/m)'],['PM Spectrum, Hs = 1 m (P = ' num2str(P_Hs(1),5) ' W/m)'],['PM Spectrum, Hs = 2 m  P = ' num2str(P_Hs(2),5) ' W/m)'],['PM Spectrum, Hs = 3 m  (P = ' num2str(P_Hs(3),5) ' W/m)']);
legend(['Measured Wave Spectrum (Pb = ' num2str(Pb_meas,5) ' W)'],['PM Spectrum, Hs = 1 m (Pb = ' num2str(Pb_Hs(1),5) ' W)'],['PM Spectrum, Hs = 2 m  Pb = ' num2str(Pb_Hs(2),5) ' W)'],['PM Spectrum, Hs = 3 m  (Pb = ' num2str(Pb_Hs(3),5) ' W)']);

mean_datenum(nn) = (t_datenum(end)-t_datenum(1))/2;
S_heave{nn} = PxxHeave;
ff{nn} = f;
heave_timeseries{nn} = heave;


Pb_history(nn) = Pb_meas;
date_history(nn) = t_datenum(idx(1));
datestr_history{nn} = datestr(t_datenum(idx(1)),'mmmm-dd-yyyy HH:MM:SS');


else
figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
plot(f,PxxHeave,f,S_Hs(1,:),f,S_Hs(2,:),f,S_Hs(3,:));
set(gca,'FontSize',14);
xlabel('Frequency (Hz)');
ylabel('m^2/Hz');
foo = axis; foo(2) = .4; foo(4) = 15; axis(foo); clear foo;
legend(['Measured Wave Spectrum (P = ' num2str(P_measured,5) ' W/m)'],['PM Spectrum, Hs = 1 m (P = ' num2str(P_Hs(1),5) ' W/m)'],['PM Spectrum, Hs = 2 m  P = ' num2str(P_Hs(2),5) ' W/m)'],['PM Spectrum, Hs = 3 m  (P = ' num2str(P_Hs(3),5) ' W/m)']);
title(['Wave conditions between ' datestr(t_datenum(idx(1)),'mmmm-dd-yyyy HH:MM:SS') ' and '  datestr(t_datenum(idx(end)),'mmmm-dd-yyyy HH:MM:SS') ' (PDT)']);  


end


print('-dpdf',[path 'WaveConditions_' datestr(t_datenum(idx(1)),'mmmm-dd-yyyy HH_MM') '-'  datestr(t_datenum(idx(end)),'mmmm-dd-yyyy-HH_MM') ' (PDT)']);


close all

end