%wavedir.m
% compute wave direction based on velocity and pressure
% Meg Palmsten, USGS St. Pete, 3/00

if exist ('noverlap') == 0, noverlap = 256; end % # to overlap on fft


inpathfile = [inpath, datfil];
	
mnu = mean(urr);
mnv = mean(vrr);

% compute cross-spectra between velocities and pressure

[Spu, freq] = csd(px, urr,  blk, sampfreq, 512, noverlap);
[Spv, freq] = csd(px, vrr,  blk, sampfreq, 512, noverlap);

% compute wave direction

phase = atan2(Spv, Spu);
dirs = (phase).*(180/pi);

% compute wave height from velocities

[Guu, Gvv, hgt_uv] = waveht_uv(urr, vrr, zdep, shgtsens, sampfreq);

% compute ratio of hsig pressure to hsig velocity

rhsig = hgt/hgt_uv;

% angular frequency

lastfreq = max(find(freq<=0.20));
omega = 2.*pi.*freq;
fomega = (fix(omega*10000))/10000;

% determine direction at rtp
if rtp == 0
	frfp = 0; %to avoid divide by zero error
end
if rtp>0
	rfp = 2.*pi./rtp;
	frfp = (fix(rfp*10000))/10000;
end
ind = find(fomega == frfp );
dirtp = dirs(ind);

% plot

set(figh,'nextplot','replacechildren','visible','off','selectionhighlight','off')
axh = newplot;
set(axh,'visible','off')

axes('position', [.1, .57, .8, .38]);
semilogy(freq(2:lastfreq), Snn(2:lastfreq))
ylabel('m^2/s');
title('Depth Corrected Pressure Spectra')
ylim([10^-2, 10^2])

axes('position', [.1, .1, .8, .38]);
plot(freq(2:lastfreq), dirs(2:lastfreq))
ylim([-200, 200]);
title('Depth Corrected Directional Spectra');
xlabel('Frequency');
hgts = num2str(hgt);
dirtps = num2str(dirtp);
rhsigs = num2str(rhsig);
string = ['Hsig(p) = ', hgts,'   Dir = ', dirtps];
string2 = ['Hsig(p) : Hsig(u,v) = ', rhsigs];
text(.12, -250, string);
text(.12, -270, string2);

bstr = num2str(bnum);
name = [hdfil(1:2),hdfil(7:8),' ', bstr];

suptitle2(name) 

psfil = [hdfil(1:2), hdfil(7:8),'_', bstr, 'dir.ps'];
disp(['Creating: ', psfil])
eval(['print -dpsc ', imagepath, psfil]);

% collect array of data for file

alldat(1,bnum) = bnum;
alldat(2,bnum) = year;
alldat(3,bnum) = month;
alldat(4,bnum) = day;
alldat(5,bnum) = hour;
alldat(6,bnum) = minute;
alldat(7,bnum) = second;
alldat(8,bnum) = zdep;
alldat(9,bnum) = mnu;
alldat(10,bnum) = mnv;
alldat(11,bnum) = hgt;
alldat(12,bnum) = rhsig;
alldat(13,bnum) = rtp;
alldat(14,bnum) = dirtp;
alldat(15,bnum) = length(rep); %from badness2uv.m
alldat(16,bnum) = rmsdelu;	%from badness2uv.m
alldat(17,bnum) = rmsdelv;	%from badness2uv.m
alldat(18,bnum) = rmsdelup; 	%from badness2uv.m
alldat(19,bnum) = perrmsdelu;	%from badness2uv.m
alldat(20,bnum) = perrmsdelv;	%from badness2uv.m
alldat(21,bnum) = perrmsdelup; 	%from badness2uv.m
%end