function [emithist] = ppemithistgen(x, y, xp, yp, intens, dx, dy, dxp, dyp)
%
% function [emithist] = ppemithistgen(x, y, xp, yp, intens, dx, dy, dxp, dyp)
%
%    PPEMITHISTGEN - generate emittance plot from pepperpot data

    xnew = [] ; ynew = [] ; xpnew = [] ; ypnew = [] ;
%    intens = round(intens) ;
%    cumint = cumsum(intens) ;
    cumint = round(cumsum(intens)) ;
    lowint = [0; cumint(1:end-1)] ;
    intens = cumint - lowint ;

    if size(dx) == 1
        dx = dx.*(ones(length(x),1)) ;
    end

    if size(dy) == 1
        dy = dy.*(ones(length(x),1)) ;
    end

    if size(dxp) == 1
        dxp = dxp.*(ones(length(x),1)) ;
    end

    if size(dyp) == 1
        dyp = dyp.*(ones(length(x),1)) ;
    end

    for i = 1:length(x)
        if intens(i) > 0
            if i == 1
                xnew(1:cumint(i)) = x(i) - (dx(i)./2) + (rand([intens(i) 1]).*dx(i)) ;
                ynew(1:cumint(i)) = y(i) - (dy(i)./2) + (rand([intens(i) 1]).*dy(i)) ;
                xpnew(1:cumint(i)) = xp(i) - (dxp(i)./2) + (rand([intens(i) 1]).*dxp(i)) ;
                ypnew(1:cumint(i)) = yp(i) - (dyp(i)./2) + (rand([intens(i) 1]).*dyp(i)) ;
%                xpnew(1:cumint(i)) = xp(i) + (randn([intens(i) 1]).*dxp(i)) ;
%                ypnew(1:cumint(i)) = yp(i) + (randn([intens(i) 1]).*dyp(i)) ;
            else
                xnew((cumint(i-1) + 1):cumint(i)) = x(i) - (dx(i)./2) + (rand([intens(i) 1]).*dx(i)) ;
                ynew((cumint(i-1) + 1):cumint(i)) = y(i) - (dy(i)./2) + (rand([intens(i) 1]).*dy(i)) ;
                xpnew((cumint(i-1) + 1):cumint(i)) = xp(i) - (dxp(i)./2) + (rand([intens(i) 1]).*dxp(i)) ;
                ypnew((cumint(i-1) + 1):cumint(i)) = yp(i) - (dyp(i)./2) + (rand([intens(i) 1]).*dyp(i)) ;
%                xpnew((cumint(i-1) + 1):cumint(i)) = xp(i) + (randn([intens(i) 1]).*dxp(i)) ;
%                ypnew((cumint(i-1) + 1):cumint(i)) = yp(i) + (randn([intens(i) 1]).*dyp(i)) ;
            end
        end
    end

    emithist  = [xnew' ynew' xpnew' ypnew'] ;

    return