function [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d,ppemithist] = ...
    ppanalyse41(ppimfname,ppbgfname,cent,imtilt,sp2sp,pplim,width,zlen,rotang,energy,plt,bitdiv,beamview,histpart)
% Function PPANALYSE41.M
%
%    [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d,ppemithist] =
%       ppanalyse41(ppimfname,ppbgfname,cent,imtilt,sp2sp,pplim,width,zlen,rotang,energy,plt,bitdiv,beamview,histpart)
%
%    The function PPANALYSE41 rather ham-fistedly collects together a
%    series of subfunctions used to calculate emittance data from pepperpot 
%    images using the 41x41 mesh screen.  All input variables except BITDIV
%    and HISTPART are required to make PPANALYSE41 run.
%
%    [ppposdat,ppemitdat,ppim] = ppanalyse41(ppimfname,ppbgfname,cent, ...
%           imtilt,sp2sp,pplim,width,zlen,rotang,energy,plt)
%       - calculate spot-weighted emittance data.
%
%    [ppposdat,ppemitdat,ppim] = ppanalyse41(...,bitdiv) - additionally
%       specify divider for grayscale values.
%
%    [ppposdat,ppemitdat,ppim] = ppanalyse41(...,bitdiv,beamview) -
%       additionally specify orientation of pepperpot image.
%
%    [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d] = ppanalyse41(...)
%       - calculate individual pixel emittance data.
%
%    [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d,ppemithist] = ...
%           ppanalyse41(...,[],histpart)
%       - calculate individual pixel emittance data and reproduce
%       individual particle data, for HISTPART particles.
%
%    It is possible to specify any number of output variables, between 1
%    and 6, but this changes the behaviour of PPANALYSE41: specifying 3 or
%    fewer output variables only analyses the spot-weighted emittance;
%    specifying more than 3 output variables also runs the single pixel
%    emittance analysis.  The recommended configurations are as follows:
%
%    [ppposdat,ppemitdat,ppim] = ppanalyse41(...) - outputs analysed image
%    data (PPPOSDAT) and emittance data (PPEMITDAT) for the spot-weighted
%    image data, plus the data image itself (PPIM) with any background
%    image subtraction.
%
%    [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d] = ppanalyse41(...)
%    - in addition to the above, also calculates and outputs the single
%    pixel emittance data as a 5-D matrix (PPPIXEMIT5D) containing all the
%    data, and a 2-D matrix (PPPIXEMIT2D) with just the relevant pixel data
%    for plotting.
%
%    [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d,ppemithist] = ppanalyse41(...) 
%    - creates X, Y, X' and Y' data for individual particles, based on the
%    intensity data from PPPIXEMIT2D, and outputs an [N x 4] matrix
%    (PPEMITHIST) of N particles.  Emittance plots are also generated: this
%    is the best option for producing emittance plots for presentations.
%
%    The required input variables are as follows:
%
%        PPIMFNAME - the filename of the data image.
%
%        PPBGFNAME - the filename of the background image for background
%                    subtraction - specify as [] if no background
%                    subtraction required.
%
%        CENT - the centre of the image.
%
%        IMTILT - the angle of tilt of the screen with respect to the
%                 camera.
%
%        SP2SP - the spot-to-spot spacing for the screen, in pixels.
%
%        PPLIM - the cut to apply to the data
%
%        WIDTH - the width of the search window used to calculate the mean
%                weighted spot position.
%
%        ZLEN - the longitudinal distance from tungsten screen to phosphor
%               screen.
%
%        ROTANG - the rotation angle of the pepperpot assembly.
%
%        ENERGY - the energy of the beam in eV.
%
%        PLT - a 1D array of plotting data: PLT = [2 2 1 1 1 3 1 1 1 2 3]
%           PLT(1:2) is passed to PEPPERPLOTFUN for the background and data images; 
%           PLT(3) = 1 plots the hot data image;
%           PLT(4) = 1 plots the reconstructed spot image; 
%           PLT(5) = 1 plots the parts of the image above the cut level; 
%           PLT(6) = 1 plots the spot intensity data in 3D,
%                  = 2 plots the intensity data as a contour plot,
%                  = 3 plots both; 
%           PLT(7) = 1 plots the emittance plots;
%           PLT(8) = 1 plots the intensity profile of the central spot;
%           PLT(9) = 1 plots the quiver plots; 
%           PLT(10) = 1 plots the individual pixel emittance data;
%                   = 2 plots it in colour;
%           PLT(11) = 1 plots the rehistogrammed particle data (ppemithist);
%                   = 2 plots the corresponding emittance plots;
%                   = 3 plots them both.
%        It is not necessary to specify PLT(9:11).
%
%    Optional input variables:
%
%        BITDIV - factor to divide grayscale bit values by, to account for
%                 PCO 14-bit images.  Does not have to be specified; the
%                 default is 4.
%
%        BEAMVIEW - selects the orientation of the data image loaded by
%                   PEPPERPLOTFUN.  Emittance data should be calculated using
%                   beam's-eye-view data; camera's-eye-view data is easier
%                   for viewing.  Setting BEAMVIEW = ('CAM') loads and displays
%                   the image as it was viewed by the camera; BEAMVIEW = ('BEAM')
%                   flips the image left-to-right, to show the
%                   beam's-eye-view image.  BEAMVIEW = ('CAM') is the default
%                   value.
%
%        HISTPART - specifies the number of particles to produce for the
%                   matrix PPEMIHIST using the rehistogramming method. Does
%                   not have to be specified; the default is 10,000
%                   particles.
%
%    The output variables are as follows:
%
%        PPPOSDAT - position data for pepperpot data image, including
%                   intensity information, mostly from PPEMITMAT.
%
%        PPEMITDAT - a [10 x 3] structure array with the following 3
%                    fields:
%                    .name - the type of emittance data ie. 'xemit',
%                            'yprms' etc.
%                    .avespot - corresponding data for average spot
%                               calculation
%                    .singlepix - corresponding data for single pixel
%                                 calculation
%                    Emittance data is taken straight from PPEMITMAT.
%
%        PPIM - pepperpot image data, straight from PEPPERPLOTFUN, then
%               with the background image subtracted (if specified).
%
%        PPPIXEMIT5D - 5-d matrix containing emittance data for every
%                      single pixel with a spot (as defined by WIDTH).
%
%        PPPIXEMIT2D - reshaped version of PPPIXEMIT5D, containing only
%                      those pixels above threshold (PPLIM) in an N x 5
%                      matrix.
%
%        PPEMITHIST - histogram data generated using the data from
%                     PPPIXEMIT2D and randomly assigning X and Y data using
%                     the 3mm hole spacing.
%

    if nargin < 11
        error('PPAnalyse:incorrectSyntax', ...
            ['At least 11 input variables are required: [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d,ppemithist] = '...
            '\n ppanalyse41(ppimfname,ppbgfname,cent,imtilt,sp2sp,pplim,width,zlen,rotang,energy,plt)']) ;
    end

    if nargin > 14
        error('PPAnalyse:incorrectSyntax', ...
            ['Can only specify 14 input variables: [ppposdat,ppemitdat,ppim,pppixemit5d,pppixemit2d,ppemithist] = '...
            '\n ppanalyse41(ppimfname,ppbgfname,cent,imtilt,sp2sp,pplim,width,zlen,rotang,energy,plt,bitdiv,beamview,histpart)']) ;
    end

    if ( ~ischar(ppimfname) )
        error('Main image filename must be a strings') ;
    end

    if ( ~isempty(ppbgfname) ) && ( ~ischar(ppbgfname) )
        error('Background image filename must be a strings') ;
    end

    if length(plt) < 8
        error('PLT variable must contain at least 8 parameters (see the help)') ;
    end

    if length(plt) < 9
        plt(9) = 0 ;
    end

    if length(plt) < 10
        plt(10) = 1 ;
    end

    if length(plt) < 11
        plt(11) = 2 ;
    end

    if ( nargin < 14 ) || ( isempty(histpart) )
        histpart = 1e4 ;
    end

    if ( nargin < 13 ) || ( isempty(beamview) )
        beamview = ('cam') ;
    end

    if ( nargin < 12 ) || ( isempty(bitdiv) )
        bitdiv = 4 ;
    end

    if ( isempty(width) )
        width = ( round( sp2sp./4 ) ).*4 - 1 ;
    end

    plotcol = (['r+'; 'g+'; 'b+'; 'c+'; 'm+'; 'y+'; 'k+'; ...
        'r*'; 'g*'; 'b*'; 'c*'; 'm*'; 'y*'; 'k*'; ...
        'rx'; 'gx'; 'bx'; 'cx'; 'mx'; 'yx'; 'kx'; ...
        'ro'; 'go'; 'bo'; 'co'; 'mo'; 'yo'; 'ko'; ...
        'rs'; 'gs'; 'bs'; 'cs'; 'ms'; 'ys'; 'ks'; ...
        'rd'; 'gd'; 'bd'; 'cd'; 'md'; 'yd'; 'kd'; ...
        'r^'; 'g^'; 'b^'; 'c^'; 'm^'; 'y^'; 'k^'; ...
        'rv'; 'gv'; 'bv'; 'cv'; 'mv'; 'yv'; 'kv'; ...
        'r>'; 'g>'; 'b>'; 'c>'; 'm>'; 'y>'; 'k>'; ...
        'r<'; 'g<'; 'b<'; 'c<'; 'm<'; 'y<'; 'k<'; ...
        'rp'; 'gp'; 'bp'; 'cp'; 'mp'; 'yp'; 'kp'; ...
        'rh'; 'gh'; 'bh'; 'ch'; 'mh'; 'yh'; 'kh']) ;

%    Load calibration data

%    cur = cd ;
%    cd([matlabroot '/work/pepperpot/cal_data']) ;
    load ppmk2_caldat.mat ;
%    cd([cur]) ;

%    Load main image data and plot hot image

    [ppim,stdm,pphist,ppminmax] = pepperplotfun(ppimfname,0,plt(1:2),bitdiv,beamview) ;
    ppmin = ppminmax(1) ; ppmax = ppminmax(2) ;

    if pplim < ppmin
        pplow = ppmin ;
    else
        pplow = pplim ;
    end

    pphi = pplow + 200 ;

%    if ppmax > 2.*pplim
%        pphi = ppmax ;
%    else
%        pphi = 2.*pplim ;
%    end

    if plt(3) == 1
        figure, imagesc(ppim,[pplow pphi]), axis image, hold on ; colormap(hot) ;
        hold off ;
    end

%    Load background image data

    if ( ~isempty(ppbgfname) )
        [ppimbg,stdmbg,pphistbg,ppbgminmax] = pepperplotfun(ppbgfname,0,plt(1:2),bitdiv) ;
        ppbgmin = ppbgminmax(1) ; ppbgmax = ppbgminmax(2) ;
        ppimc = ppim - ppimbg ;
        ppbgcut = ( ppimc > 0 ) ;
        ppim = ppimc.*ppbgcut ;

        figure, imagesc(ppim,[0 max(max(ppim))]), axis image ;
        title(['\bf Scaled Image Data with background subtracted']) ;
    end

%    Refind central hole position to centre grid

    centwidth = ( round( sp2sp./4 ) ).*4 - 1 ;

    [cent1,stdxy,spint] = ppspot(ppim,cent,centwidth,pplim) ;
    [cent2,stdxy,spint] = ppspot(ppim,cent1,centwidth,pplim) ;
    [cent3,stdxy,spint] = ppspot(ppim,cent2,centwidth,pplim) ;
    [cent4,stdxy,spint] = ppspot(ppim,cent3,centwidth,pplim) ;

    if spint > 0
        cent = round( cent4 ) ;
    end

    disp(' ') ;
    disp(['New centre position found at: [' num2str(cent) '] ']) ;
    disp(' ') ;

%    Reconstruct hole positions

%    [ppimspots] = ppspotreco41(ppcalmat,cent,sp2sp,imtilt) ;
%    [ppimspotsb] = ppspotreco41(ppcalmat,cent,sp2sp,imtilt,1.02) ;
    [ppimspots] = ppmakegrid41(cent,sp2sp,imtilt) ;
%    [ppimspotsb] = ppmakegrid41(cent,sp2sp,imtilt) ;
    [ppimspotsb] = ppmakegrid41(cent,(sp2sp.*1.02),imtilt) ;

    vsp = size(ppimspots,1) ;
    hsp = size(ppimspots,2) ;

%    Find and plot image spots

    [ppimanal,ppimstdxyp] = ppallspots(ppim,ppimspots(:,:,1:2),ppimspotsb(:,:,1:2),width,pplim) ;

    if plt(4) == 1

        figure, imagesc(ppim,[pplow pphi]), axis image, hold on ;
        colormap(gray) ;
        for i = 1:vsp
            for j = 1:hsp
                plot(ppimspots(i,j,1),ppimspots(i,j,2),'xc') ;
                if ppimanal(i,j,5) > 0
                    plot(ppimanal(i,j,1),ppimanal(i,j,2),'xr') ;
                end
            end
        end
        hold off

    end

    if plt(5) == 1

        ppimc = ppim > pplim ;
        figure, imagesc(ppimc,[0 1]), axis image, hold on ;
        colormap(gray) ;
%        for i = 1:vsp
%            for j = 1:hsp
%                plot(ppimspots(i,j,1),ppimspots(i,j,2),'xr') ;
%                plot(ppimanal(i,j,1),ppimanal(i,j,2),'xc') ;
%            end
%        end
        hold off

    end

    phospix = sp2sp./3 ;
    invintense = squeeze(1./ppcalmat(:,:,3)) ;

%    Calculate and display resolution

    posres = 1./phospix ;
    angres = ( tan(1./(phospix.*zlen)) ).*1e3 ;

    disp(' ') ;
    disp(['Camera resolution: ' num2str(posres.*1e3) ' microns per pixel']) ;
    disp(['Angular resolution: ' num2str(angres) ' mrad']) ;
    disp(' ') ;

%    Main hole positions in (X,Y) and (r,theta,phi)

    ppimspotsoff(:,:,1) = ppimspots(:,:,1) - ppimspots(21,21,1) ;
    ppimspotsoff(:,:,2) = ppimspots(:,:,2) - ppimspots(21,21,2) ;
    ppimdata = pprtpconv(ppimspotsoff,zlen,phospix) ;

%    Spot positions in (X,Y) and (r,theta,phi)

    ppimdata(:,:,6:10) = pprtpconv(ppimanal(:,:,3:4),zlen,phospix) ;

%    Intensity data

    ppimdata(:,:,11) = ppimanal(:,:,5) ; %.*invintense ;

%    Calculate emittance data from [r,theta,phi] and plot beam intensities

    [ppemitdata,ppposdat] = ppemitmat(ppimdata(:,:,3:5), ppimdata(:,:,8:10), rotang) ;
    ppposdat(:,:,end+1) = ppimdata(:,:,11) ;

%    ppposdat(:,:,[end+1 end+2]) = ppimstdxyp ;

    intmat = ppposdat(:,:,13) ;
    intzer = (intmat > 0 ) ;
    intmax = max(max(intmat)) ;
    intmat = ( ( intmat + ( intmax./4.5)  ).*intzer ) ;

    if ( plt(6) == 1 ) || ( plt(6) == 3 ) 

        figure, surfc(ppposdat(:,:,1),ppposdat(:,:,2),ppposdat(:,:,13)) ;
        axis tight ; axis square ; hold on ; colormap(hot) ;
        xlabel(['Horizontal position, x (mm)']), ylabel(['Vertical position, y (mm)']), ...
            zlabel(['Integrated spot intensity']), ...
            title(['\bf Beam profile for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        hold off ;

    end

    if plt(6) > 1

%        figure, contourf(ppposdat(:,:,1),ppposdat(:,:,2),ppposdat(:,:,13)) ;
        figure, contourf(ppposdat(:,:,1),ppposdat(:,:,2),intmat) ;
        axis tight ; axis square ; hold on ; colormap(hot) ;
        set(gca,'Color',[0 0 0]) ;
        xlabel(['Horizontal position, x (mm)']), ylabel(['Vertical position, y (mm)']), ...
            zlabel(['Integrated spot intensity']), ...
            title(['\bf Beam profile for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        hold off ;

    end

%    Calculate emittances

    [emioutx, betax, gammax, xrms, xprms] = ppemitcalc(ppemitdata(:,:,1), ppemitdata(:,:,3), ppimdata(:,:,end), energy, [], ['X']) ;
    [emiouty, betay, gammay, yrms, yprms] = ppemitcalc(ppemitdata(:,:,2), ppemitdata(:,:,4), ppimdata(:,:,end), energy, [], ['Y']) ;

    eminames = {['xemit'], ['yemit'], ['xrms'], ['yrms'], ['xprms'], ['yprms'], ...
        ['betax'], ['betay'], ['gammax'], ['gammay']} ;

    shortemit = {emioutx, emiouty, xrms, yrms, xprms, yprms, betax, betay, gammax, gammay} ;

    ppemitdat = struct('name', eminames, 'avespot', shortemit) ;

%    Plot emittance data

    if plt(7) == 1

        figure, hold on ;
        for i = 1:vsp
            for j = 1:hsp
                if ppimdata(i,j,11) ~= 0
%                    plot(ppemitdata(i,j,1),ppemitdata(i,j,3),'xk') ;
                    plot(ppemitdata(i,j,1),ppemitdata(i,j,3),plotcol(i,:)) ;
                end
            end
        end
        xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
            title(['\bf Horizontal emittance, \epsilon_x, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        text(0.05,0.9,['\epsilon_x = ' num2str(emioutx,'%0.4g') '\pi mm mrad'],'Units','normalized') ;
        hold off ;

        figure, hold on ;
        for i = 1:vsp
            for j = 1:hsp
                if ppimdata(i,j,11) ~= 0
%                    plot(ppemitdata(i,j,2),ppemitdata(i,j,4),'xk') ;
                    plot(ppemitdata(i,j,2),ppemitdata(i,j,4),plotcol(j,:)) ;
                end
            end
        end
        xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
            title(['\bf Vertical emittance, \epsilon_y, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        text(0.05,0.9,['\epsilon_y = ' num2str(emiouty,'%0.4g') '\pi mm mrad'],'Units','normalized') ;
        hold off ;

    end

%    Count good rows and columns and calculate angular resolution in mrad

    goodspots = ( ppposdat(:,:,13) > 0 ) ;
    goodxdat = sum( sum(goodspots) > 0 ) ;
    goodydat = sum( sum(goodspots') > 0 ) ;

%    Plot intensity profile for central spot

    csph = round(hsp./2) ;
    cspv = round(vsp./2) ;

    if plt(8) == 1
        pp3dspot(ppim,squeeze(ppimanal(cspv,csph,1:2)),width,1,pplim,squeeze(ppimanal(cspv,csph,1:2))) ;
        title(['\bf Central spot intensity distribution for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        view(-45,25) ;
    end

%    Plot quiver plot of intensity and emittance

    if plt(9) == 1

        numholes = numel(ppposdat(:,:,1)) ;
        quivmat = zeros(numholes,5) ;
        quivmat(:,1) = reshape(ppposdat(:,:,1),numholes,1) ;
        quivmat(:,2) = reshape(ppposdat(:,:,2),numholes,1) ;
        quivmat(:,3) = reshape(ppposdat(:,:,6),numholes,1) ;
        quivmat(:,4) = reshape(ppposdat(:,:,7),numholes,1) ;
        quivmat(:,5) = reshape(ppposdat(:,:,13),numholes,1) ;
        numpart = sum(quivmat(:,5)) ;

        inds = find(quivmat(:,5) <= 0) ;
        quivmat(inds,:) = [] ;

        xpMean = sum(quivmat(:,3).*quivmat(:,5))./numpart ;
        ypMean = sum(quivmat(:,4).*quivmat(:,5))./numpart ;

        x = quivmat(:,1) ;
        y = quivmat(:,2) ;
        u = quivmat(:,3) - xpMean ;
        v = quivmat(:,4) - ypMean ;

        hjet = colormap(jet) ;
        hjet = hjet + 9 ;
        hjet = hjet./10 ;
        hjet(1,:) = [1 1 1] ;

%        figure, contourf(ppposdat(:,:,1),ppposdat(:,:,2),ppposdat(:,:,end)) ;
        figure, contourf(ppposdat(:,:,1),ppposdat(:,:,2),intmat) ;
        axis tight ; axis square ; hold on ; colormap(hjet) ;
        quiver(x,y,u,v,'-k','filled') ;
        xlabel(['Horizontal position, x (mm)']), ylabel(['Vertical position, y (mm)']), ...
            title(['\bf Beam profile and emittance for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        hold off ;

    end

    ppimstdxyrtp = pprtpconv(ppimstdxyp,0,phospix) ;
%    ppimstdxy = ppimstdxyrtp(:,:,5).*1000 ;
%    ppposdat(:,:,end+1) = ppimstdxyrtp(:,:,3) ;

%    Find x' and y' data for each pixel with each spot

    if nargout > 3

        disp('Calculating emittance data for all pixels... ') ;
        disp(' ') ;

        wbh = waitbar(0,'Calculating emittance data for all pixels... ','Name','Single Pixel Progress') ;

        pppixemit5d = zeros(vsp,hsp,width,width,7) ;

        cutmat = ( ppim > pplim ) ;
        immat = ( ppim - pplim ).*cutmat ;

        for i = 1:vsp
            for j = 1:hsp

%                spotxydat = ppspotallxy(ppim,ppimspots(i,j,1:2),ppimanal(i,j,1:2),width,pplim) ;
                spotxydat = ppspotallxy(immat,ppimspots(i,j,1:2),ppimanal(i,j,1:2),width) ;

                spotrtpdat = pprtpconv(spotxydat(:,:,1:2),zlen,phospix) ;

                holertpdat = zeros(size(spotrtpdat)) ;
                holertpdat(:,:,3) = ppimdata(i,j,3) ;
                holertpdat(:,:,4) = ppimdata(i,j,4) ;
                holertpdat(:,:,5) = ppimdata(i,j,5) ;

                [spotalldat] = ppemitmat(holertpdat(:,:,3:5), spotrtpdat(:,:,3:5), rotang) ;
                spotalldat(:,:,end+1) = spotxydat(:,:,end) ;

                pppixemit5d(i,j,:,:,1:5) = spotalldat(:,:,:) ;
                pppixemit5d(i,j,:,:,6) = j ;
                pppixemit5d(i,j,:,:,7) = i ;

                spotnum = ( (i - 1).*hsp ) + j ;
                waitbar((spotnum./(hsp.*vsp)),wbh) ;

            end
        end

        close(wbh) ;

%    Reshape pixel emittance data array to N x 7, containing only pixels
%    above intensity threshold.
%    The columns are ordered as follows:
%        1: X position of hole (mm) - X
%        2: Y position of hole (mm) - Y
%        3: X' angle from hole to pixel (mrad) - X'
%        4: Y' angle from hole to pixel (mrad) - Y'
%        5: Pixel intensity (counts) - I
%        6: Horizontal grid position (column number) - for plotting
%        7: Vertical grid position (row number) - for plotting

        disp('Reshaping emittance data and removing empty pixels... ') ;
        disp(' ') ;

        leng2d = numel(pppixemit5d)./7 ;
        pppixemit2d = reshape(pppixemit5d,leng2d,7) ;
        inds = find(pppixemit2d(:,5) <= 0) ;
        pppixemit2d(inds,:) = [] ;

%    Recalculate emittance values

        disp('Recalculated emittance values ... ') ;

        [emioutx, betax, gammax, xrms, xprms] = ppemitcalc(pppixemit2d(:,1), pppixemit2d(:,3), pppixemit2d(:,5), energy, [], ['X']) ;
        [emiouty, betay, gammay, yrms, yprms] = ppemitcalc(pppixemit2d(:,2), pppixemit2d(:,4), pppixemit2d(:,5), energy, [], ['Y']) ;

        disp(' ') ;

        fullemit = {emioutx, emiouty, xrms, yrms, xprms, yprms, betax, betay, gammax, gammay} ;

        ppemitdat = struct('name', eminames, 'avespot', shortemit, 'singlepix', fullemit) ;

%    Plot reshaped emittance data

        if plt(10) == 1

            figure, hold on ;
            plot(pppixemit2d(:,1),pppixemit2d(:,3),'xk') ;
            xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
                title(['\bf Horizontal emittance, \epsilon_x, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
            text(0.05,0.9,['\epsilon_x = ' num2str(emioutx,'%0.4g') '\pi mm mrad'],'Units','normalized') ;
            hold off ;

            figure, hold on ;
            plot(pppixemit2d(:,2),pppixemit2d(:,4),'xk') ;
            xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
                title(['\bf Vertical emittance, \epsilon_y, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
            text(0.05,0.9,['\epsilon_y = ' num2str(emiouty,'%0.4g') '\pi mm mrad'],'Units','normalized') ;
            hold off ;

        end

        if plt(10) == 2

            figure, hold on ;
            for i = 1:length(pppixemit2d)
                plot(pppixemit2d(i,1),pppixemit2d(i,3),plotcol(pppixemit2d(i,7),:)) ;
            end
            xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
                title(['\bf Horizontal emittance, \epsilon_x, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
            text(0.05,0.9,['\epsilon_x = ' num2str(emioutx,'%0.4g') '\pi mm mrad'],'Units','normalized') ;
            hold off ;

            figure, hold on ;
            for i = 1:length(pppixemit2d)
                plot(pppixemit2d(i,2),pppixemit2d(i,4),plotcol(pppixemit2d(i,6),:)) ;
            end
            xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
                title(['\bf Vertical emittance, \epsilon_y, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
            text(0.05,0.9,['\epsilon_y = ' num2str(emiouty,'%0.4g') '\pi mm mrad'],'Units','normalized') ;
            hold off ;

        end

%    Plot reshaped emittance data as histogram

        if plt(10) > 0

            figure, hold on ;
            ppemithistplot(pppixemit2d(:,1),pppixemit2d(:,3),pppixemit2d(:,5),[],16) ;
            xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
                title(['\bf Horizontal emittance, \epsilon_x, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
            hold off ;

            figure, hold on ;
            ppemithistplot(pppixemit2d(:,2),pppixemit2d(:,4),pppixemit2d(:,5),[],16) ;
            xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
                title(['\bf Vertical emittance, \epsilon_y, for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
            hold off ;

        end

    end

%    Create histogram of emittance data by generating X and Y information

    if nargout > 5

        if ( nargin < 13 ) || ( isempty(histpart) )
            histpart = 1e4 ;
        end

        disp('Creating randomly generated data for histogramming... ') ;
        disp(' ') ;

        intsum = sum(pppixemit2d(:,5)) ;
        intdiv = intsum./histpart ;

        intens = (pppixemit2d(:,5)./intdiv) ;
        [ppemithist] = ppemithistgen(pppixemit2d(:,1),pppixemit2d(:,2),pppixemit2d(:,3), ...
            pppixemit2d(:,4),intens,3,3,angres,angres) ;

%    Plot regenerated emittance data

        if ( plt(11) == 1 ) || ( plt(11) == 3 ) 

            figure, hold on ;
            plot(ppemithist(:,1),ppemithist(:,3),'xk') ;
            xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
                title(['\bf Horizontal emittance, \epsilon_x, for regenerated ' ...
                num2str(energy./1e3) 'keV Pepperpot image']) ;
            hold off ;

            figure, hold on ;
            plot(ppemithist(:,2),ppemithist(:,4),'xk') ;
            xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
                title(['\bf Vertical emittance, \epsilon_y, for regenerated ' ...
                num2str(energy./1e3) 'keV Pepperpot image']) ;
            hold off ;

        end

%    Plot reshaped emittance data as histogram

        if ( plt(11) > 1 )

            figure, hold on ;
            ppemithistplot(ppemithist(:,1),ppemithist(:,3)) ;
            xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
                title(['\bf Horizontal emittance, \epsilon_x, for regenerated ' ...
                num2str(energy./1e3) 'keV Pepperpot image']) ;
            hold off ;

            figure, hold on ;
            ppemithistplot(ppemithist(:,2),ppemithist(:,4)) ;
            xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
                title(['\bf Vertical emittance, \epsilon_y, for regenerated ' ...
                num2str(energy./1e3) 'keV Pepperpot image']) ;
            hold off ;

        end

    end

    return


