function varargout=nmearead(filespec)
%NMEAREAD  Read NMEA files.
%   RES=NMEAREAD(FILESPEC) reads a NMEA file FILESPEC and returns the data
%   of the GPGGA records as an array RES. FILESPEC can be a string with
%   a single filename, or FILESPEC can be a directory with wildcards or FILESPEC
%   can be a cell array of strings with filenames. In case of multiple files
%   files are read sequentially. Output is an array RES with 8 columns, with in
%   column 1 the time TK in seconds since the start of the day, in columns 2-4
%   the lattitude (rad), longitude (rad) and ellipsoidal height (m), and in 
%   columns 5-7 the fix quality, number of satellites involved in the fix and
%   horizontal dillution, and in column 8, if applicable, the latency of the
%   DGPS corrections.
%
%   NMEAREAD(FILESPEC) reads a NMEA file FILESPEC, without output argument,
%   plots the data.
%
%   Examples:
%       nmearead('XPNMEA.GGA')                      % Single file, plot data
%       nmearead('i:/offline/fugro/nmea/*XP*')      % Wildcards, plot data
%       res=nmearead('i:/offline/fugro/nmea/*XP*'); % Wildcard, no plot
%       res=nmearead({'FILE1.GGA', 'FILE2.GGA'});   % Cell array with filenames
%
%   See also NMEACHECKSUM and NMEAGPGGA.
%
%   Copyright 2010, TU Delft, DEOS.

%   Created:    4 March 2010 by Hans van der Marel
%   Modified:   

% if input is a string then convert to cell array , also set plot title

if ~iscell(filespec)
  if exist(filespec) == 2
     % filespec is a file, convert to cell array  
     fnames={filespec};
  else
     % filespec could be a directory with wildcards
     dirstruct=dir(filespec);
     dirname=fileparts(filespec);
     kk=0;
     for k=1:length(dirstruct)
       if dirstruct(k).isdir ==0
          kk=kk+1; 
          fullname=fullfile(dirname,dirstruct(k).name);
          fnames(kk)={fullname};     
       end
     end
  end
  plottitle=filespec;
else
   % filespec is a cell array, copy to fnames
   fnames=filespec;
   plottitle=[ fnames{1} ' - ' fnames{end} ];
end
plottitle=strrep(plottitle,'\','/');

% initialize the output array, pre-allocate storage and set counter (j) 

j=0;
res=zeros(1000,8);

% loop over the files

for k=1:length(fnames)

  % Get the filename

  fname=fnames{k};

  % read NMEA messages from file

  fid=fopen(fname);
  while 1
    tline = fgetl(fid);
    if ~ischar(tline), break, end
    % disp(tline)
    msgtype=tline(1:6);
    % process the NMEA messages
    switch msgtype
      case '$GPGGA' 
         [tk,plh,fq,nsat,hdop,dgpst]=nmeagpgga(tline);
        % store the results
         if j >= size(res,1)
          % preallocate storage
          res=[res ; zeros(1000,8) ];   
        end
        j=j+1;
        res(j,1:8)=[tk plh fq nsat hdop dgpst];
      otherwise
        disp(['Unsupported NMEA message type' msgtype ]);
    end
  end
  fclose(fid);

  % done reading file

end

% done reading all files, trim the results array  res
% Trim the data arrays

res(j+1:end,:)=[];

% check the number of output arguments, if no output arguments make a plot

if nargout > 0

  varargout={res};

else

  % make a plot 
  
  % compute mean latitude, longitude and height, and standard deviations

  mlat=mean(res(:,2));
  mlon=mean(res(:,3));
  mhell=mean(res(:,4));

  flat=6370000;
  flon=flat*cos(mlat);

  slat=std(flat*(res(:,2)-mlat));
  slon=std(flon*(res(:,3)-mlon));
  shell=std(res(:,4));

  % make the plot

  subplot(3,1,1)
  plot(res(:,1)/86400,flat*(res(:,2)-mlat),'b',res(:,1)/86400,flon*(res(:,3)-mlon),'r')
  legend(sprintf('lattitude (mean=%6.4f deg, std=%3.2f m)',mlat*180/pi,slat),sprintf('longitude (mean=%6.4f deg, std=%3.2f m)',mlon*180/pi,slon))
  ylabel('\Delta\phi \Delta\lambda [m]')
  %a=axis;ay=max(abs(a(3)),abs(a(4)));a(3:4)=[-ay ay];axis(a);
  datetick('x',15)
  title(plottitle)
  subplot(3,1,2)
  plot(res(:,1)/86400,res(:,4)-mhell)
  ylabel('\Deltahell [m]')
  %a=axis;ay=max(abs(a(3)),abs(a(4)));a(3:4)=[-ay ay];axis(a);
  legend(sprintf('height (mean=%3.2f m, std=%3.2f m)',mhell,shell))
  datetick('x',15)
  subplot(3,1,3)
  h=plot(res(:,1)/86400,[res(:,6) res(:,7)*10 res(:,5) res(:,8)]);
  set(h(1),'LineWidth',2);
  set(h(2),'LineWidth',2);
  hold on;
  h=plot(res(:,1)/86400,[res(:,6) res(:,7)*10 ]);
  set(h(1),'LineWidth',2);
  set(h(2),'LineWidth',2);
  hold off;
  a=axis;a(3)=0;a(4)=min(a(4),30);axis(a);
  legend('nsat','10*hdop','fq','\Deltat(DGPS)')
  ylabel('nsat hdop*10 fq \Deltat [s]')
  datetick('x',15)

end

return

