function [emiout, beta, gamma, xrms, xprms] = ppemitcalc(xin, xpin, intin, energyin, massin, xory)
%
% function [emiout, beta, gamma, xrms, xprms] = ppemitcalc(x, xp, intense, energy, mass, xory)
%
%    PPEMITCALC - Calculate beam emittance for Pepperpot data from X, X'
%    and intensity values
%
%    ppemitcalc(x, xp, intense) - displays the RMS normalised emittance of
%    the pepperpot image from the specified X and X' values of each
%    particle.  The array INTENSE gives the intensity information at each
%    point; the X, X' and INTENSE arrays can be either 1D or 2D arrays, but
%    must be the same size, each referring to a point on the pepperpot
%    image.  The X values are assumed to be in [mm] and the X' values in
%    [mrad] - the result is therefore given in [pi mm mrad].  The particles
%    are assumed to be H-minus ions travelling at 35keV.
%
%    [emiout] = ppemitcalc(x, xp, intense) - as above, but the emittance
%    value is output to the variable EMIOUT.
%
%    [emiout] = ppemitcalc(x, xp, intense, energy) - also specifies the
%    beam energy in eV.  The ENERGY array can either be a single number,
%    specifying the mean beam energy, or an array of energies for
%    individual particles, the same size as X and XP.  The default value is
%    35 keV.
%
%    [emiout] = ppemitcalc(x, xp, intense, energy, mass) - also specifies
%    the mass of the particles in the beam, in units of the atomic mass.
%    The MASS array can either be a single number, specifying the mass of
%    the particle in question, or an array of masses for individual
%    particles, the same size as X and XP.  The default value is the mass
%    of the H-minus.
%
%    [emiout] = ppemitcalc(x, xp, intense, energy, mass, xory) - makes it
%    possible to specify whether one is inputting X or Y data. XORY must be
%    a string: either 'X' or 'Y'.  This makes no difference to the
%    emittance calculation itself, only the results displayed.  The default
%    value is 'X'.
%
%    [emiout, beta, gamma] = ppemitcalc(x, xp, intense) - returns the emittance of the
%    beam along with the mean Beta and Gamma values (relativistic, not
%    Twiss).
%
%    [emiout, beta, gamma, xrms, xprms] = ppemitcalc(x, xp, intense) - returns the rms
%    X and X' values.

    if ( nargin < 3 ) | ( nargin > 6 )
        error('Incorrect number of input variables: syntax is [emiout] = ppemitcalc(x, xp, intense, energy, mass, xory)') ;
    elseif nargout > 5
        error('Too many output variables: syntax is [emiout, beta, gamma, xrms, xprms] = ppemitcalc(x, xp, intense, energy, mass, xory)') ;
    end

%    xsize = size(x) ;
%    xpsize = size(xp) ;
%    intsize = size(intense) ;

    if ( size(xin) ~= size(xpin) )
        error('PPEmitCalc:incorrectSyntax', ...
            ['Input variables x and xp are different sizes: must contain same number '...
            '\n of data points ie. pairs of x and xp values']) ;
    elseif ( size(intin) ~= size(xin) )
        error('Intensity array must be of equal size to x and xp arrays') ;
    end

%    Check for beam parameters; default is H-minus at 35keV

    c  = 299792458 ;        % Speed of light in m/s
    ma = 1.66053886e-27 ;   % Atomic mass in kg
    m0 = 1.0083736 * ma ;   % Mass of H-minus
    qe = -1.60217653e-19 ;  % Electron charge in C
    q0 = qe ;               % H-minus is singularly charged
    E0 = 35.0e3 ;           % Beam energy is 35 keV after ion source [eV]

%    Check if variables have been input, otherwise use default, and
%    calculate energy spread

    if nargin < 6 || isempty(xory)
        xory = ['X'] ;
    end

    if ~ischar(xory)
        error('xory variable must be a string') ;
    end

    if nargin < 5 || isempty(massin)
        massc = m0 ;
    else
        masssize = size(massin) ;
        mass = reshape(massin,[],1) ;
        massc = mass.*ma ;
        if ( masssize ~= size(xin) ) & ( masssize(1) ~= 1 ) & ( masssize(2) ~= 1 )
            error('Mass variable must either be single number or an array of equal size to x and xp arrays') ;
        end
    end

    if nargin < 4 || isempty(energyin)
        energy = E0 ;
        dE = 0 ;
    else
        ensize = size(energyin) ;
        energy = reshape(energyin,[],1) ;
        meanE = mean(energy) ;
        maxE = max(energy) ; minE = min(energy) ;
        dE = (maxE - minE)./meanE ;    % Energy spread [fraction of E0]
        if ( size(energyin) ~= size(xin) ) & ( ensize(1) ~= 1 ) & ( ensize(2) ~= 1 )
            error('Energy variable must either be single number or an array of equal size to x and xp arrays') ;
        end
    end

%    Calculate velocities and energy spreads

    Gam  = 1 - ((qe.*energy)./(massc.*(c.^2))) ;    % Corresponding Lorentz factor Gamma []
    BG   = sqrt( (Gam.^2) - 1 ) ;                   % Corresponding beta-gamma factor []
    Beta = sqrt( 1 - (Gam.^-2) ) ;                  % Corresponding Normalized velocity []
    dG = ( Gam - 1 ).*dE ;                          % Gamma spread [fraction of Gam]

%    Initialise variables for sum data

    xSum = 0 ; xpSum = 0 ; emittance = 0 ;
    x2Sum = 0 ; xp2Sum = 0 ; xxpSum = 0 ;
    xMean = 0 ; xpMean = 0 ;

    intense = reshape(intin,[],1) ;
    numpart = sum(intense) ;
    intgood = intense > 0 ;

%    Generate sums and mean values

    x = reshape(xin,numel(xin),1) ;
    xp = reshape(xpin,numel(xpin),1) ;

    xMean = sum(x.*intense)./numpart ;
    xpMean = sum(xp.*intense)./numpart ;

    xSum = sum( (x - xMean).*intense ) ;
    xpSum = sum( (xp - xpMean).*intense ) ;
    x2Sum = sum( ((x - xMean).^2).*intense ) ;
    xp2Sum = sum( ((xp - xpMean).^2).*intense ) ;
    xxpSum = sum( (x - xMean).*(xp - xpMean).*intense ) ;
    x2 = x2Sum./numpart ;
    xp2 = xp2Sum./numpart ;
    xxp = xxpSum./numpart ;

    emittance = sqrt( (xp2.*x2) - (xxp.^2) ) ;
    xrms = sqrt(x2) ;
    xprms = sqrt(xp2) ;
    xmax = max( x - xMean ) ;
    xdiam = max( ( x - xMean ).*intgood ) - min( ( x - xMean ).*intgood ) ;
    xpmax = max( xp - xpMean ) ;
    xpdiam = max( ( xp - xpMean ).*intgood ) - min( ( xp - xpMean ).*intgood ) ;

%    Calculate beta and gamma to normalise the emittance

    beta = mean(Beta) ;
    gamma = mean(Gam) ;
    betagamma = mean(BG) ;

%    Normalise the emittance

    emiout = (emittance.*beta.*gamma) ;

%    Calculate Twiss parameters and single particle emittances

    Tbeta = x2./emittance ;
    Tgamma = xp2./emittance ;
    Talpha = -sqrt((Tbeta.*Tgamma) - 1) ;

%    Calculate covariance matrix and give Alpha the correct sign

    covmat = cov(x,xp,1) ;
    if covmat(2) < 0
        Talpha = -Talpha ;
    end

%    Recalculate emittance

    emirecalc = ((Tgamma.*x2) + (2.*Talpha.*xxp) + (Tbeta.*xp2))./2 ;
    emirecalcnorm = (emirecalc.*beta.*gamma) ;
%    emiarrun = ((Tgamma.*x2arr) + (2.*Talpha.*xxparr) + (Tbeta.*xp2arr)) ;
%    emiarr = (emiarrun.*beta.*gamma) ;
%    emisort = sort(emiarr) ;
%    emi90 = emisort(round(length(emisort).*0.9)) ;
%    emi100 = max(emiarr) ;

%    Display all data

    disp(' ')
    disp(['Normalised RMS ' xory ' emittance is: ' num2str(emiout,8) ' [pi mm mrad].'])
    disp(' ')
    disp(['Un-normalised RMS ' xory ' emittance is: ' num2str(emittance,8) ' [pi mm mrad].'])
    disp(['Mean Beta-Gamma is: ' num2str(betagamma,8) '.'])
    disp(['RMS ' xory ' size is: ' num2str(xrms,8) ' [mm].'])
    disp(['RMS ' xory ''' size is: ' num2str(xprms,8) ' [mrad].'])
    disp(['Total ' xory ' width is: ' num2str(xdiam,8) ' [mm].'])
    disp(['Total ' xory ''' width is: ' num2str(xpdiam,8) ' [mrad].'])
    disp(['Ratio of Total to RMS size (' xory ') is : ' num2str((xdiam./xrms),8) '.'])
    disp(['Ratio of Total to RMS size (' xory ''') is : ' num2str((xpdiam./xprms),8) '.'])
    disp(' ')
    disp(['Courant-Snyder (Twiss) Parameters: '])
    disp(['Alpha = ' num2str(Talpha) ' ; Beta = ' num2str(Tbeta) ' ; Gamma = ' num2str(Tgamma) ' .'])
    disp(' ')
%    disp(['Recalculated Normalised RMS emittance is: ' num2str(emirecalcnorm,8) ' [pi mm mrad].'])
%    disp(' ')

%    disp(['Normalised 90% emittance is: ' num2str(emi90,8) ' [pi mm mrad].'])
%    disp(['Normalised 100% emittance is: ' num2str(emi100,8) ' [pi mm mrad].'])
%    disp(' ')

    return

