% 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) ;
    [ppimspotss] = ppspotreco(ppcalmat,midp,rulen,imtilt,invperscale) ;

    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

    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(ppimspotss(i,j,1),ppimspotss(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) ;

%    break

%    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 = zeros(25,25,5) ;

    [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

%    break

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

%    Calculate r,theta,phi position for holes and spots
%    ppimdata matrix is ordered as follows:
%
%    Horz distance of hole to centre
%    Vert distance of hole to centre
%    r distance from centre to hole
%    theta of line from centre to hole
%    phi of line from centre to hole
%    Horz distance from hole to spot centre
%    Vert distance from hole to spot centre
%    r distance from hole to spot centre
%    theta of line from hole to spot centre
%    phi of line from hole to spot centre
%    Normalised light intensity (compared to central spot)

    ppimdata = zeros(25,25,11) ;

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

    ppimdata(:,:,1) = (ppimspots(:,:,1) - ppimspots(13,13,1))./phospix ;            %   X
    ppimdata(:,:,2) = (ppimspots(:,:,2) - ppimspots(13,13,2))./(-phospix) ;         %   Y
    ppimdata(:,:,3) = sqrt( ppimdata(:,:,1).^2 + ppimdata(:,:,2).^2 + zlen.^2 ) ;   %   r
    ppimdata(:,:,4) = atan2( ppimdata(:,:,2),ppimdata(:,:,1) ) ;                    %   theta
    ppimdata(:,:,5) = acos( zlen./ppimdata(:,:,3) ) ;                               %   phi

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

    ppimdata(:,:,6:7) = ppimanal(:,:,3:4)./phospix ;                                %   X (& Y)
%    ppimdata(:,:,6:7) = ppimanal_35(:,:,3:4)./phospix ;                             %   X (& Y)
    ppimdata(:,:,7) = ppimdata(:,:,7).*(-1) ;                                       %   Y
    ppimdata(:,:,8) = sqrt( ppimdata(:,:,6).^2 + ppimdata(:,:,7).^2 + zlen.^2 ) ;   %   r
    ppimdata(:,:,9) = atan2( ppimdata(:,:,7),ppimdata(:,:,6) ) ;                    %   theta
    ppimdata(:,:,10) = acos( zlen./ppimdata(:,:,8) ) ;                              %   phi
    ppimdata(:,:,11) = ppimanal(:,:,5) ; %.*invintense ;                            %   Intensity
%    ppimdata(:,:,11) = ppimanal_35(:,:,5) ; %.*invintense ;                         %   Intensity

    figure, bar3(ppimdata(:,:,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

%    Check to make sure r,theta,phi matrices are calculated correctly

    recalctest = 1 ;

    if recalctest

        ppimrecalc = zeros(25,25,5) ;
        ppimrecalc(:,:,1) = ppimdata(:,:,3).*cos(ppimdata(:,:,4)).*sin(ppimdata(:,:,5)) ;
        ppimrecalc(:,:,2) = ppimdata(:,:,3).*sin(ppimdata(:,:,4)).*sin(ppimdata(:,:,5)) ;
        ppimrecalc(:,:,3) = ppimdata(:,:,3).*cos(ppimdata(:,:,5)) ;
        ppimrecalc(:,:,4) = ppimdata(:,:,8).*cos(ppimdata(:,:,9)).*sin(ppimdata(:,:,10)) ;
        ppimrecalc(:,:,5) = ppimdata(:,:,8).*sin(ppimdata(:,:,9)).*sin(ppimdata(:,:,10)) ;

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

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

    end

    break

%    Rotate coordinates by 101.25 degrees

    rotang = ( ( 90 + 11.25 )./360 ).*( 2.*pi ) ;
    pprotdata = ppimdata ;
%    pprotdata(:,:,4) = pprotdata(:,:,4) - rotang ;
%    pprotdata(:,:,9) = pprotdata(:,:,9) - rotang ;

    rotcalc = 1 ;

    if rotcalc

        ppimrecalc = zeros(25,25,5) ;
        ppimrecalc(:,:,1) = pprotdata(:,:,3).*cos(pprotdata(:,:,4)).*sin(pprotdata(:,:,5)) ;
        ppimrecalc(:,:,2) = pprotdata(:,:,3).*sin(pprotdata(:,:,4)).*sin(pprotdata(:,:,5)) ;
        ppimrecalc(:,:,3) = pprotdata(:,:,3).*cos(pprotdata(:,:,5)) ;
        ppimrecalc(:,:,4) = pprotdata(:,:,8).*cos(pprotdata(:,:,9)).*sin(pprotdata(:,:,10)) ;
        ppimrecalc(:,:,5) = pprotdata(:,:,8).*sin(pprotdata(:,:,9)).*sin(pprotdata(:,:,10)) ;

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

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

    end

%    Calculate r,theta,phi for rotated positions and calculate angles for
%    emittance
%    ppemit matrix is ordered as follows:
%
%    Horz distance of hole to centre (after rotation)
%    Vert distance of hole to centre (after rotation)
%    r distance from centre to hole
%    theta of line from centre to hole
%    phi of line from centre to hole
%    Horz distance from hole to spot centre (after rotation - x)
%    Vert distance from hole to spot centre (after rotation - y)
%    Horz angle from hole to spot centre (after rotation - x')
%    Vert angle from hole to spot centre (after rotation - y')
%    r distance from hole to spot centre
%    theta of line from hole to spot centre
%    phi of line from hole to spot centre
%    Normalised light intensity (compared to central spot)


    ppemit = pprotdata(:,:,1:7) ;
    ppemit(:,:,1) = pprotdata(:,:,3).*cos(pprotdata(:,:,4)).*sin(pprotdata(:,:,5)) ;
    ppemit(:,:,2) = pprotdata(:,:,3).*sin(pprotdata(:,:,4)).*sin(pprotdata(:,:,5)) ;
    ppemit(:,:,6) = pprotdata(:,:,8).*cos(pprotdata(:,:,9)).*sin(pprotdata(:,:,10)) ;
    ppemit(:,:,7) = pprotdata(:,:,8).*sin(pprotdata(:,:,9)).*sin(pprotdata(:,:,10)) ;
    ppemit(:,:,8) = atan( cos(pprotdata(:,:,9)).*tan(pprotdata(:,:,10)) ).*1000 ;
    ppemit(:,:,9) = atan( sin(pprotdata(:,:,9)).*tan(pprotdata(:,:,10)) ).*1000 ;
    ppemit(:,:,10:13) = pprotdata(:,:,8:11) ;

    ppangtest = ppemit(:,:,8:9) ;
    ppangtest(:,:,1) = (atan2(pprotdata(:,:,6),zlen)).*1000 ;
    ppangtest(:,:,2) = (atan2(pprotdata(:,:,7),zlen)).*1000 ;
    ppangtest(:,:,3) = ppangtest(:,:,1) - ppemit(:,:,8) ;
    ppangtest(:,:,4) = ppangtest(:,:,2) - ppemit(:,:,9) ;
    xpdiff = sum(sum(ppangtest(:,:,3))) ;
    ypdiff = sum(sum(ppangtest(:,:,4))) ;

    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') ;
%                plot(ppemit(i,j,1),ppangtest(i,j,1),'xr') ;
            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') ;
%                plot(ppemit(i,j,2),ppangtest(i,j,2),'xr') ;
            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(ppemit(:,:,1), ppemit(:,:,8), ppemit(:,:,13), 16.5e3) ;
    [emiouty, betay, gammay, yrms, yprms] = ppemitcalc(ppemit(:,:,2), ppemit(:,:,9), ppemit(:,:,13), 16.5e3) ;
%    [emioutx, betax, gammax, xrms, xprms] = ppemitcalc(ppemit(:,:,1), ppemit(:,:,8), ppemit(:,:,13), 35e3) ;
%    [emiouty, betay, gammay, yrms, yprms] = ppemitcalc(ppemit(:,:,2), ppemit(:,:,9), ppemit(:,:,13), 35e3) ;

    break

    ppemit = zeros(25,25,5) ;
    ppemit(:,:,1:2) = ppimdata(:,:,1:2) ;
    ppemit(:,:,3) = (atan2(ppimdata(:,:,3),zlen)).*1000 ;
    ppemit(:,:,4) = (atan2(ppimdata(:,:,4),zlen)).*1000 ;
    ppemit(:,:,5) = ppimdata(:,:,5) ;

    [emioutx, betax, gammax, xrms, xprms] = ppemitcalc(ppemit(:,:,1), ppemit(:,:,3), ppemit(:,:,5), 16.5e3) ;
    [emiouty, betay, gammay, yrms, yprms] = ppemitcalc(ppemit(:,:,2), ppemit(:,:,4), ppemit(:,:,5), 16.5e3) ;

    return


