function [ Z ] = Gen2DDirichletKern( xiVec,etaVec,Tx,Ty,Px,Py,N1,N2,M1,M2,Oj,Ol  )
% Generates a 2D shifted Dirichlet Kernel.
%   xiVec - vector of xi values to evaluate over
%
%   etaVec - vector of eta values to evaluate over
%
%   Tx - x unit cell periodicity, i.e., x-distance between unit cells in
%       terms of pitch, must be an integer
%
%   Tx - x unit cell periodicity, i.e., x-distance between unit cells in
%       terms of pitch, must be an integer
%
%   Px,Py - the pixel pitch in the x and y directions respectively
%
%   N1,N2 - starting and ending unit cell position (in integers) along x.
%       Depends on the number of unit cells in the x direction.
%
%   M1,M2 - starting and ending unit cell position (in integers) along y
%       Depends on the number of unit cells in the y direction.
%
%   Oj,Ol - the shift away from the (0,0) pixel within a unit cell along x
%       and y respectively, in integer units.

[XI, ETA] = meshgrid(xiVec,etaVec);

ArgXi = Tx*Px.*pi.*XI;
ArgEta = Ty*Py.*pi.*ETA;

% ArgXi(ArgXi == 0) = 10^-25;
% ArgEta(ArgEta == 0) = 10^-25;

DKxiT1 = sin((N1+N2+1).*ArgXi);
DKxiT2 = sin(ArgXi);
% Use L'Hopital's rule at 0/0
DKxiT1(abs(DKxiT2) < 10^-10) = (N1+N2+1).*cos((N1+N2+1).*ArgXi(abs(DKxiT2) < 10^-10));
DKxiT2(abs(DKxiT2) < 10^-10) = cos(ArgXi(abs(DKxiT2) < 10^-10));
DKxi = exp(1i*(N1-N2+2*Oj/Tx).*ArgXi).*DKxiT1./DKxiT2;

DKetaT1 = sin((M1+M2+1).*ArgEta);
DKetaT2 = sin(ArgEta);
DKetaT1(abs(DKetaT2) < 10^-10) = (M1+M2+1).*cos((M1+M2+1).*ArgEta(abs(DKetaT2) < 10^-10));
DKetaT2(abs(DKetaT2) < 10^-10) = cos(ArgEta(abs(DKetaT2) < 10^-10));
DKeta = exp(1i*(M1-M2+2*Ol/Ty).*ArgEta).*DKetaT1./DKetaT2;
Z = DKxi.*DKeta;
end

