close all
clear all


foo = ls;
j = 0;
for i = 1:size(foo,1)
  if findstr('.raw',foo(i,:))
    j=j+1;
    filenames{j} = foo(i,:);
  end
end
clear foo j;

displacements = [];
time = [];
for i = 1:length(filenames)
  foo = load(filenames{i});

  YR = str2double(filenames{i}(1:4));
  MO = str2double(filenames{i}(6:7));
  DA = str2double(filenames{i}(9:10));
  HR = str2double(filenames{i}(12:13));
  MN = str2double(filenames{i}(15:16));
  SE = 0;
  t0 = datenum(YR,MO,DA,HR,MN,SE);
  dt = (0.5/24)/length(foo);
  t = t0:dt:t0+.5/24-dt;
  
  time = [time t];
  displacements = [displacements; foo];
end
clear i foo YR MO DA HR MN SE T0 dt t;

%At this point all the files in the directory have been put sequentially
%into the displacements array and a cooresponding time array of matlab
%datenum's has been created. Gaps in the files will show up as gaps in the
%data.


Hs = 4*std(displacements(:,2));
TD_Pow = sum(abs(displacements(:,2).^2))/length(time)

StartTime = datenum(2010,8,12,20,27,0);
StopTime = datenum(2010,8,12,21,04,0);
idx2 = find(time<=StopTime & time >= StartTime);
Hs2 = 4*std(displacements(idx2,2));
leg2 = ['With Skirt, before tow -> Hs = ' num2str(Hs2/100,2) ' m'];
TD_Pow2 = sum(abs(displacements(idx2,2).^2))/length(idx2)

StartTime = datenum(2010,8,12,21,10,0);
StopTime = datenum(2010,8,12,21,42,0);
idx3 = find(time<=StopTime & time >= StartTime);
Hs3 = 4*std(displacements(idx3,2));
leg3 = ['With Skirt, after tow -> Hs = ' num2str(Hs3/100,2) ' m'];
TD_Pow3 = sum(abs(displacements(idx3,2).^2))/length(idx3)

StartTime = datenum(2010,8,12,22,38,0);
StopTime = datenum(2010,8,12,23,05,0);
idx4 = find(time<=StopTime & time >= StartTime);
Hs4 = 4*std(displacements(idx4,2));
leg4 = ['Without Skirt -> Hs = ' num2str(Hs4/100,2) ' m'];
TD_Pow4 = sum(abs(displacements(idx4,2).^2))/length(idx4)


figure('PaperSize',[11 17],'PaperType','tabloid','Position',[503 73 838 1057],'PaperPosition',[.5 .5 10,16]);

subplot(4,1,1)
plot(time-7/24,displacements(:,2))
title('WaveRider Heave Displacements - August 12, 2010');
datetick('x');
xlabel('Time of Day, Pacific Time');
ylabel('Heave Displacement (cm)');
leg = ['Entire Record -> Hs = ' num2str(Hs/100,2) ' m'];
legend(leg);

subplot(4,1,2)
plot(time(idx2)-7/24,displacements(idx2,2))
datetick('x');
xlabel('Time of Day, Pacific Time');
ylabel('Heave Displacement (cm)');
legend(leg2);

subplot(4,1,3)
plot(time(idx3)-7/24,displacements(idx3,2))
datetick('x');
xlabel('Time of Day, Pacific Time');
ylabel('Heave Displacement (cm)');
legend(leg3);

subplot(4,1,4)
plot(time(idx4)-7/24,displacements(idx4,2))
datetick('x');
xlabel('Time of Day, Pacific Time');
ylabel('Heave Displacement (cm)');
legend(leg4);

print -dpdf ElevationHistory-Aug12-2010


window = length(displacements(:,2))/512;
overlap = [];  %Defaults to 50%
Fs = 1.28;  %Sampling Frequency
[Pxx_all,f] = pwelch(displacements(:,2),window,overlap,[],Fs);
[Pxx_2,f2] = pwelch(displacements(idx2,2),window,overlap,[],Fs);
[Pxx_3,f3] = pwelch(displacements(idx3,2),window,overlap,[],Fs);
[Pxx_4,f4] = pwelch(displacements(idx4,2),window,overlap,[],Fs);

FD_Pow = sum(Fs*Pxx_all)/(2*(length(f)-1))
FD_Pow2 = sum(Fs*Pxx_2)/(2*(length(f2)-1))
FD_Pow3 = sum(Fs*Pxx_3)/(2*(length(f3)-1))
FD_Pow4 = sum(Fs*Pxx_4)/(2*(length(f4)-1))
    


figure('Position',[225 377 1235 725],'PaperOrientation','landscape','PaperPosition',[.5 .5 10,7.5]);
plot(f,Fs*Pxx_all/10000,f2,Fs*Pxx_2/10000,f3,Fs*Pxx_3/10000,f4,Fs*Pxx_4/10000);  %Plotted as m^2/Hz per Hz
legend(leg,leg2,leg3,leg4);
xlabel('Frequency f (Hz)');
ylabel('S(f) (m^2/Hz)');
title('PSD comparisons of WaveRider Heave Motions - August 12, 2010');

print -dpdf PSD_Comparisons-Aug12-2010

