%plotdirspec
%program to read AXYS data from ECS, 2JUN01, Peter H. Dahl on board Melville
%modified by PHD 6MAY03 to replace web-site program
clear

%  get into the DIRSPEC directory for reading data              
cd DIRSPEC  
pth=['200106032230.DIRSPEC'];  %3 JUNE 2001 2230 UTC good demonstration  
fid = fopen(pth,'r');

%Now the program will begin scaning the input file string for data
while 1
   header = fgetl(fid);
   if findstr(header, ' ETC ') > 0, break,end
end
%The header is now read through and data scanning can begin.
%Need to use fgetl to start after the header

for i = 1:1860			%As specified for data
   dirdat = fgetl(fid);
   data(i,:) = str2num(dirdat);
   %data(:,i) = sscanf(dirdat, '%f');		%alternative into columns
end

lastline = fgetl(fid);						%reads last three data points	
data = reshape(data',1,1860*8);				%changes from 1860,8 __ 1,1860*8
addpad = str2num(lastline);					%pads on last three data
data = [data, addpad];
data = reshape(data,121,123);				%creates the right size matrix
fclose('all');

cd ..  %return to original directory of seabreeze program

%Given spectral and directional ranges
f = .03:0.005:(0.005*122)+.03;
theta = 0:3:360;

clf
xx=find(f >= .09); %remove waves of periods > 11 s

test=data(:,xx);
sc2=10*log10(max(test(:)));
sc=[sc2-30 sc2];  %plot 20 dB of dynamic range
subplot(2,1,1)
imagesc((f(xx)),theta,10*log10(data(:,xx)),sc)
axis([.1 .6 0 360])
set(gca,'Xtick',[.1 .2 .3 .4 .5 .6])


h=xlabel('FREQUENCY (Hz)');set(h,'fontsize',14,'color','k');
h=ylabel('WAVE DIRECTION FROM (deg)');set(h,'fontsize',14,'color','k');
set(gca,'fontsize',14);


subplot(2,1,2)

SD=sum(data)*3;
h=semilogy(f,SD,'b');set(h,'linewidth',2)
hold

h=xlabel('FREQUENCY (Hz)');set(h,'fontsize',14,'color','k');
h=ylabel('SPECTRAL DENSITY m^{2}/Hz');set(h,'fontsize',14,'color','k');
set(gca,'fontsize',14);
set(gca,'Xtick',[.1 .2 .3 .4 .5 .6 ])
set(gca,'Ytick',[.01 .02 .04 .1 .2 .4 .8 1.6])
ymax=1.1*max(SD);
axis([.1 .6 0 ymax])


x=find(f >= .09);
H3=4*sqrt(trapz(f(x),SD(x)));
%put this in the title
subplot(2,1,1)
month=str2num(pth(6));
if month == 5
title([pth(7:8),'MAY ',pth(9:12),' UTC, H_{1/3} =',num2str(H3,2),' m' ])
else
title([pth(7:8),'JUNE ',pth(9:12),' UTC, H_{1/3} =',num2str(H3,2),' m' ])
end