function [xnew, xpnew, ynew, ypnew] = emitcor(xdat, ydat, xaxis, yaxis, xyhist)
%
% function [xnew, xpnew, ynew, ypnew] = emitcor(xdat, ydat, xaxis, yaxis, xyhist)
%
%    EMITCOR - Create correlated X-Y beam distribution based on X-X' and
%    Y-Y' emittance plots.
%
%    function [xnew, xpnew, ynew, ypnew] = emitcor(xdat, ydat, xaxis, yaxis, xyhist)
%
%    Using the EMITCONV function creates 1D arrays of position and angle
%    data in each plane, with each row corresponding to X-X' (or Y-Y') data
%    for a single particle.  However, since X and Y data must be created
%    seperately, there is no X-Y correlation, other than that randomly
%    produced by the ordering of the entries in the arrays.
%
%    EMITCOR takes pairs of X-X' and Y-Y' values and matches them together
%    to create a correlated X-Y distribution, based upon a histogram of X-Y
%    data.  The input arguments are as follow:
%
%        XDAT - 2D array of size (N x 2) with horizontal position (column
%        1) and angle (column 2) data for each particle;
%        YDAT - 2D array of size (N x 2) with vertical position (column
%        1) and angle (column 2) data for each particle;
%        XAXIS - 1D array of length (L) with X position data for each
%        column of the histogram data;
%        YAXIS - 1D array of length (M) with Y position data for each 
%        row of the histogram data;
%        XYHIST - an (M x L) 2D array of histogram data ie. number of
%        particles for each phase space cell.  The size of XYHIST must
%        match the lengths of XDAT and YDAT.
%
%    Given that the data in XDAT and YDAT should correspond to the same
%    number of particles for which one wishes to create a distribution, the
%    sizes of XDAT and YDAT must match.  The output arguments are as
%    follows:
%
%        XNEW - 1D array of size (N) with horizontal position data for each
%        particle;
%        XPNEW - 1D array of size (N) with horizontal angle data for each
%        particle;
%        YNEW - 1D array of size (N) with vertical position data for each
%        particle;
%        YPNEW - 1D array of size (N) with vertical angle data for each
%        particle;
%
%    Each row now corresponds to 4D data for a single particle, where (N)
%    is the number of particles.

    if nargin ~= 5
        error('Incorrect number of input variables: syntax is [xnew, xpnew, ynew, ypnew] = emitcor(xdat, ydat, xaxis, yaxis, xyhist)') ;
    elseif nargout ~= 4
        error('Incorrect number of output variables: syntax is [xnew, xpnew, ynew, ypnew] = emitcor(xdat, ydat, xaxis, yaxis, xyhist)') ;
    elseif ( ndims(xdat) ~= 2 ) | ( ndims(ydat) ~= 2 )
        error('Incorrect number of dimensions for X and Y data') ;
    end

    x = xdat(:,1) ; xp = xdat(:,2) ;
    y = ydat(:,1) ; yp = ydat(:,2) ;
    numpart = length(x) ;

    histsize = size(xyhist) ;
    datsize = [length(yaxis) length(xaxis)] ;
    numhist = sum(sum(xyhist)) ;

    if ( length(x) ~= length(y) )
        error('X and Y data are of different lengths - must be of equal size to correspond to the same particle set') ;
    elseif ( histsize(1) ~= datsize(1) ) | ( histsize(2) ~= datsize(2) )
        error('The sizes of XDAT and YDAT must match the histogram dimensions') ;
    end

%    Create new particle positions

    xt = [] ; xpt = [] ;
    yt = [] ; ypt = [] ;
    xl = x ; xpl = xp ;
    yl = y ; ypl = yp ;
    xrange = [] ; yrange = [] ;

    xdif = (max(xaxis) - min(xaxis))./(length(xaxis) - 1) ;
    ydif = (max(yaxis) - min(yaxis))./(length(yaxis) - 1) ;
    xmin = xaxis - (xdif./2) ; xmax = xaxis + (xdif./2) ;
    ymin = yaxis - (ydif./2) ; ymax = yaxis + (ydif./2) ;

    fracdist = xyhist./numhist ;
%    fracdist = (xyhist./numhist).*0.8 ;
    partsum = 0 ;
    ysumdist = 0 ;

    for i = 1:length(yaxis)

        xsumdist = 0 ;

        for j = 1:length(xaxis)

            xsumdist = sum(fracdist(i,1:j)) + ysumdist ;
            partsumdist = round(numpart.*xsumdist) ;
            reqparts = partsumdist - partsum ;

            if ( reqparts ~= 0 )

                xlo = xmin(j) ; xhi = xmax(j) ;
                ylo = ymin(i) ; yhi = ymax(i) ;

                if j == 1
                    [xidx] = find( (xl > -inf) & (xl < xhi), reqparts ) ;
                elseif i == length(xaxis)
                    [xidx] = find( (xl >= xlo) & (xl < inf), reqparts ) ;
                else
                    [xidx] = find( (xl >= xlo) & (xl < xhi), reqparts ) ;
                end

                if i == 1
                    [yidx] = find( (yl > -inf) & (yl < yhi), reqparts ) ;
                elseif i == length(yaxis)
                    [yidx] = find( (yl >= ylo) & (yl < inf), reqparts ) ;
                else
                    [yidx] = find( (yl >= ylo) & (yl < yhi), reqparts ) ;
                end

                if length(yidx) < length(xidx)
                    xidx((length(yidx) + 1):end) = [] ;
                elseif length(xidx) < length(yidx)
                    yidx((length(xidx) + 1):end) = [] ;
                end

                xex = xl(xidx) ; xpex = xpl(xidx) ;
                yex = yl(yidx) ; ypex = ypl(yidx) ;
                xt = [xt xex'] ; xpt = [xpt xpex'] ;
                yt = [yt yex'] ; ypt = [ypt ypex'] ;
                xl(xidx) = [] ; xpl(xidx) = [] ;
                yl(yidx) = [] ; ypl(yidx) = [] ;

            end

            partsum = partsum + reqparts ;

        end

        ysumdist = xsumdist ;

    end

    leftpart = length(xl) ;
    disp(['Number of left over particles = ' num2str(leftpart) ]) ;
    disp(' ') ;
    remordidx = randperm(leftpart) ;
    ylrp = yl(remordidx) ;
    yplrp = ypl(remordidx) ;
    xt = [xt xl'] ; xpt = [xpt xpl'] ;
    yt = [yt ylrp'] ; ypt = [ypt yplrp'] ;

    xnew = xt' ; xpnew = xpt' ;
    ynew = yt' ; ypnew = ypt' ;

    return

