function out=ApDm2d(x,y,M1,M2,N1,N2,U)
%
% function out=ApDm2d(x,y,M1,M2,N1,N2,U)
%
% compute 2d aperture filtered Dirichlet function
%   out=sum_{k1=-L1}^L1 sum_{k2=-L2}^L2 U(k1,k2) Dm(x+k1,y_k2,M1,N1,M2,N2)
%      =sum_{k1=-L1}^L1 sum_{k2=-L2}^L2 U(k1,k2) 
%          Dm(x+k1,M1,N1) Dm(y+k2,M2,N2)
% over mod N1,N2 one period where U is sized (L2,L1) (odd)
% and x vales are in (1..N1) and y values in (1..N2)
%
% uses the routine Dm2d.m
%

% written by D.G. Long at BYU 30 Dec. 2010
% note: this code is designed to merely illustrate the computation.  No
% claims are made for its efficiency
%

if prod(size(U)) == 1  % special case of ideal aperture
  out=U(1)*Dm2d(x,y,M1,M2,N1,N2);
else                   % other apertures
  L=(size(U)-1)/2;     % get aperture array size (assumed to be small)
  k1=-L(2):L(2); k2=-L(1):L(1);
  [kk1, kk2]=meshgrid(k1,k2);
  X=mod(x-kk1(1)+N1,N1);
  Y=mod(y-kk2(1)+N2,N2);
  out=U(1,1)*Dm2d(X,Y,M1,M2,N1,N2);
  for l=2:prod(size(U))
    if (U(l)~=0) % speed up for binary apertures
      X=mod(x-kk1(l)+N1,N1);
      Y=mod(y-kk2(l)+N2,N2);
      out=out+U(l)*Dm2d(X,Y,M1,M2,N1,N2);
    end
  end
end