function out = readNdbcSpec(buoy,dates,varargin)

%READNDBCSPEC  Reads NDBC spectral data (1D & 2D)
%   OUT = readNdbc2DSpec(BUOY,DATES) reads NDBC spectral data from 
%   BUOY files from (or between and including) DATE(S). OUT is a structure
%   with dates found, frequencies, directions, energy (variance) density, 
%   and calculated Hs.
%
%   OUT = readNdbc2DSpec(BUOY,DATES,METHOD) selects the METHOD for 
%   calculating the directional distribution function:
%        1 - LHM (Longuet-Higgins et al., 1963) (DEFAULT)
%        2 - LHM with weighting to aviod negative energies
%        3 - MLM Maximum Likelihood Method (Oltman-Shay and Guza, 1984)
%        4 - MEM Maximum Entropy Method (Lygre and Krogstad, 1986)
%
%   OUT = readNdbc2DSpec(BUOY,DATES,METHOD,DIRRES) allows for desired 
%   direction resolution.
%
%   Input:
%     BUOY = num or 'text'; (ex. buoy=42001 | buoy='../ndbc/42001')
%     DATE(S) = yyyymmddHH | [yyyymmddHH yyyymmddHH];
%     METHOD = 1 | 2 | 3 | 4
%     DIRRES = num
%   Output:
%     OUT = struct('t',dates_found,'f',frequencies,'d',directions,...
%           'S',2D_variance_density,'E',1D_variance_density,'Hs',Hs);
%        size(OUT.E) = [length(OUT.f) length(OUT.t)]
%        size(OUT.S) = [length(OUT.f) length(OUT.d) length(OUT.t)]
%        2D variance density has units of m^2/Hz/degr
%
%   !!!ASSUMPTIONS!!!:
%     1 - Beyond the input BUOY number, file name(s) MUST 
%        contain the NDBC designations as follows:
%           BUOYw* = Spectral wave density data
%           BUOYd* = Spectral wave (alpha1) direction data
%           BUOYi* = Spectral wave (alpha2) direction data
%           BUOYj* = Spectral wave (r1) direction data
%           BUOYk* = Spectral wave (r2) direction data
%        The '*' indicates that any text may occur after the NDBC
%        designations - examples: 42001w.dat, 42001w2004.blah, 
%        42001wmydata.txt). For any of these BUOY=42001 would be fed to
%        readNdbcSpec.
%     2 - Files contain the header line.
%     3 - Date form is YYYY MM DD HH. I've come across YY MM DD HH and 
%        YYYY MM DD HH mm in the 'Standard Meteorological' data reports but
%        not in the Spectral data.
%
%   out = readNdbc2DSpec(buoy,dates,{method,dirres});
%
% Dave Thompson (dthompson@usgs.gov)

%% Check inputs
if nargin==4
   opt = cell2mat(varargin(1));
   dirres = cell2mat(varargin(2));
elseif nargin==3
   opt = cell2mat(varargin(1));
   dirres = 1;
elseif nargin==2
   opt = 1;
   dirres = 1;
end
      
%% Check for data files.
% Get all the file names that make up the spectral data. The letter
% designations should be as they come from NDBC (i.e. w, d, i, j, k).
tok = ['w';'d';'i';'j';'k'];
fname = num2str(buoy);
direc = char(regexp(fname,'.+\/|.+\\','match'));
fnames = [];
for ii=1:length(tok)
   tmp = strtrim(evalc(['dir ',fname,tok(ii),'*']));
   if isempty(findstr(tmp,'not found')) 
      disp([' found ',direc,tmp])
      fnames = strvcat(fnames,[direc,tmp]);
   end
end
if isempty(fnames)
   disp([' ',fname,' file(s) not found!!'])
   out = [];
   return
end
%% Read in the data.
% E = spectral wave density data  -> E = (f,dates)
% D = spectral wave direction data -> D = (f,(alpha1,alpha2,r1,r2),dates)
% Edates = dates found in wave density data
% Ddates = dates found in wave direction data (BUOYd* file only)
% Efreqs = freqs found in wave density data
% Dfreqs = freqs found in wave direction data (BUOYd* file only)
% I made this thing store and check for differences in dates and 
% frequencies between the density and one of the direction files because 
% I've come across a cases where they didn't match up!  Ideally this is 
% overkill. The code is NOT checking for differences between the 4
% direction data files since I've never encountered that.
for mm=1:size(fnames,1) % Main loop is for each data file -> fnames counter
   jj = 1;              % dates counter
   fid = fopen(fnames(mm,:),'r');

   while 1  % Read file line-by-line.
      line = fgetl(fid);

      % Break if we've reached the end of the file.
      if ~ischar(line)
         break
      end

      % Get the frequencies listed in the header.
      % mm==2 line assumes that the direction files (d, i, j, k) all
      % have the same frequencies.
      if regexp(line,'^\D')
         dnum = find(isletter(line)==1,1,'last');
         if mm==1 % working on density data file
            Efreqs = sscanf(line(dnum+1:end),'%f');
         else   % working on the direction data files
            Dfreqs = sscanf(line(dnum+1:end),'%f');
         end
         % If you come across different forms of the date in the data make
         % the adjustment here.
         %if dnum==11
         %   dform = 'yymmddHH';
         %elseif dnum==13
            dform = 'yyyymmddHH';
         %elseif dnum==16
         %   dform = 'yyyymmddHHMM';
         %end
         continue
      end

      % Find the data of interest.
      ddate = str2num(line(isspace(line(1:dnum))==0));
      while ddate>=dates(1) && ddate<=dates(end)
         % mm==2 line assumes that the direction files (d, i, j, k) all
         % have the same dates.
         if mm==1 % working on density data file
            Edates(jj) = datenum(line(1:dnum),dform);
         else   % working on the direction data files
            Ddates(jj) = datenum(line(1:dnum),dform);
         end
         tmp = sscanf(line(dnum+1:end),'%f');
         tmp(tmp==999) = nan;
         if mm==1 % working on density data file
            E(:,jj) = tmp;
         else     % working on the direction data files
            D(:,mm-1,jj) = tmp;
         end
         jj = jj + 1;   % advance dates counter
         line = fgetl(fid);
         % Break if we've reached the end of the file.
         if ~ischar(line)
            disp([' Reached EOF before finding last date requested in ',...
               fnames(mm,:),' !!'])
            break
         end
         ddate = str2num(line(isspace(line(1:dnum))==0));
      end
   end

   if ~exist('E','var')
      disp(' Requested data not found in density data!')
      out = [];
      return
   end
   if mm==1
      if setdiff(diff(str2num(datestr(Edates(:,1),'HH'))),[1 -23])
         disp([' Did not find all dates requested in ',fnames(mm,:),' !!'])
      end
   else
      if setdiff(diff(str2num(datestr(Ddates(:,1),'HH'))),[1 -23])
         disp([' Did not find all dates requested in ',fnames(mm,:),' !!'])
      end
   end
   fclose(fid);
end

%% Calculate things.
if exist('D','var') % working with 2D spectra
   % Check for differences in dates and frequencies between the energy and 
   % direction data, tell me about it, and fix it.
   if sum(Edates-Ddates)~=0
      disp(' Energy and Direction data dates don''t match!!')
      [dates,idd,ide] = intersect(Ddates,Edates);
      D = D(:,:,idd);
      E = E(:,ide);
   else
      dates = Edates;
   end

   if length(Dfreqs)-length(Efreqs)~=0
      disp(' Energy and Direction data have different number of frequencies!!')
      [freqs,idd,ide] = intersect(Dfreqs,Efreqs);
      D = D(idd,:,:);
      E = E(ide,:);
   elseif sum(Efreqs-Dfreqs)~=0
      disp(' Energy and Direction frequencies don''t match!!')
      [freqs,idd,ide] = intersect(Dfreqs,Efreqs);
      D = D(idd,:,:);
      E = E(ide,:);
   else
      freqs = Efreqs;
   end
   % Calculate the directional distribution function.
   dirs = 0:dirres:359;
   if opt==1   % Longuet-Higgins et al. (1963) method.
      for ii=1:length(dates)      % for each date
         for jj=1:length(freqs)   % for each frequency
            for kk=1:length(dirs) % for each direction
               DD(jj,kk,ii) = dirres*(1/180)*(1/2 + ...
                  D(jj,3,ii)/100*cosd(dirs(kk)-D(jj,1,ii)) + ...
                  D(jj,4,ii)/100*cosd(2*(dirs(kk)-D(jj,2,ii))));
            end
         end
      end
   elseif  opt==2 % LHM with weighting to aviod negative energies.
      for ii=1:length(dates)      % for each date
         for jj=1:length(freqs)   % for each frequency
            for kk=1:length(dirs) % for each direction
               DD(jj,kk,ii) = dirres*(1/180)*(1/2 + ...
                  (2/3)*(D(jj,3,ii)/100)*cosd(dirs(kk)-D(jj,1,ii)) + ...
                  (1/6)*(D(jj,4,ii)/100)*cosd(2*(dirs(kk)-D(jj,2,ii))));
            end
         end
      end
   elseif opt==3 % MLM
      for ii=1:length(dates)      % for each date
         for jj=1:length(freqs)   % for each frequency
            DD(jj,:,ii) = (1/dirres)*mlm(D(jj,3,ii)/100,...
               D(jj,4,ii)/100,D(jj,1,ii),D(jj,2,ii),dirs);
         end
      end
   elseif opt==4  % MEM
      for ii=1:length(dates)      % for each date
         for jj=1:length(freqs)   % for each frequency
            DD(jj,:,ii) = (1/dirres)*mem(D(jj,3,ii)/100,...
               D(jj,4,ii)/100,D(jj,1,ii),D(jj,2,ii),dirs);
         end
      end
   end
   % Finally calculate the 2D spectra.
   for ii=1:length(dates)     % for each date
      for jj=1:length(dirs)   % for each direction
         S(:,jj,ii) = E(:,ii).*DD(:,jj,ii);
      end
   end

   % Calculate Hs - S needs to be interpolated onto an equal interval freq
   % space for this.
   fi = [freqs(1):min(diff(freqs)):freqs(end)]';
   for ii=1:size(S,3)
      Hs(ii) = 4*sqrt(trapz(trapz(interp2(dirs,freqs,S(:,:,ii),dirs,fi))*...
         (fi(2)-fi(1))*dirres));
   end
else % Working with 1D spectra.
   dates = Edates;
   freqs = Efreqs;
   dirs = []; S = [];
   % Calculate Hs - E needs to be interpolated onto an equal interval freq
   % space for this.
   fi = [freqs(1):min(diff(freqs)):freqs(end)]';
   for ii=1:size(E,2)
      Hs(ii) = 4*sqrt(trapz(interp1(freqs,E(:,ii),fi))*(fi(2)-fi(1)));
   end
end
   
%% Output.
out = struct('t',dates','f',freqs,'d',dirs','S',S,'E',E,'Hs',Hs');

return
%% Some extra stuff...
% Use this to plot the spectra - this chunk 'should' run properly on OUT.
% Works for both 1D and 2D spectra... plots both for 2D data.
% For a series, all plot limits are set to the max value of the series.
figure
if ~isempty(out.S)
   subplot(1,2,1)
end
he = semilogy(out.f,out.E(:,1),'r','LineWidth',2);
% semilogy(out.f,sum(out.S(:,:,1),2)) % should fall over top of out.E(:,ii)
grid on
axis([0 max(out.f) 0 max(max(out.E))])
xlabel('frequency (Hz)'); ylabel('E (m^2/Hz)')
set(gca,'XMinorTick','on')
het = title(datestr(out.t(1),'mm/dd/yy HHMM'),'FontSize',12,...
   'FontWeight','bold');
if ~isempty(out.S)
   subplot(1,2,2)
   % initialize polar plot
   % fix TTickLabels to get rid of mmpolars +/- numbering scheme
   for ii=1:12
      tmp = -30*ii + 300;
      tmp(tmp<0) = tmp + 360;
      ttick(ii) = cellstr([num2str(tmp),'\circ']);
   end
   h = mmpolar([0 2*pi], [min(out.f) max(out.f)],...
      'RTickUnits','Hz','FontWeight','bold','TTickLabel',ttick);
   delete(h)
   hold on
   
   % build spectral mesh and convert for polar type plot
   [Ndg,Nfg] = meshgrid([0;out.d],out.f);
   [Npx,Npy] = pol2cart((90-Ndg)*pi/180,Nfg);
   
   maxS = max(max(max(out.S)));
   c = maxS;
   for ii=1:11
      c = [c(1)/2; c];
   end
   cmap = jet(length(c)-1); 
   
   tmp = log10(out.S(:,:,1)); 
   tmp = [tmp(:,1) tmp];
   hs = pcolor(Npx*2,Npy*2,tmp); shading flat
   colormap(cmap)
   caxis(log10([c(1) c(end)]))
   set(gca,'DefaultTextUnits','normalized')
   hst = text(.5,1.1,datestr(out.t(1),'mm/dd/yy HHMM'),'FontSize',12,...
      'FontWeight','bold','HorizontalAlignment','center');
   
   ca = colorbar('Location','SouthOutside');
   set(ca,'XTick',log10(c),'XTickLabel',num2str(c,'%.5f'));
   axes(ca)
   rotateticklabel(ca);
end

if length(out.t)>1
   warning off
   disp(' [Enter] to continue')
   pause
   for ii=2:length(out.t)
      set(he,'YData',out.E(:,ii))
      set(het,'String',datestr(out.t(ii),'mm/dd/yy HHMM'))
      
      if ~isempty(out.S)
         tmp = log10(out.S(:,:,ii)); 
         tmp = [tmp(:,1) tmp];
         set(hs,'CData',tmp)
         set(hst,'String',datestr(out.t(ii),'mm/dd/yy HHMM'))
      end
      
      disp(' [Enter] to continue')
      pause
   end
   warning on
end