 function [xi,yi,hc] = maxr(a)
%MAXR   Find interpolated maximizer(s) and max value(s)
%       for (each column of) a.
%
%               [xi,yi,hc] = maxr(a)
%
%  Calls max() followed by qint() for quadratic interpolation.
%
   [m,n] = size(a);
   if m==1, a=a'; t=m; m=n; n=t; end;
   [y,x] = max(a);
   xi=x;    % vector of maximizer locations, one per col of a
   yi=y;    % vector of maximum values, one per column of a
   if nargout>2, hc = zeros(1,n); end
   for j=1:n,   % loop over columns
     if x(j)>1, % only support interior maxima
       if x(j)<m,
         [xdelta,yij,cj] = qint(a(x(j)-1,j),y(j),a(x(j)+1,j));
         xi(j) = x(j) + xdelta;
         if nargout>2, hc(j) = cj; end
         if (nargout>1), yi(j) = yij; end
       end;
     end;
   end;