function [ppemitdata,ppalldata] = ppemitmat(holertp,spotrtp,rotang)
% Function PPEMITMAT.M
%
%    [ppemitdata,ppalldata] = ppemitmat(holertp,spotrtp,rotang)
%
%    PPEMITMAT creates a matrix, PPEMITDATA, that contains all the
%    necessary coordinate data for calculating emittance, using data
%    extracted from Pepperpot images: predicted hole postions (in both
%    Cartesian [X,Y] coordinates and spherical polar [r,theta,phi]
%    coordinates), measured spot positions (again in both Cartesian [X,Y]
%    coordinates and spherical polar [r,theta,phi] coordinates) and
%    emittance angles (in Cartesian [X',Y'] coordinates).
%
%    As before, the spherical polar coordinates are defined such that theta
%    is the azimuthal angle from the x-axis (ie. around the z-axis) and phi
%    the polar angle (ie. angle away from the z-axis).  The angle theta=0
%    corresponds to the x-z plane and phi=0 to the z-axis.  The units for
%    the spherical polar coordinates must be (mm,rads,rads).  All the lines
%    in [r,theta,phi] should originate on the tungsten grid.
%
%    [PPEMITDATA] = PPEMITMAT(HOLERTP,SPOTRTP) converts the [r,theta,phi]
%    information in the two input matrices, HOLERTP and SPOTRTP, into X and
%    Y emittance data contained in the matrix PPEMITDATA.  HOLERTP must
%    contain [r,theta,phi] data for the predicted hole positions; SPOTRTP
%    must contain [r,theta,phi] data for the measured spot positions;
%    PPEMITDATA contains the [X,X'] and [Y,Y'] emittance data.  The
%    dimensions/arrangement of the matrices are as follows:
%
%        HOLERTP - An [M x N x 3] matrix of hole locations, in
%                  [r,theta,phi], relative to the central hole position.
%                  The grid must be M holes by N holes (usually 25 x 25);
%                  [M,N,1] gives the r-distance to the centre of hole
%                  (M,N), in mm; [M,N,2] and [M,N,3] give the theta and phi
%                  angles, respectively, in radians.
%
%        SPOTRTP - An [M x N x 3] matrix of measured spot positions, in
%                  [r,theta,phi], relative to the measured hole position.
%                  The grid must be M holes by N holes (usually 25 x 25);
%                  [M,N,1] gives the r-distance to the centre of spot
%                  (M,N), in mm; [M,N,2] and [M,N,3] give the theta and phi
%                  angles, respectively, in radians.
%
%        PPEMITDATA - An [M x N x 4] matrix, where (M,N) corresponds to the
%                     spot position as before; [M,N,1:2] gives the
%                     predicted [X,Y] position of each hole, in (mm);
%                     [M,N,3:4] gives the [X',Y'] emittance angle of each
%                     measured spot, in (mrad).
%
%    [PPEMITDATA] = PPEMITMAT(HOLERTP,SPOTRTP,ROTANG) also rotates the
%    cooordinate system of both the hole and spot positions by an angle
%    ROTANG radians, to account for any rotational misalignment of the
%    pepperpot setup.
%
%    [PPEMITDATA,PPALLDATA] = PPEMITMAT(HOLERTP,SPOTRTP,ROTANG) produces a
%    second output matrix, PPALLDATA, that contains all the information
%    from PPEMITDATA, HOLERTP and SPOTRTP, with any additional rotation,
%    into a single [M x N x 12] matrix that is structured in the following
%    way:
%
%        1: X distance from grid centre to hole, in mm (after rotation) - X
%        2: Y distance from grid centre to hole, in mm (after rotation) - Y
%        3: r distance from grid centre to hole, in mm
%        4: theta of line from grid centre to hole, in rads (after rotation)
%        5: phi of line from grid centre to hole, in rads
%        6: X distance from hole to spot centre, in mm (after rotation)
%        7: Y distance from hole to spot centre, in mm (after rotation)
%        8: X' angle from hole to spot centre, in mrad (after rotation) - X'
%        9: Y' angle from hole to spot centre, in mrad (after rotation) - Y'
%        10: r distance from hole to spot centre, in mm
%        11: theta of line from hole to spot centre, in rads (after rotation)
%        12: phi of line from hole to spot centre, in rads
%
%    As such, PPEMITDATA = PPALLDATA(:,:,[1:2,8:9]).
%
%    When specifying the input matrices, it is important to note that the
%    definition for the origins of the HOLERTP and SPOTRTP matrices is
%    different.  All the lines in HOLERTP must originate on the centre of
%    the tungsten grid ie. measurements are relative to the position of the
%    central hole.  In SPOTRTP, however, the line for each spot originates
%    at the corresponding hole on the tungsten grid, NOT the grid centre
%    ie. measurements are relative to the corresponding grid hole.
%

    if nargin < 2
        error('Need at least 2 input variables: [ppemitdata,ppalldata] = ppemitmat(holertp,spotrtp)') ;
    elseif nargin > 3
        error('Can only specify 3 input variables: [ppemitdata,ppalldata] = ppemitmat(holertp,spotrtp,rotang)') ;
    elseif nargout < 1
        error('Must specify at least 1 output variable: [ppemitdata] = ppemitmat(holertp,spotrtp,rotang)') ;
    elseif nargout > 2
        error('Can only specify 2 output variables: [ppemitdata,ppalldata] = ppemitmat(holertp,spotrtp,rotang)') ;
    end

    if nargin < 3
        rotang = 0 ;
    end

%    Make sure holertp and spotrtp are the same size

    hmsize = size(holertp) ;
    spmsize = size(spotrtp) ;

    if ( hmsize(1) ~= spmsize(1) ) | ( hmsize(2) ~= spmsize(2) ) | ( hmsize(3) ~= spmsize(3) ) | ( hmsize(3) ~= 3 )
        error('holertp and spotrtp must be [M x N x 3] matrices of the same size')
    end

%    Create new matrix for all data

    ppalldata = zeros(hmsize(1),hmsize(2),12) ;

%    Rotate and insert data from input matrices

    ppalldata(:,:,3) = holertp(:,:,1) ;             % r of holes
    ppalldata(:,:,4) = holertp(:,:,2) - rotang ;    % theta of holes after rotation
    ppalldata(:,:,5) = holertp(:,:,3) ;             % phi of holes

    ppalldata(:,:,10) = spotrtp(:,:,1) ;             % r of spots
    ppalldata(:,:,11) = spotrtp(:,:,2) - rotang ;    % theta of spots after rotation
    ppalldata(:,:,12) = spotrtp(:,:,3) ;             % phi of spots

%    Calculate rotated X & Y for hole positions and spot offsets

    ppalldata(:,:,1) = ppalldata(:,:,3).*cos(ppalldata(:,:,4)).*sin(ppalldata(:,:,5)) ;
    ppalldata(:,:,2) = ppalldata(:,:,3).*sin(ppalldata(:,:,4)).*sin(ppalldata(:,:,5)) ;

    ppalldata(:,:,6) = ppalldata(:,:,10).*cos(ppalldata(:,:,11)).*sin(ppalldata(:,:,12)) ;
    ppalldata(:,:,7) = ppalldata(:,:,10).*sin(ppalldata(:,:,11)).*sin(ppalldata(:,:,12)) ;

%    Calculate emittance angles from theta and phi

    ppalldata(:,:,8) = atan( cos(ppalldata(:,:,11)).*tan(ppalldata(:,:,12)) ).*1000 ;
    ppalldata(:,:,9) = atan( sin(ppalldata(:,:,11)).*tan(ppalldata(:,:,12)) ).*1000 ;

%    Create ppemitdata matrix

    ppemitdata = ppalldata(:,:,1:2) ;
    ppemitdata(:,:,3:4) = ppalldata(:,:,8:9) ;

    return

