% PPCALIB
%
%    PPCALIB - Calibrate pepperpot images

    clear all ; close all ;

    [ppimbg,stdmbg,pphistbg,ppbgminmax] = pepperplotfun('E:\pepperpot\04_05_06\500ms_bglight_10.tif',0,[2 0]) ;

    ppbgmin = ppbgminmax(1) ; ppbgmax = ppbgminmax(2) ;
    load E:\pepperpot\cal_data\pixcal.mat ;
    load E:\pepperpot\cal_data\ppcaldat_good.mat ;

    lpos = [207 1795] ;
    rpos = [1945 1900] ;

%    [midpos,imtilt,rulen] = pptiltfind(ppim) ;
    [midpos,imtilt,rulen] = pptiltfind(ppimbg,lpos,rpos) ;
    midp = [1115 980] ;

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

    figure, imagesc(ppimbg,[ppbgmin ppbgmax]), axis image, hold on ;
    colormap(gray) ;
%    plot(midpos(1),midpos(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppimspots(i,j,1),ppimspots(i,j,2),'xc') ;
        end
    end
    hold off

    [ppim,stdm,pphist,ppminmax] = pepperplotfun('E:\pepperpot\04_05_06\500ms_16-5kV_50.tif',0,[2 0]) ;
    [ppim_35,stdm_35,pphist_35,ppminmax_35] = pepperplotfun('E:\pepperpot\04_05_06\500ms_35kV_50.tif',0,[2 0]) ;
    ppmin = ppminmax(1) ; ppmax = ppminmax(2) ;

    figure, imagesc(ppim,[ppmin ppmax]), axis image, hold on ;
    colormap(gray) ;
%    plot(midpos(1),midpos(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppimspots(i,j,1),ppimspots(i,j,2),'xc') ;
        end
    end
    hold off

    figure, imshow(ppim,[ppmin 150]) ;

    figure, imagesc(ppim,[ppmin 150]), axis image, hold on ;
    colormap(gray) ;
%    plot(midpos(1),midpos(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppimspots(i,j,1),ppimspots(i,j,2),'xc') ;
            plot(ppimspotsb(i,j,1),ppimspotsb(i,j,2),'xr') ;
        end
    end
    hold off

    figure, imagesc(ppim_35,[ppmin 150]), axis image, hold on ;
    colormap(gray) ;
%    plot(midpos(1),midpos(2),'xr') ;    % Mark central point
    for i = 1:25
        for j = 1:25
            plot(ppimspots(i,j,1),ppimspots(i,j,2),'xc') ;
            plot(ppimspotsb(i,j,1),ppimspotsb(i,j,2),'xr') ;
        end
    end
    hold off

    figure, imagesc(ppim,[ppmin+4 250]), axis image, hold on ; colormap(hot) ;

%    pplim = 110 ;
    pplim = 105 ;
    width = 51 ;

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

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            plot(ppimanal(i,j,3),ppimanal(i,j,4),'xk') ;
        end
    end
    hold off

    figure, imagesc(ppim,[ppmin 150]), axis image, hold on ;
    colormap(gray) ;
    for i = 1:25
        for j = 1:25
            plot(ppimspots(i,j,1),ppimspots(i,j,2),'xr') ;
            plot(ppimanal(i,j,1),ppimanal(i,j,2),'xc') ;
        end
    end
    hold off

    [ppimanal_35] = ppallspots(ppim_35,ppimspots(:,:,1:2),ppimspotsb(:,:,1:2),width,pplim) ;

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            plot(ppimanal_35(i,j,3),ppimanal_35(i,j,4),'xk') ;
        end
    end
    hold off

    figure, imagesc(ppim_35,[ppmin 150]), axis image, hold on ;
    colormap(gray) ;
    for i = 1:25
        for j = 1:25
            plot(ppimspots(i,j,1),ppimspots(i,j,2),'xr') ;
            plot(ppimanal_35(i,j,1),ppimanal_35(i,j,2),'xc') ;
        end
    end
    hold off

    rupix = rulen./80 ;
    sppix = (rulen.*mlrulratio)./72 ;
    phospix = sppix.*perscale ;
    zlen = 13 ;
    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) ;
    ppimdata_35 = pprtpconv(ppimspotsoff,zlen,phospix) ;

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

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

%    Intensity data

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

    figure, bar3(ppimdata(:,:,11)) ;
    figure, bar3(ppimdata_35(:,:,11)) ;

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            plot(ppimdata(i,j,1),ppimdata(i,j,2),'xk') ;
        end
    end
    hold off

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            plot(ppimdata(i,j,6),ppimdata(i,j,7),'xk') ;
        end
    end
    hold off

%    Rotate coordinates by 101.25 degrees

    rotang = ( ( 90 + 11.25 )./360 ).*( 2.*pi ) ;

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

    [ppemitdata,ppemit] = ppemitmat(ppimdata(:,:,3:5), ppimdata(:,:,8:10), rotang) ;
    [ppemitdata_35,ppemit_35] = ppemitmat(ppimdata_35(:,:,3:5), ppimdata_35(:,:,8:10), rotang) ;

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            if ppemit(i,j,9) ~= 0
                plot(ppemit(i,j,1),ppemit(i,j,8),'xk') ;
            end
        end
    end
    xlabel(['Horizontal position, x (mm)']), ylabel(['Horizontal divergence, x'' (mrad)']), ...
        title(['\bf Horizontal phase space distribution, x-x'', for 16.5keV Pepperpot image']) ;
    hold off

    figure, hold on ;
    for i = 1:25
        for j = 1:25
            if ppemit(i,j,9) ~= 0
                plot(ppemit(i,j,2),ppemit(i,j,9),'xk') ;
            end
        end
    end
    xlabel(['Vertical position, y (mm)']), ylabel(['Vertical divergence, y'' (mrad)']), ...
        title(['\bf Vertical phase space distribution, y-y'', for 16.5keV Pepperpot image']) ;
    hold off

    [emioutx, betax, gammax, xrms, xprms] = ppemitcalc(ppemitdata(:,:,1), ppemitdata(:,:,3), ppimdata(:,:,end), 16.5e3) ;
    [emiouty, betay, gammay, yrms, yprms] = ppemitcalc(ppemitdata(:,:,2), ppemitdata(:,:,4), ppimdata(:,:,end), 16.5e3) ;
    [emioutx_35, betax_35, gammax_35, xrms_35, xprms_35] = ...
        ppemitcalc(ppemitdata_35(:,:,1), ppemitdata_35(:,:,3), ppimdata_35(:,:,end), 35e3) ;
    [emiouty_35, betay_35, gammay_35, yrms_35, yprms_35] = ...
        ppemitcalc(ppemitdata_35(:,:,2), ppemitdata_35(:,:,4), ppimdata_35(:,:,end), 35e3) ;

    return


