close all
clear all

path = 'CDIP WaveData\';
wave_samplefreq = 1.28;

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



timedate = newdata(:,1);  
north = newdata(:,2)/100;  %heave displacement in meters
west = newdata(:,3)/100;
heave = newdata(:,4)/100;
clear newdata;
t = 0:length(heave)-1;
t = t/wave_samplefreq;

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

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

TotalEnergy = sum(PxxHeave*(f(2)-f(1)));

%Compute Pierson Moskowitz spectrum with same total energy.
alpha = 8.1e-3;
beta = .74;
g = 9.81;
w = 2*pi*f;
%for Hs = 1:-.1:.1;
%  w0 = sqrt(0.21*g/Hs);
%  S = (alpha*g^2)*exp(-beta*(w0./w).^4)./(w.^5);
%  sum(S(2:end))
%  if(sum(S(2:end)) < TotalEnergy)
%      PM_Hs = Hs;
%      break;
%  end
%end

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


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 19:30 and 22:15 (GMT) on January 25th, 2011');


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','PM Spectrum, Hs = 1 m','PM Spectrum, Hs = 2 m','PM Spectrum, Hs = 3 m');

print('-dpdf',[path 'Feb17-2011_WaveConditions']);
print('-dpng',[path 'Feb17-2011_WaveConditions']);
print('-dill',[path 'Feb17-2011_WaveConditions']);


figure
plot(t,heave,tt,zt);