function [histmat, xout, yout, custmap] = ppemithistplot(xdat, ydat, intens, xbins, ybins, xmin, xmax, ymin, ymax, ptype)
%
% function [histdat, xout, yout, custmap] = ppemithistplot(xdat, ydat, intens, xbins, ybins, xmin, xmax, ymin, ymax, ptype)
%
%    PPEMITHISTPLOT - Create and plot 2 dimensional histogram of emittance
%    data.
%
%    ppemithistplot(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] = ppemithistplot(xdat, ydat) - returns the 2D matrix
%    containing the histogram information.
%
%    [histdat, xout, yout] = ppemithistplot(xdat, ydat) - also returns the
%    X and Y plot data along with the histogram information.
%
%    [histdat, xout, yout, custmap] = ppemithistplot(xdat, ydat) - also
%    returns the custom colour map used to create the shading on the
%    histogram plots, based on 'jet(64)'.
%
%    ppemithistplot(xdat, ydat, intens) - provides intensity data for each
%    of the x and y positions.  The default is 1.
%
%    ppemithistplot(xdat, ydat, intens, 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.
%
%    ppemithistplot(xdat, ydat, intens, 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 PPEMITHISTPLOT
%    calculate the ranges automatically based upon the input data.
%
%    It is possible to leave any of these last 7 input variables
%    unspecified, such as:
%
%        ppemithistplot(xdat, ydat, intens, xbins, [], xmin, xmax, [], [])
%
%    PPEMITHISTPLOT will then use the default values for any empty
%    variables.  The example above would allow one to specify all the X
%    histogramming data without specifying the Y variables.  The default
%    values are:
%
%        INTENS: 1 (one particle per XDAT/YDAT entry)
%        XBINS: 41 (one bin per row/column of holes in the tungsten mesh)
%        YBINS: 32
%        XMIN: -60 (maximum distance of bottom/left hole from centre)
%        XMAX: 60 (maximum distance of top/right hole from centre)
%        YMIN: minimum value in YDAT dataset
%        YMAX: maximum value in YDAT dataset
%
%    The PPEMITHISTPLOT function uses the CONTOURF 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.  To specify the plot type, use:
%
%        ppemithistplot( ... , ptype)
%
%    ...where PTYPE can take the following values:
%
%        PTYPE = 1: plot data using 16 level CONTOURF (Default)
%        PTYPE = 2: plot data using 10 level CONTOURF
%        PTYPE = 3: plot data using SURFC
%        PTYPE = 4: plot data as image pixels using IMAGESC

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

    if ( nargin < 10 ) || isempty(ptype)
        ptype = 1 ;
    end
    if ( nargin < 9 ) || isempty(ymax)
        ymax = max(ydat) ;
    end
    if ( nargin < 8 ) || isempty(ymin)
        ymin = min(ydat) ;
    end
    if ( nargin < 7 ) || isempty(xmax)
%        xmax = max(xdat) ;
        xmax = 60 ;
    end
    if ( nargin < 6 ) || isempty(xmin)
%        xmin = min(xdat) ;
        xmin = -60 ;
    end
    if ( nargin < 5 ) || isempty(ybins)
%        ybins = 32 ;
        ybins = 16 ;
    end
    if ( nargin < 5 ) || isempty(xbins)
        xbins = 41 ;
    end
    if ( nargin < 3 ) || isempty(intens)
        intens = ones(length(xdat),1) ;
    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 ) ;

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

    [iidx] = find ( intens > 0 ) ;
    xdat = xdat ;
    ydat = ydat ;

    for yloop = 1:ybins

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

        for xloop = 1:xbins

            xlo = xedges( xloop ) ;
            xhi = xedges( xloop + 1 ) ;

            for datloop = 1:length(xdat)

                if yloop == ybins

                    if xloop == xbins

                        if ( (ydat(datloop) >= ylo) && (ydat(datloop) <= yhi) ) && ( (xdat(datloop) >= xlo) && (xdat(datloop) <= xhi) )
                            histmat(yloop,xloop) = histmat(yloop,xloop) + intens(datloop) ;
                        end

                    else

                        if ( (ydat(datloop) >= ylo) && (ydat(datloop) <= yhi) ) && ( (xdat(datloop) >= xlo) && (xdat(datloop) < xhi) )
                            histmat(yloop,xloop) = histmat(yloop,xloop) + intens(datloop) ;
                        end

                    end

                else

                    if xloop == xbins

                        if ( (ydat(datloop) >= ylo) && (ydat(datloop) < yhi) ) && ( (xdat(datloop) >= xlo) && (xdat(datloop) <= xhi) )
                            histmat(yloop,xloop) = histmat(yloop,xloop) + intens(datloop) ;
                        end

                    else

                        if ( (ydat(datloop) >= ylo) && (ydat(datloop) < yhi) ) && ( (xdat(datloop) >= xlo) && (xdat(datloop) < xhi) )
                            histmat(yloop,xloop) = histmat(yloop,xloop) + intens(datloop) ;
                        end

                    end

                end

            end

        end

    end

    custmap = jet(64) ;
    custmap(1,:) = [1 1 1] ;

    switch ptype
        case 1
            contourf(xgrid,ygrid,histmat,16), axis tight ;
        case 2
            contourf(xgrid,ygrid,histmat), axis tight ;
        case 3
            surfc(xgrid,ygrid,histmat), axis tight, shading interp ; view(0,90) ;
        case 4
            imagesc([xmin xmax],[ymin ymax],histmat,[0 max(max(histmat))]), axis tight ;
        otherwise
            contourf(xgrid,ygrid,histmat,16), axis tight ;
    end

    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
