function [ppposdat,ppemitdat,pppixemit5d,pppixemit2d,ppemithist] = ...
    ppanalyse(ppbgfname,lpos,rpos,midpos,ppimfname,pplim,width,zlen,rotang,energy,plt)
% Function PPANALYSE.M
%
%    [ppposdat,ppemitdat,pppixemit5d,pppixemit2d,ppemithist] =
%       ppanalyse(ppbgfname,lpos,rpos,midpos,ppimfname,pplim,width,zlen,rotang,energy,plt)
%
%    The function PPANALYSE rather ham-fistedly collects together a series
%    of subfunctions used to calculate emittance data from pepperpot
%    images.  The required input variables are as follows:
%
%        PPBGFNAME - the filename of the background image.
%
%        LPOS, RPOS - the left- and right-hand ruler positions for the
%                     background image.
%
%        MIDPOS - the central position calculated from the background
%                 image.
%
%        PPIMFNAME - the filename of the data image.
%
%        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(1:2) is passed to
%               PEPPERPLOTFUN for the background and data images; PLT(3)
%               is passed to PPTILTFIND; PLT(4) = 1 plots the hot data
%               image; PLT(5) = 1 plots the reconstructed spot image;
%               PLT(6) = 1 plots the 3D intensity data; PLT(7) = 1 plots
%               the emittance plots; PLT(8) = 1 plots the intensity profile
%               of the central spot.
%
%    The output variables are as follows:
%
%        PPPOSDAT - position data for pepperpot data image, including
%                   intensity information.
%
%        PPEMITDAT - emittance information, straight from PPEMITCALC.
%
%        PPPIXEMIT5D - 5-d matrix containing emittance data for ever 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', ...
            ['11 input variables are required: [ppposdat,ppemitdat,pppixemit5d,pppixemit2d] = '...
            '\n ppanalyse(ppbgfname,lpos,rpos,midpos,ppimfname,pplim,width,zlen,rotang,energy,plt)']) ;
    end

    if ( ~ischar(ppbgfname) ) | ( ~ischar(ppimfname) )
        error('Background and main image filenames must be strings') ;
    end

    if length(plt) ~= 8
        error('PLT variable must contain 8 parameters (see the help)') ;
    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']) ;

%    Load background image data

    if (ppbgfname(end-2:end) == 'mat')
        load(ppbgfname,'ppimbg','stdmbg','pphistbg','ppbgminmax') ;
    else
        [ppimbg,stdmbg,pphistbg,ppbgminmax] = pepperplotfun(ppbgfname,0,plt(1:2)) ;
    end

    ppbgmin = ppbgminmax(1) ; ppbgmax = ppbgminmax(2) ;

%    Load calibration data

    cur = cd ;
    cd([matlabroot '/work/pepperpot/cal_data']) ;
    load ppmk1_caldat.mat ;
    cd([cur]) ;
    vsp = size(ppcalmat,1) ;
    hsp = size(ppcalmat,2) ;

%    Find image tilt and reconstruct hole positions

    [midp,imtilt,rulen] = pptiltfind(ppimbg,lpos,rpos,plt(3)) ;

    [ppimspots] = ppspotreco(ppcalmat,midpos,rulen,imtilt) ;
    [ppimspotsb] = ppspotreco(ppcalmat,midpos,rulen,imtilt,1.02) ;

%    Load main image data and plot hot image

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

    if plt(4) == 1
        figure, imagesc(ppim,[90 250]), axis image, hold on ; colormap(hot) ;
    end

%    Find and plot image spots

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

    if plt(5) == 1

        figure, imagesc(ppim,[ppmin 150]), 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

%        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

    rupix = rulen./80 ;
    sppix = (rulen.*mlrulratio)./72 ;
    phospix = sppix.*perscale ;
    invintense = squeeze(1./ppcalmat(:,:,3)) ;

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

    ppimspotsoff(:,:,1) = ppimspots(:,:,1) - ppimspots(13,13,1) ;
    ppimspotsoff(:,:,2) = ppimspots(:,:,2) - ppimspots(13,13,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 ;

    if plt(6) == 1

%        figure, bar3(ppimdata(:,:,11)) ;
        figure, surfc(flipud(ppimdata(:,:,11))) ;% colormap(jet) ;
        axis tight ; axis square ;
%        figure, surf(ppimdata(:,:,6),ppimdata(:,:,7),ppimdata(:,:,11)) ;

    end

%    Calculate emittance data from [r,theta,phi]

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

    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 phase space distribution, x-x'', for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        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 phase space distribution, y-y'', for ' num2str(energy./1e3) 'keV Pepperpot image']) ;
        hold off

    end

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

    eminames = {['emioutx']; ['betax']; ['gammax']; ['xrms']; ['xprms']; ['emiouty']; ...
        ['betay']; ['gammay']; ['yrms']; ['yprms']} ;

    ppemitdat = ([emioutx, emiouty, xrms, yrms, xprms, yprms, betax, betay, gammax, gammay]) ;

%    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

    ppimstdxyrtp = pprtpconv(ppimstdxyp,zlen,phospix) ;
    ppimstdxy = ppimstdxyrtp(:,:,5).*1000 ;

%    Plotting test for pixel emittance data - plotting central spot

    if 0

        spotxydat = ppspotallxy(ppim,ppimspots(cspv,csph,1:2),ppimanal(cspv,csph,1:2),width,pplim) ;
        spotrtpdat = pprtpconv(spotxydat(:,:,1:2),zlen,phospix) ;
        holertpdat = zeros(size(spotrtpdat)) ;
        holertpdat(:,:,3) = ppimdata(cspv,csph,3) ;
        holertpdat(:,:,4) = ppimdata(cspv,csph,4) ;
        holertpdat(:,:,5) = ppimdata(cspv,csph,5) ;

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

        minijet = colormap(jet(200)) ;
        minijet(1,:) = 1 ;

        figure, hold on ;
        surfc(spotrtpdat(:,:,1),spotrtpdat(:,:,2),spotxydat(:,:,3)) ; %,'EdgeColor','interp') ;
        colormap(minijet) ;
        axis tight ; axis square ;
        hold off ;

        figure, hold on ;
        surfc(spotalldat(:,:,3),spotalldat(:,:,4),spotalldat(:,:,5)) ; %,'EdgeColor','interp') ;
        colormap(minijet) ;
        axis tight ; axis square ;
        hold off ;

    end

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

    if nargout > 2

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

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

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

                spotxydat = ppspotallxy(ppim,ppimspots(i,j,1:2),ppimanal(i,j,1:2),width,pplim) ;
                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,:,:,:) = spotalldat(:,:,:) ;

            end
        end

%    Reshape pixel emittance data array to N x 5, containing only spots
%    above intensity threshold

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

    end

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

    if nargout > 4

        intens = round(pppixemit2d(:,5)./10) ;
        [ppemithist] = ppemithistgenxy(pppixemit2d(:,1),pppixemit2d(:,2),pppixemit2d(:,3), ...
            pppixemit2d(:,4),intens,3,3) ;

    end

    return


