function rel = next_hdg(ref, hdg)
%  1.  NAME
%
%  next_hdg - Get nearest heading value
%
%  2.  SYNOPSIS
%
%  rel = next_hdg(ref, hdg )
%
%  3.  DESCRIPTION
%
%  This Matlab function calculates nearest heading value relative to the
%  reference heading value ref. This function is required when
%  calculating statistics of heading values that cross the boundary
%  {360,0}. Its vailidity depends on the fact that the difference between
%  both heading values is less than +-180, since the relative heading is
%  given by the reference heading plus the heading change in the
%  direction of the shortest angle difference.
%
%  Input parameters may be vectors of the same length to convert several
%  heading values at once.
%
%  If heading values are stored in an array, one dimension must be of the
%  same length as the reference values. The other dimension is assumed to
%  be a list off successive heading values. Each of them is recalculated
%  relatively to the previous value, the first value being calculated
%  relatively to its corresponding reference element.
%
%  4.  PARAMETERS
%
%     ref
%        Value for reference heading.
%
%     hdg
%        Value of current heading.
%
%  5.  RETURNS
%
%     new
%        The heading calculated relatively to the reference.
%
%  6.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

% Check length:
if length(ref) ~= size(hdg,1),
  error('Both arguments must have the same length.');
end
% Check size of ref:
if size(ref,1) < size(ref,2),
  ref = ref';
end

% Convert angles to a 0..360 range:
hdg = mod(hdg, 360);
rel = ones(size(hdg)) * Inf;

%--- Get relative heading:
ii = 1:size(hdg,2);
for i=ii,
  d1 = hdg(:,i) - ref(:);
  d2 = 360 - abs(d1(:));
  % Select angle change in direction of shortest difference:
  ih = find( (abs(d1) <= d2) & (rel(:,i) == Inf));
  if ~isempty(ih),
	 rel(ih(:),i) = ref(ih(:)) + d1(ih(:));
  end
  ih = find( (d1 <= 0) & (rel(:,i) == Inf));
  if ~isempty(ih),
	 rel(ih(:),i) = ref(ih(:)) + d2(ih(:));
  end
  ih = find(rel(:,i) == Inf);
  if ~isempty(ih),
	 rel(ih(:),i) = ref(ih(:)) - d2(ih(:));
  end
  ref = rel(:,i);
end

return;

