function is_within = angle_within( a, a_min, a_max, epsilon )

  if nargin < 3
    error( [ mfilename ' needs at least 3 input arguments.' ] );
  end
  
  if ~exist( 'epsilon', 'var' )
    epsilon = 1e-10;
  end
  
  % Check the sizes

  tmp = sortrows( unique( [ size( a ); size( a_min ); size( a_max ) ], 'rows' ) );
  
  if size( tmp, 1 ) > 2
    error( [ mfilename ': inconsistent input argument sizes (each argument must be either [ 1 1 ] or desired output size).' ] );
  end

  if size( tmp, 1 ) > 1

    if ~all( tmp( 1,: ) == [ 1 1 ] )
    error( [ mfilename ': inconsistent input argument sizes (each argument must be either [ 1 1 ] or desired output size).' ] );

    else
      desired_size = tmp( 2,: );
      
      % Replicate scalars if necessary
      if numel( a ) == 1,	a     = repmat( a,     desired_size );      end
      if numel( a_min ) == 1,   a_min = repmat( a_min, desired_size );      end
      if numel( a_max ) == 1,   a_max = repmat( a_max, desired_size );      end

    end

  end
  

  % Do the actual comparison of angles:
  % Does "a" lie within the arc from "a_min" to "a_max"?
  
  is_within = ( mod( a - a_min + epsilon, 2 * pi ) < mod( 2 * epsilon + a_max - a_min, 2 * pi ) );
  