function hh = myquiver3(p,v,normals,headsize,headratio,fore,aft,varargin)
%MYQUIVER3 3-D quiver plot.
%   This differs from the standard quiver3 in the following ways:
%   1.  The origins and velocities are given as row vectors, or matrices of
%   the same. 
%   2.  It takes an additional argument NORMALS, which should follow the
%   velocity argument.  The arrowheads will be drawn in a plane perpendicular
%   to the corresponding normal vector.
%   3.  It takes an additional argument HEADSIZE specifying the size of the
%   arrowheads as a proportion of the size of the shaft.  The handle
%   returned by quiver3 has a MaxHeadSize attribute, but I find there
%   appears to be a maximum size that it allows, irrespective of the value
%   requested.
%   4.  It takes an additional argument AXISHANDLE.
%   5.  It does not return a quiver object.  Instead, it returns an N*3
%   array of handles to line objects.  The three handles in each row are
%   the shaft and two barbs of a single arrow.

  % fprintf( 1, 'myquiver3\n' );

    error(nargchk(2,inf,nargin,'struct'));
    if (nargin < 3) || isempty(normals)
        normals = findPerpVector( v );
    end
    if (nargin < 4) || isempty(headsize), headsize = 0.5; end
    if (nargin < 5) || isempty(headratio), headratio = 0.6; end
    if (nargin < 6) || isempty(fore), fore = 1; end
    if (nargin < 7) || isempty(aft), aft = 1; end

    ends = p+fore*v;
    starts = p-aft*v;

    h1 = line( [starts(:,1), ends(:,1)]', ...
               [starts(:,2), ends(:,2)]', ...
               [starts(:,3), ends(:,3)]', ...
               varargin{:} );

    if headsize == 0
        if nargout>0, hh = h1; end
    else
        headlength = headsize*(fore+aft);
        barbstart = ends - v*headlength;
        barboffset = cross( v, normals, 2 )*(headlength*headratio);  % Depends on normals being unit vectors.
        h2 = line( [barbstart(:,1) + barboffset(:,1), ends(:,1)]', ...
                   [barbstart(:,2) + barboffset(:,2), ends(:,2)]', ...
                   [barbstart(:,3) + barboffset(:,3), ends(:,3)]', ...
                   varargin{:} );

        h3 = line( [barbstart(:,1) - barboffset(:,1), ends(:,1)]', ...
                   [barbstart(:,2) - barboffset(:,2), ends(:,2)]', ...
                   [barbstart(:,3) - barboffset(:,3), ends(:,3)]', ...
                   varargin{:} );
        if nargout>0, hh = [h1,h2,h3]; end
    end
end
