function [histmat, xout, yout, custmap] = hist2dplot(data, xbins, ybins, xmin, xmax, ymin, ymax)
%
% function [histdat, xout, yout, custmap] = hist2dplot([xdat, ydat], xbins, ybins, xmin, xmax, ymin, ymax)
%
%    HIST2DPLOT - Create and plot 2 dimensional histogram
%
%    hist2dplot([xdat, ydat]) - plots a 2D histogram with 32 horizontal and
%    vertical bins and automatically calculates the horizontal and vertical
%    ranges based on the input dataset.  XDAT contains the X position and
%    YDAT the Y position for the data to be histogrammed.  X, Y and Z label
%    data for the plot must be added afterwards.
%
%    [histdat] = hist2dplot([xdat, ydat]) - returns the 2D matrix
%    containing the histogram information.
%
%    [histdat, xout, yout] = hist2dplot([xdat, ydat]) - also returns the X
%    and Y plot data along with the histogram information.
%
%    [histdat, xout, yout, custmap] = hist2dplot([xdat, ydat]) - also
%    returns the custom colour map used to create the shading on the
%    histogram plots, based on 'jet(64)'.
%
%    hist2dplot([xdat, ydat], xbins, ybins) - specifies the number of
%    horizontal and vertical bins for the histogram.  The default is 32.
%    The horizontal and vertical range is calculated automatically.
%
%    hist2dplot([xdat, ydat], xbins, ybins, xmin, xmax, ymin, ymax) -
%    specifies the horizontal and vertical ranges along with the number of
%    horizontal and vertical bins.  The default is to let HIST2DPLOT
%    calculate the ranges automatically based upon the input data.
%
%    The HIST2DPLOT function uses the SURFC function to plot the data,
%    since this provides the best colourmap and interpolation facility.
%    However, any of the 3D plotting functions can be used (with the
%    exception of BAR3) to plot the data, since the syntax and input
%    structure is the same.

    if ( nargin > 7 )
        error('Too many input variables: syntax is hist2dplot([xdat, ydat], xbins, ybins, xmin, xmax, ymin, ymax)') ;
    elseif ( nargin == 4 ) | ( nargin == 5 ) | ( nargin == 6 )
        error('Must specify X and Y min and max for histogram: syntax is hist2dplot([xdat, ydat], xbins, ybins, xmin, xmax, ymin, ymax)') ;
    elseif ( nargin == 2 )
        error('Must specify number of bins in X and Y: syntax is hist2dplot([xdat, ydat], xbins, ybins)') ;
    elseif ( nargout > 4 )
        error('Too many output variables: syntax is [histdat, xout, yout, custmap] = hist2dplot([xdat, ydat])') ;
    elseif ( nargout == 2 )
        error('Must specify both X and Y data: syntax is [histdat, xout, yout] = hist2dplot([xdat, ydat])') ;
    end

    [nrow, ncol] = size(data) ;

    if ncol ~= 2
        error('Input data must specify X and Y data in two columns') ;
    end

    xdat = data(:,1) ; ydat = data(:,2) ;

    if ( nargin < 7 )
        ymax = max(ydat) ;
    end
    if ( nargin < 6 )
        ymin = min(ydat) ;
    end
    if ( nargin < 5 )
        xmax = max(xdat) ;
    end
    if ( nargin < 4 )
        xmin = min(xdat) ; 
    end
    if ( nargin < 3 )
        xbins = 32 ; ybins = 32 ;
    end

    xdiff = (xmax - xmin)./(xbins - 1) ;
    ydiff = (ymax - ymin)./(ybins - 1) ;
    xdiff2 = xdiff./2 ; ydiff2 = ydiff./2 ;
    xedges = linspace( xmin-xdiff2, xmax+xdiff2, xbins + 1 ) ;
    yedges = linspace( ymin-ydiff2, ymax+ydiff2, ybins + 1 ) ;
    xgrid = linspace( xmin, xmax, xbins ) ;
    ygrid = linspace( ymin, ymax, ybins ) ;
%    xedges = linspace( xmin, xmax, (xbins + 1) ) ;
%    yedges = linspace( ymin, ymax, (ybins + 1) ) ;
%    xdiff = ( ( xedges(2) - xedges(1) ) / 2 ) ;
%    ydiff = ( ( yedges(2) - yedges(1) ) / 2 ) ;
%    xgrid = ( xedges + xdiff ) ; xgrid(end) = [] ;
%    ygrid = ( yedges + ydiff ) ; ygrid(end) = [] ;

    histmat = zeros( length(ygrid), length(xgrid) ) ;

    for yloop = 1:ybins

        ylo = yedges( yloop ) ;
        yhi = yedges( yloop + 1 ) ;

        if yloop == ybins
            [yidx] = find ( (ydat >= ylo) & (ydat <= yhi) ) ;
        else
            [yidx] = find ( (ydat >= ylo) & (ydat < yhi) ) ;
        end

        numx = xdat(yidx) ;

        if ( ~isempty(numx) )

            xhist = histc( numx, xedges ) ;
            xhist(end-1) = xhist(end-1) + xhist(end) ;
            xhist(end) = [] ;
            histmat(yloop,:) = xhist ;

        end

    end

    custmap = jet(64) ;
    custmap(1,:) = [1 1 1] ;
    surfc(xgrid,ygrid,histmat), axis tight, shading interp ; view(0,90) ;
%    contourf(xgrid,ygrid,histmat,16) ;
    colormap(custmap) ; hold on ;
    xlines = [min(xdat) ymin max(xdat) ymin; ...
              min(xdat) ymax max(xdat) ymax] ;
    ylines = [xmin min(ydat) xmin max(ydat); ...
              xmax min(ydat) xmax max(ydat)] ;
    plot(xlines(:,1),xlines(:,2),'-k') ;
    plot(xlines(:,3),xlines(:,4),'-k') ;
    plot(ylines(:,1),ylines(:,2),'-k') ;
    plot(ylines(:,3),ylines(:,4),'-k') ;
    hold off ;
    
    xout = xgrid ; yout = ygrid ;

    return
