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
     if kk < 1
        disp(['No files found in ' dirname ' with filespec ' filespec ]);
        return
     end
     % sort filenames alphabetically
     fnames=sort(fnames);
  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));

  % remove height outliers
  
  hell=res(:,4);
  mhell2=median(hell);
  ok=find(hell > mhell2-10 & hell < mhell2+10 );
  hell(hell <= mhell2-10)=mhell2-10;  
  hell(hell >= mhell2+10)=mhell2+10;  
  mhell=mean(hell(ok));
  shell=std(hell(ok));
  notok=1:length(hell);
  notok(ok)=[];
    
  % 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);
  accept=find(res(:,5)==5);
  plot(res(accept,1)/86400,hell(accept),'g.')
  ylabel('hell [m]')
  mhella=mean(hell(accept));
  shella=std(hell(accept));
  legend(sprintf('height (mean=%3.2f m, std=%3.2f m, FQ=FLOAT)',mhella,shella))
  datetick('x',15)
  title(plottitle)
  subplot(3,1,2)
  %plot(res(:,1)/86400,res(:,4)-mhell,'b.')
  plot(res(ok,1)/86400,hell(ok),'b.')
  %ylabel('\Deltahell [m]')
  ylabel('hell [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))
  hold on
  plot(res(notok,1)/86400,hell(notok),'rx')
  plot(res(accept,1)/86400,hell(accept),'g.')
  a=axis;;a(3:4)=[mhell2-10.3 mhell2+10.3];axis(a);
  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),40);axis(a);
  legend('nsat','10*hdop','fq','\Deltat(DGPS)')
  ylabel('nsat hdop*10 fq \Deltat [s]')
  datetick('x',15)



%  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);
  accept=find(res(:,5)==5);

  % analyze data with acceptable solutions
  t=res(accept,1);
  h=hell(accept);
  lat=res(accept,2);
  lon=res(accept,3);
  % find the gaps
  dt=diff(t);
  interv=[find(dt > 30*60) ; length(t) ];
  i2=0;
  for ii=1:length(interv)
     i1=i2+1;
     i2=interv(ii);
     mh=mean(h(i1:i2));
     sh=std(h(i1:i2));
     mlat=mean(lat(i1:i2))*180/pi;
     mlon=mean(lon(i1:i2))*180/pi;
     mt=mean(t(i1:i2));
     fprintf('%7.0f %4d-%4d  %4d  %9.4f %9.4f  %7.2f %7.2f\n',mt,i1,i2,i2-i1+1,mlat,mlon,mh,sh)
  end

end

return

