% readAdcpDSpec.m

% Hacked from plotCountourSpec2.m

% Read in reformatted RDI DSpec files - code snippet for Jeff Hanson.


% file names for processing ADCP - a list of specInfo Files
%[Dname, VSPname, InfoName] = textread('FileList', '%s %s %s');
[InfoName] = textread('FileList', '%s');

maxSf = 0;
d2r=pi/180;
minFreqDirSpec = 0.05;
maxFreq = 0.32;
minFreq = 0.04;
% ///////////////// MAIN LOOP ////////////////////////

for ii=1:length(InfoName)
% for ii=1:2

  % read specInfo file
[adcpID, projectDir, deployID, deployNum, deltad, ndirs, dfrdi, nfreqs, dir1, wavePlotDir, ccw, depthOffset, daysPerPlot ] = ...
   textread(char(InfoName(1,:)),'%s %s %s %d %d %d %f %d %d %d %d %f %d');
% ndirs = 360 / deltad;

dname=char(strrep(InfoName, 'SpecInfo', 'DSpec'));

dirs=[dir1 : deltad : dir1+(ndirs-1)*deltad ];

freqs=[dfrdi*0.5: dfrdi : (nfreqs - 0.5) * dfrdi]';
%i=find(freqs >= minFreq & freqs <= maxFreq);            % ??? add find for > min freq ???
%idir=find(freqs >= minFreqDirSpec & freqs <= maxFreq);            % ??? add find for > min freq ???

%dname=char(Dname(ii,:));
dpos=findstr(dname,'DSpec');
yrpos=findstr(dname,'.200')+1;
dotpos=findstr(dname,'.');
baseName = dname(dpos(end):dotpos(end)-1);
starttime =  dname(yrpos(end):dotpos(end-1)-1);      % one less dotpos since .UTC is appended to date

year=str2num(starttime(1:4));
mon=str2num(starttime(5:6));
day=str2num(starttime(7:8));
time=str2num(starttime(9:12));
hr=str2num(starttime(9:10));
minute=str2num(starttime(11:12));
sdate=datenum([year mon day hr minute 0]);

% ------------  RDI Dspec, Vspec, Spec, Pspec  ----------------------------
rdiRdir = load(dname);         % load rdi directional spectra (units are per deg) (freq,deg)
% rotate the directions
dirsR=[dir1 : deltad : dir1+(ndirs-1)*deltad ];

% The following arranges the array to start directions from wavePlotDir and go CW
if (wavePlotDir < min(dirsR))
   j=find(dirsR >= wavePlotDir + 360);
   dirs(1:length(j)) = dirsR(j) - 360;
   rdi360(:,1:length(j)) = rdiRdir(:,j);        % get f-d spec
   i=find(dirsR < wavePlotDir +360);
   dirs(length(j)+1:length(j)+length(i)) = dirsR(i);
   rdi360(:,length(j)+1:length(j)+length(i)) = rdiRdir(:,i);        % get f-d spec
else
   j=find(dirsR >= wavePlotDir);
   dirs(1:length(j)) = dirsR(j);
   rdi360(:,1:length(j)) = rdiRdir(:,j);        % get f-d spec
   i=find(dirsR < wavePlotDir);
   dirs(length(j)+1:length(j)+length(i)) = dirsR(i);
   rdi360(:,length(j)+1:length(j)+length(i)) = rdiRdir(:,i);        % get f-d spec
end
i=find(dirs >= wavePlotDir);
dirs(i) = dirs(i) - 360;

rdi360 = (rdi360 ./ 1000000);         % mm^2 to m^2
rdi360 = rdi360 ./ ndirs;

% get indicies for different frequency ranges, subranges for energy banding & "Z-tests"
i=find(freqs >= 0.04 & freqs <= 0.34);

Hmo=4*sqrt(dfrdi*sum(sum(rdi360(i,:))));

% ----------------------- Plot -----------------------------

if (1==0)  % change if you want plots

%  Frequency-Direction Spectra - single contour figure
figure
% contour
h = axes('Position',[0 0 1 1],'Visible','off');
axes('Position',[0.1 0.15 0.8 0.75])
i=find(freqs >= 0.04 & freqs <= maxFreq);

  contour(dirs, freqs(i), rdi360(i,:),20);
colorbar

set(gca,'View',[0 90])
%set(gca,'YLim',[0 0.7])
ylabel('Frequency (Hz)','FontSize', 12)
xlabel('Direction (TN)','FontSize', 12)

%hmo=sprintf('%4.2fm',Hv);
%txt=['F-D Spectra,       Hmo=' hmo];
txt=['ADCP F-D Spectra, ' char(adcpID) ':  ' datestr(sdate,0)];
title(txt,'FontSize', 14)
grid on

set(gcf,'CurrentAxes',h);
str = ['Energy (m^{2}/Hz/Deg)'];
text(0.95, 0.4, str, 'Rotation',90, 'FontSize', 12);

% delete(get(0,'CurrentFigure'));
end  % if 1==0

end   % main ii loop
