ar all
close all
clear all

rho = 1025;
g = 9.81;

path = 'CDIPData\';
wave_samplefreq = 1.28;

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

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

%First sample: 20111108000200      

first_samp_datenum = datenum(2011,11,11,11,0,0);

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

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

idx = 1:length(newdata);

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 = 0.5*rho*g*trapz(f,Integrand)


for i = 1:length(Hs)
  Integrand = S_Hs(i,:)'./f;
  Integrand(1) = 0; 
  P_Hs(i) = 0.5*rho*g*trapz(f,Integrand);
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,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)']);
foo = axis;
foo(2) = .35;
axis(foo);

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 = (t_datenum(end)-t_datenum(1))/2;
S_heave = PxxHeave;
ff = f;
heave_timeseries = heave;

df = f(2)-f(1);
PV_measured = 0;
for j = 1:length(f)
   PV_measured = PV_measured + f(j)*sqrt(2*PxxHeave(j)*df);
end
PV_measured = PV_measured*0.5*pi*rho*g;


for i = 1:length(Hs)
    PV_Hs(i) = 0;
    for j = 1:length(f)
      i 
      j
      PV_Hs(i) = PV_Hs(i) + f(j)*sqrt(2*S_Hs(i,j)*df)
    end
end
PV_Hs = PV_Hs*0.5*pi*rho*g;

figure('Position',[196  79  1004  598],'PaperOrientation', 'landscape','PaperPositionMode','manual','PaperType','usletter','PaperPosition',[0 0 11 8.5]); 
subplot(2,1,1)
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,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)']);
foo = axis;
foo(2) = .35;
axis(foo);
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(PV_measured,5) ' W/m^3)'],['PM Spectrum, Hs = 1 m (P = ' num2str(PV_Hs(1),5) ' W/m^3)'],['PM Spectrum, Hs = 2 m  P = ' num2str(PV_Hs(2),5) ' W/m^3)'],['PM Spectrum, Hs = 3 m  (P = ' num2str(PV_Hs(3),5) ' W/m^3)']);
foo = axis;
foo(2) = .35;
axis(foo);










%plot(mean_datenum(18:end),P_measured(18:end));


save([path 'xy_Nov10_0900PDT'],'mean_datenum','first_samp_datenum','time','ff','S_heave','heave_timeseries');