function [emierror, emierrorfrac] = emiterrors(xin, xpin, xerrin, xperrin, intin, interrin, dispvars)
%
% function [emierror, emierrorfrac] = emiterrors(x, xp, xerr, xperr, intense, interr)
%
%    EMITERRORS - Calculate errors on unnormalised rms emittance
%
%    emiterrors(x, xp, xerr, xperr)
%    emiterrors(x, xp, xerr, xperr, intense)
%    emiterrors(x, xp, xerr, xperr, intense, interr)
%    emiterrors(x, xp, xerr, xperr, intense, interr, dispvars)
%
%    [emierror] = emiterrors(...)
%    [emierror, emierrorfrac] = emiterrors(...)
%
%    emiterrors(x, xp, xerr, xperr) - calculate absolute errors on
%    emittance for list of particles with positions X, angles XP, position
%    errors XERR and angle errors XPERR.  The X and XP arrays must be the
%    same size.  The X and XERR values are assumed to be in [mm] and the XP
%    and XPERR values in [mrad] - the result is therefore given in [pi mm
%    mrad].  The XERR and XPERR arrays can either be single values,
%    indicating an equal error on all measurements, or arrays of the same
%    size as X and XP.  Empty arrays for either value would indicate zero
%    error: this is identical to setting both values to zero.
%
%    Errors are quoted for UNNORMALISED EMITTANCE: error value must be
%    multiplied by (BetaGamma) to give normalised rms emittance error.
%
%    emiterrors(x, xp, xerr, xperr, intense) - also specifies the
%    intensity, or particle density, at the point [X, XP] with the array
%    INTENSE.  The 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).
%    The default value is 1.
%
%    emiterrors(x, xp, xerr, xperr, intense, interr) - specifies an error
%    on the intensity with the array INTERR.  The units are assumed to be
%    the same as INTENSE.  As with INTENSE, the INTERR array must either be
%    a single value - indicating an equal error on each value - or an array
%    of the same size as INTENSE.  The default value is 0 ie. no intensity
%    error.
%
%    emiterrors(x, xp, xerr, xperr, intense, interr, dispvars) - display
%    variances for each stage of summation; used for troubleshooting and
%    debugging.
%
%    [emierror] = emiterrors(...) - outputs the absolute rms emittance
%    error to the variable EMIERROR.
%
%    [emierror, emierrorfrac] = emiterrors(...) - also outputs the
%    fractional rms emittance error to the variable EMIERRORFRAC.

% 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
%       Complete calculation of emittance errors, based on DIPAC paper 1st
%       draft.  Tested with FETS pepperpot data.
%
%======================================================================

%    Parse input and output arguments

    if ( nargin < 4 ) || ( nargin > 7 )
        error('Incorrect number of input variables: syntax is emiterrors(x, xp, xerr, xperr, intense, interr, dispvars)') ;
    elseif nargout > 2
        error('Too many output variables: syntax is [emierror, emierrorfrac] = emiterrors(...)') ;
    end

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

    if ( length(xp) ~= numvals )
        error('EmitErrors: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 < 7 ) || ( isempty(dispvars) )
        dispvars = 0 ;
    end

    if ( nargin < 6 ) || ( isempty(interrin) )
        interr = zermat ;
    elseif ( length(interrin) == 1 )
        interr = onemat.*interrin ;
    else
        interr = reshape(interrin,[],1) ;
    end

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

    if ( nargin < 5 ) || ( 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

    if ( nargin == 3 )
        error('Must specify X and XP errors in pairs: emiterrors(x, xp, xerr, xperr, intense, interr)') ;
    end

    if ( isempty(xperrin) )
        xperr = zermat ;
    elseif ( length(xperrin) == 1 )
        xperr = onemat.*xperrin ;
    else
        xperr = reshape(xperrin,[],1) ;
    end

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

    if ( isempty(xerrin) )
        xerr = zermat ;
    elseif ( length(xerrin) == 1 )
        xerr = onemat.*xerrin ;
    else
        xerr = reshape(xerrin,[],1) ;
    end

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

%    Initialise variables for sum data

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

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

%    Generate sums and mean values

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

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

    emisquared = (xp2.*x2) - (xxp.^2) ;
    emittance = sqrt( emisquared ) ;
    if nargout == 0 || dispvars
        disp(' ') ;
        disp(['Un-normalised RMS emittance is: ' num2str(emittance,8) ' [pi mm mrad].']) ;
    end

    int2 = intense.^2 ;
    intsum2 = numpart.^2 ;
    intsum4 = numpart.^4 ;

%    Calculate error terms

    emitvar = zeros(6,1) ; emierror = 0 ; emierrorfrac = 0 ;

%    Error given by (insert into LaTeX): 
%    \begin{equation}
%    \begin{split}
%    &\sigma_{\varepsilon}^{2} = 
%     \frac{ \left( \sum_{i=1}^{N} \rho_{i}^{2} x_{i}^{2} \sigma_{x_{i}}^{2} 
%     + \dfrac{ x_{i}^{4} \sigma_{\rho_{i}}^{2}}{4} \right) 
%     \left( \sum_{j=1}^{N}\rho_{j} x_{j}'^{2} \right)^{2} }
%     {\varepsilon^{2} \left( \sum_{k=1}^{N}\rho_{k} \right)^{4}} \\
%     &+ \frac{ \left( \sum_{i=1}^{N} \rho_{i}^{2} x_{i}'^{2} \sigma_{x'_{i}}^{2} 
%     + \dfrac{ x_{i}'^{4} \sigma_{\rho_{i}}^{2}}{4} \right) 
%     \left( \sum_{j=1}^{N}\rho_{j} x_{j}^{2} \right)^{2} }
%     {\varepsilon^{2} \left( \sum_{k=1}^{N}\rho_{k} \right)^{4} } + \\
%     &\frac{\left( \sum_{i=1}^{N}\limits x_{i}^{2} x_{i}'^{2} \sigma_{\rho_{i}}^{2} 
%     + \rho_{i}^{2} x_{i}'^{2} \sigma_{x_{i}}^{2} + \rho_{i}^{2} x_{i}^{2} \sigma_{x_{i}'}^{2} \right)
%     \left( \sum_{j=1}^{N}\limits\rho_{j} x_{j} x_{j}' \right)^{2} }
%     {\varepsilon^{2} \left( \sum_{k=1}^{N}\rho_{k} \right)^{4} } \\
%     &- \frac{ 2 \left( \sum_{i=1}^{N} \rho_{i}^{2} x_{i}^{4} x_{i}'^{4} \sigma_{\rho_{i}}^{2} 
%     + \rho_{i}^{4} x_{i}^{2} x_{i}'^{4} \sigma_{x_{i}}^{2} 
%     + \rho_{i}^{4} x_{i}^{4} x_{i}'^{2} \sigma_{x_{i}'}^{2} \right) }
%     {\varepsilon^{2} \left( \sum_{k=1}^{N}\rho_{k} \right)^{4} } \\
%     &+ \frac{ \varepsilon^{2} \left( \sum_{i=1}^{N} \sigma_{\rho_{i}}^{2} \right) }
%     { \left( \sum_{k=1}^{N}\rho_{k} \right)^{2} }
%    \label{emitvar5}
%    \end{split}
%    \end{equation}

%    Calculate first term

    emitvar(1) = ( sum( ( int2.*x2arr.*(xerr.^2) ) + ( (x2arr.^2).*(interr.^2)./4 ) ) ).*...
                    (xp2Sum.^2)./( emisquared.*intsum4 ) ;
    if dispvars, disp(['Emittance variance term 1 = ' num2str(emitvar(1),8)]), end ;

%    Calculate second term

    emitvar(2) = ( sum( ( int2.*xp2arr.*(xperr.^2) ) + ( (xp2arr.^2).*(interr.^2)./4 ) ) ).*...
                    (x2Sum.^2)./( emisquared.*intsum4 ) ;
    if dispvars, disp(['Emittance variance term 2 = ' num2str(emitvar(2),8)]), end ;

%    Calculate third term

    emitvar(3) = ( sum( ( x2arr.*xp2arr.*(interr.^2) ) + ( int2.*xp2arr.*(xerr.^2) ) + ...
                    ( int2.*x2arr.*(xperr.^2) ) ) ).*(xxpSum.^2)./( emisquared.*intsum4 ) ;
    if dispvars, disp(['Emittance variance term 3 = ' num2str(emitvar(3),8)]), end ;

%    Calculate fourth term

    emitvar(4) = -2.*( sum( ( int2.*(x2arr.^2).*(xp2arr.^2).*(interr.^2) ) + ...
                    ( (int2.^2).*x2arr.*(xp2arr.^2).*(xerr.^2) ) + ...
                    ( (int2.^2).*(x2arr.^2).*xp2arr.*(xperr.^2) ) ) )./( emisquared.*intsum4 ) ;
    if dispvars, disp(['Emittance variance term 4 = ' num2str(emitvar(4),8)]), end ;

%    Calculate fifth term

    emitvar(5) = ( emisquared.*( sum( interr.^2 ) ) )./intsum2 ;
    if dispvars, disp(['Emittance variance term 5 = ' num2str(emitvar(5),8)]), end ;

    emitvar(6) = sum(emitvar(1:5)) ;
    if nargout == 0 || dispvars
        disp(' ')
        disp(['Emittance variance = ' num2str(emitvar(6),8)])
    end

    emierror = sqrt(emitvar(6)) ;
    if nargout == 0 || dispvars
        disp(' ')
        disp(['Un-normalised Emittance error = ' num2str(emierror,8) ' [pi mm mrad].'])
    end

    emierrorfrac = emierror./emittance ;
    if nargout == 0 || dispvars
        disp(' ') ; disp(['Fractional Emittance error = ' num2str(emierrorfrac,8)]) ;
        disp(' ')
    end

    return

