function [emiout, beta, gamma, xrms, xprms, twiss, emierror, emierrorfrac] = emitcalc(xin, xpin, intin, energyin, massin, xory, xerrin, xperrin, interrin)
%
% function [emiout, beta, gamma, xrms, xprms, twiss, emierror, emierrorfrac] = 
%            emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr, interr)
%
%    EMITCALC - Calculate beam emittance and errors from X and X' values
%
%    emitcalc(x, xp)
%    emitcalc(x, xp, intense)
%    emitcalc(x, xp, intense, energy)
%    emitcalc(x, xp, intense, energy, mass)
%    emitcalc(x, xp, intense, energy, mass, xory)
%    emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr)
%    emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr, interr)
%
%    [emiout] = emitcalc(...)
%    [emiout, beta] = emitcalc(...)
%    [emiout, beta, gamma] = emitcalc(...)
%    [emiout, beta, gamma, xrms] = emitcalc(...)
%    [emiout, beta, gamma, xrms, xprms] = emitcalc(...)
%    [emiout, beta, gamma, xrms, xprms, twiss] = emitcalc(...)
%    [emiout, beta, gamma, xrms, xprms, twiss, emierror] = emitcalc(...)
%    [emiout, beta, gamma, xrms, xprms, twiss, emierror, emierrorfrac] = emitcalc(...)
%
%    emitcalc(x, xp) - displays the RMS normalised emittance of the beam
%    from the specified X and X' values of each particle.  The X and X'
%    arrays must be the same size.  The X values are assumed to be in [mm]
%    and the X' values in [mrad] - the result is therefore given in [pi mm mrad].
%    For the normalisation, the particles are assumed to be H-minus ions
%    travelling at 35keV.
%
%    emitcalc(x, xp, intense) - also specifies the intensity, or particle
%    density, at the point [X, XP] with the array INTENSE.  INTENSE array
%    must either be the same size and shape as X and XP, or a single value
%    - usually 1 - to indicate that each value is equally weighted (this
%    would correspond to single particle data).
%
%    emitcalc(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.
%
%    emitcalc(x, xp, intense, energy, mass) - also specifies the mass of
%    the particles in the beam in eV.  The MASS array can either be a
%    single number, specifying the mass of all identical particles in the
%    beam, or an array of masses for individual particles of the same size
%    as X and XP.  The default value is the mass of the H-minus.
%
%    emitcalc(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'.
%
%    emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr) - includes
%    the calculation of emittance errors, based on the position error,
%    XERR, and the angle error, XPERR.  XERR and XPERR can either be single
%    values, if the error on each measurement is the same, or arrays of the
%    same size as X and XP.  XERR and XPERR are ABSOLUTE NOT PERCENTAGE
%    errors, so must be given in units of [mm] and [mrad] respectively.
%    Both arrays must be specified to calculate errors.
%
%    emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr, interr) -
%    also includes intensity errors in the emittance calculation.  INTERR
%    can either be a single value, if the error on each intensity
%    measurement is the same, or an array of the same size as INTENSE, and
%    must have the same units as INTENSE.
%
%    [emiout] = emitcalc(...) - outputs the normalised RMS emittance to the
%    variable EMIOUT.
%
%    [emiout, beta, gamma] = emitcalc(...) - also outputs the mean Beta and
%    Gamma values of the beam (relativistic, not Twiss).
%
%    [emiout, beta, gamma, xrms, xprms] = emitcalc(...) - outputs the rms X
%    and X' sizes of the beam.
%
%    [emiout, beta, gamma, xrms, xprms, twiss] = emitcalc(...) - outputs
%    the Twiss parameters of the beam as the structure array TWISS:
%    TWISS.BETA gives BETA_X, TWISS.ALPHA gives ALPHA_X and TWISS.GAMMA
%    gives GAMMA_X.  These are for the UNNORMALISED beam.
%
%    [emiout, beta, gamma, xrms, xprms, twiss, emierror] = emitcalc(...) -
%    outputs the absolute error on the normalised RMS emittance.  Only
%    works if XERR and XPERR (and optionally INTERR) are specified.
%
%    [emiout, beta, gamma, xrms, xprms, twiss, emierror, emierrorfrac] = emitcalc(...)
%        - outputs the fractional emittance error in addition to the
%        absolute error.

% File is the property of Simon Jolly, Imperial College London.  File
% released under the GNU public license.
% For information on the theory, see DIPAC'09 paper WEOA02.
%
% File history
%
%   17-May-2008 S. Jolly
%       File modified to calculate emittance errors, based on DIPAC paper
%       1st draft.  Tested with FETS pepperpot data.
%
%======================================================================

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

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

%    Check if variables have been input, otherwise use defaults, and
%    calculate beta, gamma and energy spread

    if ( nargin < 2 ) || ( nargin > 9 )
        error('Incorrect number of input variables: syntax is emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr, interr)') ;
    elseif nargout > 8
        error('Too many output variables: syntax is [emiout, beta, gamma, xrms, xprms, twiss, emierror, emierrorfrac] = emitcalc(...)') ;
    end

    x = reshape(xin,[],1) ;
    xp = reshape(xpin,[],1) ;
    numvals = length(x) ;
    onemat = ones(numvals,1) ;

    if ( length(xp) ~= numvals )
        error('EmitCalc:incorrectSyntax', ...
            ['Input variables x and xp are different sizes: must contain same number '...
            '\n of data points ie. pairs of x and xp values']) ;
    end

    if ( nargin < 9 ) || ( isempty(interrin) )
        interrin = 0 ;
    end

    if ( nargin == 7 )
        error('Must specify X and XP errors in pairs: emitcalc(x, xp, intense, energy, mass, xory, xerr, xperr)') ;
    end

    if ( nargin < 7 )
        xerrin = 0 ;
        xperrin = 0 ;
    elseif ( isempty(xerrin) ) || ( isempty(xperrin) )
        if ( isempty(xerrin) )
            xerrin = 0 ;
        end
        if ( isempty(xperrin) )
            xperrin = 0 ;
        end
    end

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

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

    if ( nargin < 5 ) || ( isempty(massin) )
        mass = m0 ;
    elseif ( ischar(massin) )
        switch massin
            case 'H0'
                mass = maev ;
            case 'H-'
                mass = m0 ;
            case 'e'
                mass = 510998.9 ;
            otherwise
                warning('Unknown mass selected: defaulting to H- mass') ;
                mass = m0 ;
        end
    else
        mass = reshape(massin,[],1) ;
    end

    if ~( length(mass) == 1 ) && ~( length(mass) == numvals )
        error('MASS variable must either be single valued or the same length as X and XP') ;
    end

    if ( nargin < 4 ) || ( isempty(energyin) )
        energy = E0 ;
    else
        energy = reshape(energyin,[],1) ;
    end

    if ~( length(energy) == 1 ) && ~( length(energy) == numvals )
        error('ENERGY variable must either be single valued or the same length as X and XP') ;
    else
        meanE = mean(energy) ;
        maxE = max(energy) ; minE = min(energy) ;
        dE = (maxE - minE)./meanE ;    % Energy spread [fraction of E0]
    end

    if ( nargin < 3 ) || ( isempty(intin) )
        intense = onemat ;
    elseif ( length(intin) == 1 )
        intense = onemat.*intin ;
    else
        intense = reshape(intin,[],1) ;
    end

    if ( length(intense) ~= numvals )
        error('INTENSE variable must either be single valued or the same length as X and XP') ;
    end

%    Calculate velocities and energy spreads

    Gam  = 1 + (energy./mass) ;         % 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 ;

    numpart = sum(intense) ;
    intgood = intense > 0 ;

%    Generate sums and mean values

    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 ;
    xarr = x - xMean ;
    xparr = xp - xpMean ;
    x2arr = (x - xMean).^2 ;
    xp2arr = (xp - xpMean).^2 ;
    xxparr = xarr.*xparr ;

    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) ;
    twiss.alpha = Talpha ;
    twiss.beta = Tbeta ;
    twiss.gamma = Tgamma ;

%    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) ;

%    Calculate errors if required

    if ( nargin > 6 )
        [emierrorunnorm, emierrorfrac] = emiterrors(x, xp, xerrin, xperrin, intense, interrin) ;
        emierror = ( emierrorunnorm.*beta.*gamma ) ;
        emierrortext = ([' +/- ' num2str(emierror,8)]) ;
    else
        emierror = 0 ;
        emierrorfrac = 0 ;
        emierrortext = [''] ;
    end

%    Display all data

    disp(' ')
    disp(['Normalised RMS ' xory ' emittance is: ' num2str(emiout,8) emierrortext ' [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(['Normalised 90% ' xory ' emittance is: ' num2str(emi90,8) ' [pi mm mrad].'])
%    disp(['Normalised 100% ' xory ' emittance is: ' num2str(emi100,8) ' [pi mm mrad].'])
    disp(' ')

    return

