function varargout = statupd(data, st, varargin)
%  1.  NAME
%
%  statupd - Update statistics array
%
%  2.  SYNOPSIS
%
%  [upd, hdg] = statupd(data, stat, options )
%
%  3.  DESCRIPTION
%
%  Udates the current statistics array with the given subset of data. See
%  stat(5m) for more information.
%
%  4.  PARAMETERS
%
%     data
%        Matrix (MxN) of N data for the M different variables.
%
%     stat
%        Satitistics array (Mx7).
%
%  5.  OPTIONS
%
%     '-hdg', hdg
%        Non-NaN elements in vector hdg of length M specifiy the last
%        value of heading,which is to be used to calculate the mean
%        heading. This allow to calculate heading mean values accross the
%        360/0 point.  Updated heading statistics is converted within 0
%        .. 360.
%
%  6.  RETURNS
%
%     upd
%        Statistics array updated with the data vector.
%
%     hdg optional
%        Array with last value of heading,to be used for successive calls
%        to this function. This output argument should be used if option
%        -hdg is used.  .nr ll -1
%
%  7.  SEE ALSO
%
%  flagdata(3m)stat(5m)statcalc(3m)statinit(3m)statlim(3m)stat_str(3m)statupd(3m)
%
%  8.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

%--- Defaults:
d = data;
nvars = size(d,1);
if nvars ~= size(st,1),
  error([ 'Statistic and data arrays do not have the same number of' ...
		  ' columns' ]);
end
hdg = ones(nvars,1) * NaN;

%--- Options:
l = nargin - 2;
if l > 0,
  i = 1;
  while i <= l,
	 tmp = varargin{i};
	 if isstr(tmp),
		switch tmp
		  case '-hdg',
			 hdg = varargin{i+1};
			 i = i + 2;
		  otherwise,
			 error([ 'Option ' tmp 'not implemented' ]); 
		end
	 else
		% Process other options
	 end
  end
end

%--- Calculate relative heading:
ihdg = find(~isnan(hdg));
if ~isempty(ihdg),
  d(ihdg,:) = next_hdg(hdg(ihdg), d(ihdg,:));
end

%--- Modify heading for extrema calculations:
if ~isempty(ihdg),
  data(ihdg,:) = mod(data(ihdg,:), 360);
end

%--- Indexes to stat array:
min_ = 1; max_ = 2; imin_ = 3; imax_ = 4; mn_ = 5; sd_ = 6; n_ = 7;

%--- Update extrema:
ii = 1:nvars;
for i=ii,
  jm = min(data(i,:));
  if jm < st(i, min_),
	 j = find(data(i,:) <= jm);
	 st(i, min_) = jm;
	 st(i, imin_) = j(1);
  end
  jm = max(data(i,:));
  if jm > st(i, max_),
	 j = find(data(i,:) >= jm);
	 st(i, max_) = jm;
	 st(i, imax_) = j(1);
  end
end

%--- Update statistics:
for i=ii,
  good = find(~isnan(d(i,:)));
  if ~isempty(good),
	 st(i,mn_) = st(i,mn_) + sum(d(i,good));
	 st(i,sd_) = st(i,sd_) + sum( d(i,good) .* d(i,good) );
	 st(i, n_) = st(i,n_) + length(good);
  end
end

if nargout > 0,
  varargout{1} = st;
end
if nargout > 1,
  hdg = ones(nvars,1) * NaN;
  if ~isempty(ihdg),
	 n = size(hdg,2);
	 hdg(ihdg) = d(ihdg, n);
  end
  varargout{2} = hdg;
end

return;
