function [x, xp] = emitconv(xdat, ydat, histdat, newpart)
%
% function [x, xp] = emitconv(xdat, ydat, histdat, newpart)
%
%    EMITCONV - Convert histogrammed emittance information into arrays of X
%    and X' information.
%
%    [x, xp] = emitconv(xdat, ydat, histdat) - converts histogrammed or
%    measured emittance data into arrays of X and X' data.  No assumption
%    is made about the units of the emittance data, since the units stay
%    the same when the new arrays are produced.  In order to produce the
%    necessary arrays, the input arrays must have the following structure:
%
%        XDAT - 1D array of length (M) with position data for each column
%        of the histogram data;
%        YDAT - 1D array of length (N) with angle or velocity data for each
%        row of the histogram data;
%        HISTDAT - an (N x M) 2D array of histogram data ie. number of
%        particles for each phase space cell.  The size of HISTDAT must
%        match the lengths of XDAT and YDAT.
%
%    The new particle dataset is created by first creating a linear
%    distribution of random numbers for each phase space cell, using the
%    RAND command.  Each of these numbers then corresponds to one
%    coordinate of a particle in the new dataset.
%
%    [x, xp] = emitconv(xdat, ydat, histdat, newpart) - generates a new set
%    of particles with NEWPART particles.  This allows a real dataset to be
%    expanded/reduced depending on the number of desired particles for the
%    new data, while retaining the same particle distribution.

    if ( nargin < 3 ) | ( nargin > 4 )
        error('Incorrect number of input variables: syntax is [x, xp] = emitconv(xdat, ydat, histdat, newpart)') ;
    elseif nargout ~= 2
        error('Incorrect number of output variables: syntax is [x, xp] = emitconv(xdat, ydat, histdat, newpart)') ;
    end

    numhist = sum(sum(histdat)) ;
    numpart = round(numhist) ;

    if nargin < 4
        newpart = numpart ;
    end

    histsize = size(histdat) ;
    datsize = [length(ydat) length(xdat)] ;
    
    if ( histsize(1) ~= datsize(1) ) | ( histsize(2) ~= datsize(2) )
        error('The sizes of XDAT and YDAT must match the histogram dimensions') ;
    end

%    Calculate old emittance value to compare with the new value

    xnumsum = sum(histdat) ; ynumsum = sum(histdat') ;

    xoweight = xdat.*xnumsum ;
    yoweight = ydat.*ynumsum ;
    xmean = sum(xoweight)./numhist ;
    ymean = sum(yoweight)./numhist ;
    xshift = xdat - xmean ;
    yshift = ydat - ymean ;
    xweight = xshift.*xnumsum ;
    yweight = yshift.*ynumsum ;
    xsum = sum(xweight) ;
    ysum = sum(yweight) ;
    x2sum = sum((xshift.^2).*xnumsum) ;
    y2sum = sum((yshift.^2).*ynumsum) ;
    x2 = x2sum./numhist ;
    y2 = y2sum./numhist ;
    xypl = yshift'*xshift ;
    xyweight = xypl.*histdat ;
    xysum = sum(sum(xyweight)) ;
    xy = xysum./numhist ;
    emitold = sqrt((x2.*y2) - xy.^2) ;

    disp(' ')
    disp(['Old Unnormalised emittance is ' num2str(emitold) ' [pi mm mrad].'])

%    Create new particle positions

    xt = [] ; xpt = [] ;
    xrange = [] ; yrange = [] ;

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

    fracdist = histdat./numhist ;
    partsum = 0 ;
    ysumdist = 0 ;

    for i = 1:length(ydat)

        xsumdist = 0 ;

        for j = 1:length(xdat)

            xsumdist = sum(fracdist(i,1:j)) + ysumdist ;
            partsumdist = round(newpart.*xsumdist) ;
            parts = partsumdist - partsum ;
%            parts = round(newpart.*fracdist(i,j)) ;
            xrange = xmin(j) + rand(parts,1).*xdif ;
            yrange = ymin(i) + rand(parts,1).*ydif ;
            xt = [xt xrange'] ;
            xpt = [xpt yrange'] ;
            partsum = partsum + parts ;

        end

        ysumdist = xsumdist ;

    end

    x = xt' ; xp = xpt' ;

    emitcalc(x, xp) ;

    return

