function varargout = swan_io_spectrum(varargin)
%SWAN_IO_SPECTRUM         read SWAN 1D or 2D spectrum file.
%
%    DAT = swan_io_spectrum(fname)
%    DAT = swan_io_spectrum  % launches file load GUI
%    DAT = swan_io_spectrum(INP.spec(i))
% 
% where fname is the table file name (e.g. *.s1d, *.s2d)
% where INP.spec is returned by INP = swan_io_input('INPUT')
%
% Reads an ASCII SWAN spectrum file into a struct DAT which looks like:
%
%                     filename: 'r01n1t01.sp1'
%                         path: ''
%                         name: 'r01n1t01'
%                          ext: '.sp1'
%                         time: 0
%          number_of_locations: 1
%                            x: 400000
%                            y: 6440000
%               frequency_type: 'absolute'
%        number_of_frequencies: 25
%                    frequency: [1x25 double]
%        dimension_of_spectrum: 0/1/2
%         number_of_quantities: 3
%    quantity_exception_values: [-99 -999 -9]
%               quantity_names: {'EnDens'  'NDIR'  'DSPRDEGR'}
%               quantity_units: {'J/m2/Hz'  'degr'  'degr'}
%                         data: [3x25 double]
%                       EnDens: [1x25 double]
%                         NDIR: [1x25 double]
%                     DSPRDEGR: [1x25 double]
%             data_description: [3x45 char]
%                      read_at: 'dd-mmm-yyy hh:mm:ss'
%                     iostatus: 0
%
% NOTE: SWAN_IO_SPECTRUM reads also PAR, S1D and S2D files generated by 
%       keyword TEST, as well as the HOTFile (and the new MUDFile).
%
% See also: SWAN_IO_INPUT, SWAN_IO_TABLE, SWAN_IO_BOT, SWAN_IO_GRD

%% Copyright
%   --------------------------------------------------------------------
%   Copyright (C) 2006-2009 Deltares
%       Gerben de Boer
%
%       gerben.deboer@deltares.nl	
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

%% Version
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$

% 2009 Jun 05: added option to pass data field from swan_io_input, that contains fullfilename
% 2009 Jun 04: use new matlab code-cells syntax to divide code into 'chapters'
% 2008 May 07: made it work for SPEC2D with NODATA, 
%              Improved debug listing with mod, 
%              Replaced i's with iloc, idir etc. 
%              Pre-allcoated SPEC2D for speed.
% 2008 May 13: read n_locations <, mxc, myc> after keyword LOCATION<S> (introduced in mud special for MUDF and for HOTF)) 
%              and optionally plot reshaped matrix
%              also read long quantity names
%              use fscanf to read co-ordinates order faster
%              pre-allocate 1D data blocks
% 2008 Oct 17: made 2D output [nloc x  nfreq x ndir]
%              give 2D output array name of quantity
% 2009 Apr 23: added to check for existance of file
% 2013 Jan 16: Compatibility with 1D spectral outputs of SWAN Ver. 4.8** or previous (see line 536)

%TO DO:

%% try to harmonize output with dummy dimensions of length 1 ??
%             DAT.data_description = ['1st dimension = number_of_locations          ';
%                                     '2st dimension = number_of_frequencies        ';
%                                     '3rd dimension = number_of_directions         ';
%                                     '4th dimension = number_of_iterations         '];
%                                        number_of_locations mx
%                                        number_of_locations my

%% Defaults

OPT.debug = [0 0 0]; % 1st is tree, 2nd is all lines, 3rd is pcolor of reshaped (x,y) matrix
OPT.mod   = 1000; % for OPT.debug(2)
OPT.fast  = 1; % use fscanf, this does not allow comments in between the data lines 
               % which SWAN does not add (but some external spectrum write tools do)

%% No file name specified if even number of arguments
%  i.e. 2 or 4 input parameters

   if mod(nargin,2)     == 0 
     [fname, pathname, filterindex] = uigetfile( ...
        {'*.sp*;*.s*d', 'SWAN spectrum files (*.sp*;*.s*d)'; ...
         '*.*'   ,'All Files (*.*)'}, ...
         'SWAN 1D spectrum file');
      
      if ~ischar(fname) % uigetfile cancelled
         DAT.filename   = [];
         iostat         = 0;
      else
         DAT.fullfilename = [pathname, fname];
         iostat           = 1;
      end

   %% No file name specified if odd number of arguments
   elseif mod(nargin,2) == 1 % i.e. 3 or 5 input parameters
      
      if isstruct(varargin{1})
      
      DAT = varargin{1}; % only field fullfilename required
      
      elseif ischar(varargin{1})

      DAT.fullfilename  = varargin{1};
      
      end
      
      iostat            = 1;
   end
   
   [DAT.file.path DAT.file.name DAT.file.ext] = fileparts(DAT.fullfilename);   

%% Open file

if iostat==1 %  0 when uigetfile was cancelled
             % -1 when uigetfile failed

%% check for file existence (1)

   tmp = dir(DAT.fullfilename);
   
   if length(tmp)==0
      
      if nargout==1
         error(['Error finding file: ',DAT.fullfilename]);
      else
         iostat = -1;
      end   
      
   elseif length(tmp)>0
   
      DAT.file.date     = tmp.date;
      DAT.file.bytes    = tmp.bytes;
      
%% check for file opening (2)

      fid       = fopen  (DAT.fullfilename,'r');
      
      if fid < 0
         
         if nargout==1
            error(['Error opening file: ',DAT.fullfilename]);
         else
            iostat = -2;
         end      
         
      elseif fid > 2
         
%% check for file reading (3)

         %try
         
%% Read 1st line

            rec = fgetl_no_comment_line(fid,'$');
            if ~strcmp(strtok(upper(rec)),'SWAN')
                fclose(fid);
                error('no SWAN on first line')
            end
   
            rec = fgetl_no_comment_line(fid,'$');
   
%% Read TIME encoding (optional)

            if strcmp(strtok(upper(rec)),'TIME')
   
               rec = fgetl_no_comment_line(fid,'$');
               [DAT.timecode] = strread(rec,'%d',1);
               if OPT.debug(1)
                  OPT.disp;disp(['time ',num2str(DAT.timecode)])
               end
               rec = fgetl_no_comment_line(fid,'$');
            else
               DAT.timecode = 0;
            end
            %1 : ISO-notation 19870530.153000
            %2 : (as in HP compiler) ’30?May?87 15:30:00’
            %3 : (as in Lahey compiler) 05/30/87.15:30:00
            %4 : 15:30:00
            %5 : 87/05/30 15:30:00’
            %6 : as in WAM 8705301530
            
            if     DAT.timecode==1, DAT.timefmt = 'yyyymmdd.HHMMSS';
            elseif DAT.timecode==2, DAT.timefmt = 'dd-mmm-yy HH:MM:SS';
            elseif DAT.timecode==3, DAT.timefmt = 'mm-dd-yy HH:MM:SS';
            elseif DAT.timecode==4, DAT.timefmt = 'HH:MM:SS';
            elseif DAT.timecode==5, DAT.timefmt = 'yy/mm/dd HH:MM:SS';
            elseif DAT.timecode==6, DAT.timefmt = 'yymmddHHMM';
            end
   
%% Read ITER (optional, for test output)

            if strcmp(strtok(upper(rec)),'ITER')
   
               rec = fgetl_no_comment_line(fid,'$');
               [DAT.iter] = strread(rec,'%d');
               if OPT.debug(1)
                  disp(['iter ',num2str(DAT.iter)])
               end
               rec = fgetl_no_comment_line(fid,'$');
            else
               DAT.iter = [];
            end
            
%% Read # locations and coordinates
            
            if strcmp(strtok(upper(rec)),'LOCATIONS') | strcmp(strtok(upper(rec)),'LONLAT');
               
               if strcmp(strtok(upper(rec)),'LOCATIONS')
                  x = 'x';
                  y = 'y';
               end
               if strcmp(strtok(upper(rec)),'LONLAT');
                  x = 'lon';
                  y = 'lat';
               end

               rec = fgetl_no_comment_line(fid,'$');
              %[DAT.number_of_locations]     = sscanf(rec,'%i',1);
               [DAT.number_of_locations,rec] = strtok(rec);
               DAT.number_of_locations      = str2num(DAT.number_of_locations);
               if OPT.debug(1)
                 disp(['number_of_locations ',num2str(DAT.number_of_locations)])
               end         
   
              % read improvement in SWANmud version: mxc, myc
              [DAT.mxc,rec] = strtok(rec);
               DAT.mxc      = str2num(DAT.mxc);
              [DAT.myc,rec] = strtok(rec);
               DAT.myc      = str2num(DAT.myc);
               if isempty(DAT.mxc)
               DAT.mxc      = DAT.number_of_locations;
               DAT.myc      = 1;
               end
               if OPT.debug(1) & DAT.myc >1
                 disp(['number_of_locations mxc ',num2str(DAT.mxc)])
                 disp(['number_of_locations myc ',num2str(DAT.myc)])
               end         
               DAT.(x) = repmat(nan,[1 DAT.number_of_locations]);
               DAT.(y) = repmat(nan,[1 DAT.number_of_locations]);
               
               if ~OPT.fast
               % OLD: slow, but comments are allowed between numbers with this approach
                  for iloc=1:DAT.number_of_locations
                     rec         = fgetl_no_comment_line(fid,'$');
                     numbers     = sscanf(rec,'%f',2);
                     DAT.(x)(iloc) = numbers(1);
                     DAT.(y)(iloc) = numbers(2);
                     if OPT.debug(2)
                        if mod(iloc,OPT.mod)==0
                        disp(['   iloc,',x,',',y,': ',num2str([iloc,DAT.(x)(iloc),DAT.(y)(iloc)])]);
                        end
                     end         
                  end
               elseif OPT.fast
               %  NEW:fast, but no comments allowed between numbers
                  raw = fscanf(fid,'%f',2*DAT.number_of_locations);
                  DAT.(x) = raw(1:2:end);
                  DAT.(y) = raw(2:2:end);
               end % OPT.fast
               DAT.(x) = reshape(DAT.(x),DAT.myc, DAT.mxc);
               DAT.(y) = reshape(DAT.(y),DAT.myc, DAT.mxc);
               
               DAT.(x)(DAT.(x)==-99) = nan;
               DAT.(y)(DAT.(y)==-99) = nan;
               
               if OPT.debug(3)==1 & DAT.myc>1
                  TMP = figure;
                  pcolorcorcen(DAT.(x),DAT.(y),DAT.(x),[.5 .5 .5])
                  axis equal
                  pausedisp
                  try
                     close(TMP)
                  end
               end
               
               rec = fgetl_no_comment_line(fid,'$');
               
            else
            
               fclose(fid);
               error('either keyword LOCATIONS or LONLAT required.')
               
            end
            
%% Read # of frequencies

            if strcmp(strtok(upper(rec)),'AFREQ') | ...
               strcmp(strtok(upper(rec)),'RFREQ')
   
               if     strcmp(strtok(upper(rec)),'AFREQ')
                  DAT.frequency_type = 'absolute';
               elseif strcmp(strtok(upper(rec)),'RFREQ')
                  DAT.frequency_type = 'relative';
               end
   
               rec = fgetl_no_comment_line(fid,'$');
               [DAT.number_of_frequencies] = sscanf(rec,'%i',1);
               if OPT.debug(1)
                  disp(['number_of_frequencies ',num2str(DAT.number_of_frequencies)])
               end         
   
               DAT.frequency = repmat(nan,[1 DAT.number_of_frequencies]);
   
               for ifreq=1:DAT.number_of_frequencies
                  rec                  = fgetl_no_comment_line(fid,'$');
                  numbers              = sscanf(rec,'%e',2);
                  DAT.frequency(ifreq) = numbers(1);
               end
               
               rec = fgetl_no_comment_line(fid,'$');
   
            else
   
               % When no frequencies present: file contains integral parameters form TEST command
               DAT.dimension_of_spectrum = 0;
               DAT.number_of_frequencies = 0;
               
            end
            
%% Read # of directions (only for 2D spectra)
   
            if     strcmp(strtok(upper(rec)),'NDIR') 
               DAT.dimension_of_spectrum = 2;
               DAT.direction_convention  = 'nautical';
            elseif strcmp(strtok(upper(rec)),'CDIR') 
               DAT.dimension_of_spectrum = 2;
               DAT.direction_convention  = 'cartesian';
            elseif ~(isfield(DAT,'dimension_of_spectrum'))
            % see above: When no frequencies present: file contains integral parameters form TEST command
            % or MUDFILE (only frequencies, no directions)
               DAT.dimension_of_spectrum = 1;
            end
   
            if DAT.dimension_of_spectrum==2
               rec = fgetl_no_comment_line(fid,'$');
               [DAT.number_of_directions] = sscanf(rec,'%i',1);
               if OPT.debug(1)
                  disp(['number_of_directions ',num2str(DAT.number_of_directions)])
               end            
               if DAT.number_of_directions==0
                  fclose(fid);
                  error('Number of directions = 0. Remove *DIR keyword or use > 0 directions')
               end
   
               DAT.directions = repmat(nan,[1 DAT.number_of_directions]);
   
               for idir=1:DAT.number_of_directions
                  rec                  = fgetl_no_comment_line(fid,'$');
                  numbers              = sscanf(rec,'%e',2);
                  DAT.directions(idir) = numbers(1);
               end
               rec = fgetl_no_comment_line(fid,'$');
            end
   
%% Read # of quantities: fixed at 3 for 1D spectrum ,
%                     extended to 5 for MUD version.
%                        fixed at 1 for 2D spectrum.
%                                 2 for new MUDFile.

            if strcmp(strtok(upper(rec)),'QUANT')
   
               rec = fgetl_no_comment_line(fid,'$');
               [DAT.number_of_quantities] = sscanf(rec,'%i',1);
               if OPT.debug(1)
                  disp(['number_of_quantities ',num2str(DAT.number_of_quantities)])
               end
   
              %DAT.quantity_names            = not preallocated;
              %DAT.quantity_units            = not preallocated;
               DAT.quantity_exception_values = repmat(nan,[1 DAT.number_of_quantities]);
   
               for iquant=1:DAT.number_of_quantities
                  rec                                   = fgetl_no_comment_line(fid,'$');
                 [DAT.quantity_names{iquant},rec]       = strtok(rec);
                  DAT.quantity_names_long{iquant}       = strtrim(rec);
                 
   
                  rec                                   = fgetl_no_comment_line(fid,'$');
                  DAT.quantity_units{iquant}            = strtok(rec);
   
                  rec                                   = fgetl_no_comment_line(fid,'$');
                  numbers                               = sscanf(rec,'%e',2);
                  DAT.quantity_exception_values(iquant) = numbers(1);
               end
            end
   
%% CASE: read integrated test quantities S1D

            if DAT.timecode >0
               rec = fgetl_no_comment_line(fid,'$');
               DAT.time = datenum(strtok(rec),DAT.timefmt);
               DAT.time_ftell = ftell(fid); % breadcrumb
            end
            if DAT.dimension_of_spectrum==0 &  ~isempty(DAT.iter)
            
               while ~feof(fid)
   
                  % Read data block per location
                  for iloc=1:DAT.number_of_locations
                     if iloc==1
                        rec       = fgetl_no_comment_line(fid,'$');
                        iteration = str2num(strtok(rec));
                     end
                     
                         if feof(fid)
                             break
                         end
                         start_of_key_word = strfind(rec,'Test point');
                         shift = 10;
   
                         location_index = sscanf(rec(start_of_key_word+shift:end),'%i',1);
   
                         rawdata = fscanf(fid,'%e',  DAT.number_of_quantities);
   
                         % Split block into array per quantity
                         %  where the 1st dimension is the location
                         for j=1:DAT.number_of_quantities
                            quantity_name = char(DAT.quantity_names{j});
                            array         = rawdata(j);
                            array(array==DAT.quantity_exception_values(j))=nan;
                            DAT.(quantity_name)(iloc,iteration) = array;
                         end  
                  end
   
                  DAT.data_description = ['1st dimension = number_of_locations          ';
                                          '2nd dimension = number_of_iterations         '];
   
               end         
            
%% CASE: read test points S1D
   
            elseif DAT.dimension_of_spectrum==1 &  ~isempty(DAT.iter)
               
               while ~feof(fid)
               
                  % Read data block per location
                  for iloc=1:DAT.number_of_locations
                     if iloc==1
                        rec       = fgetl_no_comment_line(fid,'$');
                        iteration = str2num(strtok(rec));
                     end
                     
                         rec               = fgetl_no_comment_line(fid,'$');
                         if feof(fid)
                             break
                         end
                         start_of_key_word = strfind(rec,'Test point');
                         shift = 10;
   
                         location_index = sscanf(rec(start_of_key_word+shift:end),'%i',1);
   
                         rawdata = fscanf(fid,'%e',  DAT.number_of_quantities*...
                                                     DAT.number_of_frequencies);
   
                         rawdata  = reshape(rawdata,[DAT.number_of_quantities ...
                                                     DAT.number_of_frequencies]);
   
                         % Split block into array per quantity
                         %  where the 1st dimension is the location
                         for j=1:DAT.number_of_quantities
                            quantity_name = char(DAT.quantity_names{j});
                            array         = rawdata(j,:);
                            array(array==DAT.quantity_exception_values(j))=nan;
                            DAT.(quantity_name)(location_index,:,iteration) = array;
                         end  
                  end
                  DAT.data_description = ['1st dimension = number_of_locations          ';
                                          '2st dimension = number_of_frequencies        ';
                                          '3rd dimension = number_of_iterations         '];
               end
   
%% CASE: read S1D (incl. MUDFile)
   
            elseif DAT.dimension_of_spectrum==1
            
                % count times
                if DAT.timecode > 0
                    rec0 = rec;
                    pos_pnt = ftell(fid);
                    while ischar(rec)
                       for iloc=1:DAT.number_of_locations
                       rec = fgetl_no_comment_line(fid,'$');
                       rawdata  = fscanf(fid,'%e', DAT.number_of_quantities*...
                                                   DAT.number_of_frequencies);
                       end
                       rec = fgetl_no_comment_line(fid,'$');
                       if ischar(rec)
                       DAT.time(end+1)       = datenum(strtok(rec),DAT.timefmt);
                       DAT.time_ftell(end+1) = ftell(fid); % breadcrumb
                       end
                    end
                    fseek(fid, pos_pnt, 'bof'); % move of the pointer to the original position
                    rec = rec0;
                else
                    DAT.time = [nan];  % mind length(DAT.time) should be 1 for allocating arrays
                end
                
               % pre-allocate for speed (makes HUGE difference)
               for j=1:DAT.number_of_quantities
                  quantity_name       = char(DAT.quantity_names{j});
                  DAT.(quantity_name) = nan  ([DAT.number_of_locations ...
                                               DAT.number_of_frequencies ...
                                               length(DAT.time)]);                
               end
               
               for it=1:length(DAT.time)
               % Read data block per location
               for iloc=1:DAT.number_of_locations
               
                  if OPT.debug(2)
                     if mod(iloc,OPT.mod)==0
                     disp(['processing station ',num2str(iloc),' of ',num2str(DAT.number_of_locations)])
                     end
                  end
                  
                  rec = fgetl_no_comment_line(fid,'$');
                  
                  if strcmp(strtok(upper(rec)),'LOCATION')
                        
                        start_of_key_word = strfind(rec,'LOCATION');
                        shift             = 8;
                        
                        location_index    = sscanf(rec(start_of_key_word+shift:end),'%i',1);
                        
                        if ~(location_index == iloc)
                           fclose(fid);
                           error(['order of data in file does not match location order for line with ',rec])
                        end
   
                        if length(DAT.time) >1
                            rawdata  = fscanf(fid,'%e', DAT.number_of_quantities*...
                                                        DAT.number_of_frequencies);
                        else
                            % >>>> START OF modification by sfl, stylianos.flampouris@gmail.com, 16 Jan 2013
                            % The original script was importing the 1D-spectra with the fscanf as
                            % applied at the lines 534-535. Some times the values of EnDens are
                            % extremely small, e.g.  0.1991E-244. In order, the length of the output 
                            % string, to be kept constant, the "E" is neglected, so for the example the
                            % exported value is 0.1991-244. When this value is imported with fscanf, 
                            % matlab interpets them as two values 0.1991 and -244 which apparantly is a 
                            % mistake.
                            %
                            % In order to overcome the problem, the length of fscanf(fid,'%e',inf) is 
                            % examined, when the length of the data is longer than the number of 
                            % frequencies times the number of quantites the data are imported as
                            % strings and are splitted into substrings according to the position of the
                            % white spaces. The code of the split of the strings is not generic.
                            %
                            % TODO note this does not work YET when multiple times are in 1 file, 
                            % due to fscanf inf whcih scans also time of next block
                            %
                            % fscanf is used twice so it's necessary to know the possition of the 
                            % pointer before the use of fscanf for first time
                            
                            pos_pnt = ftell(fid);   
                            rawdata = fscanf(fid,'%e',inf);
                            if length(rawdata)==DAT.number_of_quantities*DAT.number_of_frequencies; % 1st time use of fscanf, now the pointer is at the end of the file
                               % bussiness as usual: proceed
                            else % when the previously described problem occurs
                                fseek(fid, pos_pnt, 'bof'); % move of the pointer to the original position
                                rawdata1 = fscanf(fid,'%100c', DAT.number_of_quantities*...
                                                        DAT.number_of_frequencies); % the 100 is just a big number of characters for each line, fscanf stops at each line when there are not more values.
                                rawdata=[];
                                stp=length(rawdata1)/DAT.number_of_frequencies;
                                spc_pos = strfind(rawdata1(1:stp),' '); % find the position of the white spaces

                                for ir=1:stp:length(rawdata1)

                                    if isempty(strfind(rawdata1(ir:ir+spc_pos(2)-1),'E'))
                                        % Building the number x.xxxxE-zzz and converting it into double  
                                        dum_str = [rawdata1(ir:ir+strfind(rawdata1(ir:ir+spc_pos(2)-1),'-')-2),...
                                            'E',...
                                        rawdata1(ir+strfind(rawdata1(ir:ir+spc_pos(2)-1),'-')-1:ir+spc_pos(2)-1)];
                                        rawdata =[rawdata;str2double(dum_str) ];      
                                    else
                                        rawdata =[rawdata;str2double(rawdata1(ir:ir+spc_pos(2)-1)) ];
                                    end
                                    rawdata =[rawdata;str2double(rawdata1(ir+spc_pos(2):ir+spc_pos(end)-1)) ]; %DIR
                                    rawdata =[rawdata;str2double(rawdata1(ir+spc_pos(end):ir+stp-1)) ]; %DirSPReaD
                                end 
                            end
                            % <<<< END OF modification, 16 Jan 2013
                        end

                        if isinf(rawdata)
                            fclose(fid);
                            error([DAT.filename,' contains Inf data for location ',num2str(iloc)])
                        else
                        rawdata  = reshape(rawdata,[DAT.number_of_quantities ...
                                                    DAT.number_of_frequencies]);
                        end
   	    
                        % Split block into array per quantity
                        %  where the 1st dimension is the location
                        for j=1:DAT.number_of_quantities
                           quantity_name               = char(DAT.quantity_names{j});
                           array                       = rawdata(j,:);
                           array(array==DAT.quantity_exception_values(j))=nan;
                           DAT.(quantity_name)(iloc,:,it) = array;
                        end
   
                  elseif strcmp(strtok(upper(rec)),'NODATA')
                  
                       % disp('keyword NODATA not tested yet !!!!!!!!!!!!!!!!!!!')
                        
                        nans = nan.*DAT.frequency;
                        for j=1:DAT.number_of_quantities
                           quantity_name               = char(DAT.quantity_names{j});
                           DAT.(quantity_name)(iloc,:,it) = nans;
                        end
   
                  else
                        
                        fclose(fid);
                        error(['unknown data in record: ',rec])
   
                  end   
                        
               end % for i=1:DAT.number_of_locations
               rec = fgetl_no_comment_line(fid,'$'); % skip time row
               end % time
               
               if DAT.myc==1
                  DAT.data_description = ['1st dimension = number_of_locations          ';
                                          '2st dimension = number_of_frequencies        '];
               elseif DAT.myc>1
               
               % reshape locations into 2D matrix if present
               
                  DAT.data_description = ['1st dimension = mxc                          ';
                                          '2nd dimension = myc                          ';
                                          '3rd dimension = number_of_frequencies        '];
                  for j=1:DAT.number_of_quantities
                     quantity_name = char(DAT.quantity_names{j});
                     DAT.(quantity_name) = reshape(DAT.(quantity_name),...
                                                  [DAT.myc DAT.mxc DAT.number_of_frequencies]);
                     % is this correctly reshaped
                  end  
               end
               
%% CASE: read test points S2D
   
            elseif DAT.dimension_of_spectrum==2 &  ~isempty(DAT.iter)
   
               warning('TEST S2D does not tested yet !!!')
               
               while ~feof(fid)
               
               % Read data block per location
               
               count = 0;
               
               for iloc=1:DAT.number_of_locations
   
               if OPT.debug(2)
                  if mod(iloc,OPT.mod)==0
                  disp(['processing station ',num2str(iloc),' of ',num2str(DAT.number_of_locations)])
                  end
               end
   
               for iquant=1:DAT.number_of_quantities
               
               count = count + 1;
   
                  if iloc==1
                     rec       = fgetl_no_comment_line(fid,'$');
                     iteration = str2num(strtok(rec));
                  end
                     
                  rec = fgetl_no_comment_line(fid,'$');
                  
                  quantity_name = char(DAT.quantity_names{iquant});
                  
                  if strcmp(strtok(upper(rec)),'FACTOR')
                     rec    = fgetl_no_comment_line(fid,'$');
                     factor = sscanf(rec,'%e',2);
                     
                     data  = fscanf(fid,'%e',DAT.number_of_frequencies*...
                                             DAT.number_of_directions);
                     data  = reshape(data,  [DAT.number_of_directions  DAT.number_of_frequencies]).*factor;
   
                     data(data==DAT.quantity_exception_values(iquant))=nan;
                     
                     DAT.(quantity_name)(iloc,:,:)   = data'; % note transpose
                     
                  elseif strcmp(strtok(upper(rec)),'ZERO') % for whole block of zeros
                        
                     DAT.(quantity_name)(iloc,:,:)   = zeros([DAT.number_of_frequencies...
                                                              DAT.number_of_directions ]);
                    %warning('keyword ZERO   not fully tested yet !!!!!!!!!!!!!!!!!!!')
                  
                  elseif strcmp(strtok(upper(rec)),'NODATA')
                        
                     DAT.(quantity_name)(iloc,:,:)   = nan  ([DAT.number_of_frequencies...
                                                              DAT.number_of_directions ]);
                    %warning('keyword NODATA not fully tested yet !!!!!!!!!!!!!!!!!!!')
   
                  end
                  
                  %pcolorcorcen(data)
                  %title([quantity_name,' @ loc ',num2str(iloc),' count: ',num2str(count)])
                  %pausedisp
            
               end % for i=1:DAT.number_of_quantities
               end % for i=1:DAT.number_of_locations
               end % while
               
               % TO DO: insert reshape (...,mxc, ...myc)
               DAT.data_description = ['1st dimension = number_of_locations          ';
                                       '2st dimension = number_of_frequencies        ';
                                       '3rd dimension = number_of_directions         ';
                                       'For plotting:                                ';
                                       'd = swan_spectrum(filename);                 ';
                                       'surf(d.frequencies,d.directions,d.data)% or  ';
                                       'pcolor(d.frequency,d.directions,d.data)    ';
                                       'shading interp                               ';
                                       'set(gca,''xscale'',''log'')                      ';
                                       'xlabel(''freq [Hz]'')                          ';
                                       'ylabel(''dir [\circ]'')                        '];
   
%% CASE: read full S2D (incl HOTFile)

            elseif DAT.dimension_of_spectrum==2
            
                % count times
                if DAT.timecode > 0
                    rec0 = rec;
                    pos_pnt = ftell(fid);
                    while ischar(rec)
                       for iloc=1:DAT.number_of_locations
                       rec = fgetl_no_comment_line(fid,'$'); % FACTOR, ZERO or NODATA
                       if strcmp(strtok(upper(rec)),'FACTOR')
                          rec    = fgetl_no_comment_line(fid,'$');
                          data  = fscanf(fid,'%e',DAT.number_of_frequencies.*...
                                                  DAT.number_of_directions);
                       end
                       end
                       rec = fgetl_no_comment_line(fid,'$');
                       if ischar(rec)
                       DAT.time(end+1)       = datenum(strtok(rec),DAT.timefmt);
                       DAT.time_ftell(end+1) = ftell(fid); % breadcrumb
                       end
                    end
                    fseek(fid, pos_pnt, 'bof'); % move of the pointer to the original position
                    rec = rec0;
                else
                    DAT.time = [nan];  % mind length(DAT.time) should be 1 for allocating arrays
                end
            
               quantity_name = DAT.quantity_names{1}; % for 2D there can be only 1 quantity
               %% pre-allocate for speed (makes HUGE difference)
               DAT.(quantity_name) = nan  ([DAT.number_of_locations ...
                                            DAT.number_of_frequencies ...
                                            DAT.number_of_directions  ...
                                            length(DAT.time)]);
               for it=1:length(DAT.time)
               % Read data block per location
               for iloc=1:DAT.number_of_locations
               
                  if OPT.debug(2)
                     if mod(iloc,OPT.mod)==0
                     disp(['processing station ',num2str(iloc),' of ',num2str(DAT.number_of_locations)])
                     end
                  end
                  
                  % SWAN writes in following loop
                  %  DO 290 IX = 1, MXC
                  %    DO 280 IY = 1, MYC
                  %      INDX = KGRPNT(IX,IY)
                  %      IF (INDX.EQ.1) THEN
                  %        WRITE (NREF,220) 'NODATA'                                     40.08
                  %      ELSE
                  %        DO IS = 1, MSC
                  %          WRITE (NREF, FIX_SPEC) (NINT(ACLOC(ID,IS)/EFAC), ID=1,MDC)      40.13
                  %        ENDDO
                  %      ENDIF
                  %    ENDDO
                  %  ENDDO
           
                  rec = fgetl_no_comment_line(fid,'$');
                  if strcmp(strtok(upper(rec)),'FACTOR')
                     rec    = fgetl_no_comment_line(fid,'$');
                     factor = sscanf(rec,'%e',2);
                     
                     data  = fscanf(fid,'%e',DAT.number_of_frequencies.*...
                                             DAT.number_of_directions);
%                     DAT.number_of_frequencies
%                     DAT.number_of_directions 
%                    length(data)
                   data  = reshape(data  ,[DAT.number_of_directions DAT.number_of_frequencies]).*factor;
              
                   DAT.(quantity_name)(iloc,:,:,it)   = data'; % note transpose
                     
                  elseif strcmp(strtok(upper(rec)),'ZERO')
                        
                     DAT.(quantity_name)(iloc,:,:,it) = zeros([DAT.number_of_frequencies ...
                                                               DAT.number_of_directions]);
                  
                  elseif strcmp(strtok(upper(rec)),'NODATA')
                        
                     DAT.(quantity_name)(iloc,:,:,it) = nan  ([DAT.number_of_frequencies ...
                                                               DAT.number_of_directions]);
   
                  else
                        
                     fclose(fid);
                     error(['unknown data in record: ',rec])
                     
                  end
            
               end % for i=1:DAT.number_of_locations
               rec = fgetl_no_comment_line(fid,'$'); % skip time row
               end % time
               
               if DAT.myc==1
                  DAT.data_description = char({'1st dimension = number_of_locations          ',
                                          '2st dimension = number_of_frequencies        ',
                                          '3rd dimension = number_of_directions         ',
                                          'For plotting:                                ',
                                          'd = swan_spectrum(filename);                 ',
                                          'surf(d.frequency,d.directions,d.data)% or  ',
                                          'pcolor(d.frequency,d.directions,d.data)    ',
                                          'shading interp                               ',
                                          'set(gca,''xscale'',''log'')                      ',
                                          'xlabel(''freq [Hz]'')                          ',
                                          'ylabel(''dir [\circ]'')                        '});
               elseif DAT.myc>1
               
               % reshape locations into 2D matrix if present
               
                  DAT.data_description = char({'1st dimension = mxc                          ',
                                          '2nd dimension = myc                          ',
                                          '3rd dimension = number_of_frequencies        ',
                                          '4th dimension = number_of_directions         ',
                                          'For plotting:                                ',
                                          'd = swan_spectrum(filename);                 ',
                                          'surf(d.frequency,d.directions,d.data)% or  ',
                                          'pcolor(d.frequency,d.directions,d.data)    ',
                                          'shading interp                               ',
                                          'set(gca,''xscale'',''log'')                      ',
                                          'xlabel(''freq [Hz]'')                          ',
                                          'ylabel(''dir [\circ]'')                        '});

                  DAT.(quantity_name) = reshape(DAT.(quantity_name),...
                                               [DAT.myc DAT.mxc DAT.number_of_frequencies DAT.number_of_directions]);
               end % if DAT.myc==1          
   
            end % if DAT.dimension_of_spectrum
   
         %catch
         % 
         %   if nargout==1
         %      error(['Error reading file: ',fname])
         %   else
         %      iostat = -3;
         %   end      
         %
         %end % try
   
         fclose(fid);
   
      end % if fid <0
   
   end % if length(tmp)==0
   
end % if iostat (GUI)
   
DAT.read.with     = 'read by $Id$ by G.J. de Boer (WL | Delft Hydraulics)';
DAT.read.at       = datestr(now);
DAT.read.iostatus = iostat;

%% Function output

if nargout      ==0 | nargout==1
   varargout= {DAT};
elseif nargout==2
   varargout= {DAT, iostat};
end

%% EOF
