close all
clear all

rho = 1025;
g = 9.81;


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

path = 'CDIPData\';
wave_samplefreq = 1.28;

newdata = importdata([path 'xy_156_0e54701_noheaders.txt']);

%First sample: 20111108000200      

first_samp_datenum = datenum(2011,11,8,0,2,0)-8/24;

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

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


for nn =   18:1:(length(newdata)/4608)-1
    nn

idx_inwater = (nn-1)*4608+1;  %4608 samples per hour
idx_outofwater = (nn)*4608+1;
idx = idx_inwater:idx_outofwater;

timedate = newdata(idx,1);  
north = newdata(idx,2)/100;  %heave displacement in meters
west = newdata(idx,3)/100;
heave = newdata(idx,4)/100;
%clear newdata;
t = time(idx)-time(idx(1));
t_datenum = first_samp_datenum+time(idx)/(24*3600);

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



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

 Hs = 1:3;
 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(nn) = 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(nn) = 2*V*c0*sqrt(sum(Pow_meas))


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

%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

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


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(1),'mmmm-dd-yyyy HH:MM:SS') ' and '  datestr(t_datenum(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');
legend(['Measured Wave Spectrum (P = ' num2str(P_measured(nn),5) ' W/m    Pb = ' num2str(Pb_meas(nn),5) ' W)'],['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)']);

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

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

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,:),'linewidth',2);
set(gca,'FontSize',20);
xlabel('Frequency (Hz)');
ylabel('m^2/Hz');
legend(['Measured Wave Spectrum (Pb = ' num2str(Pb_meas(nn),4) ' W)'],['PM Spectrum, Hs = 1 m (Pb = ' num2str(Pb_Hs(1),4) ' W)'],['PM Spectrum, Hs = 2 m  Pb = ' num2str(Pb_Hs(2),4) ' W)'],['PM Spectrum, Hs = 3 m  (Pb = ' num2str(Pb_Hs(3),4) ' W)']);
title(['Wave conditions at ' datestr(t_datenum(1),'mmmm-dd-yyyy HH:MM:SS') ' (PDT)']);

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



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');

end



plot(mean_datenum(18:end),Pb_meas(18:end));


save('HourlyWaveData','mean_datenum','first_samp_datenum','time','ff','S_heave','heave_timeseries');