function out=Dm(x,M,N)
%
%  function out=Dm(x,M,N)
%
% compute M-order Dirichlet function with period N at x, which can be a vector
%
%   out=sin(pi*(2*M+1)*x/N) ./ (N*sin(pi*x/N))
%

% note: all arguments are normally integers, but routine will work for 
% non-integer values
%
% witten by D.G. Long at BYU 15 Dec 2005
%
% The following commented-out method is faster in Matlab, but it generates
% an annoying error message (division by zero) even though it works properly
%   out=sin(pi*(2*M+1)*x/N) ./ (N*sin(pi*x/N))
%   ind=find(x==0);
%   out(ind)=(2*M+1)/N;
% so, the following method is used instead
%

out=zeros(size(x))+(2*M+1)/N;
num=sin(pi*(2*M+1)*x/N);
denom=sin(pi*x/N);
ind=find(denom~=0);
out(ind)=num(ind)./(N*denom(ind));
