function y = circularPdfWrappedPoisson(r, lambda, m, varargin)
%
%< circularPdfWrappedPoisson >
%  Returns pdf value of Wrapped Poisson distribution. This distribution
%  represents the probability of th (steps of rotation) given that average
%  rotation step number is lambda, where one step is 2*pi/m. Therefore, th
%  must be equal to (2*pi/m)*r, where r is integer. In this sense, this
%  function receives r instead of th. If you want to use th, see the
%  example below.
%
%  circularPdfWrappedPoisson(r, lambda, m, varargin)
%       r:   number of circular steps
%       lambda: mean number of steps in a given period (determined by user)
%       m:   number of steps in one circle, i.e., step size is 2*pi/m
%  Optional parameter
%       tol: tolerance level
%       th:  Use this variable if you don't want to use r.
%
%  Example
%       r = 2;
%       lambda = 4;
%       m = 8;
%       y = circularPdfWrappedPoisson(r,lambda,m);
%       tol = circularDefaultTol();
%       y = circularPdfWrappedPoisson(r,lambda,m,tol);
%       th = 2*pi/m*r;
%       y = circularPdfWrappedPoisson(r,lambda,m,tol,th);
%       y = circularPdfWrappedPoisson(NaN,lambda,m,NaN,th);
%
%  Unit is radian.
%


[tol, th] = circularArgChk(varargin);
if isnan(tol)
    tol=circularDefaultTol();
end
if ~isnan(th)
    th = mod(th,2*pi);
    r = floor(th/(2*pi/m))+1;
end
if isnan(r) & isnan(th)
    warning('no input...');
    y = nan;
    return;
end
if isnan(th)
    if mod(r,1)~=0
        warning('r should be integer...');
        y = nan;
        return;
    elseif r>=m | r<0
        warning('r should be an integer in [0,m).');
        y = nan;
        return;
    end
end

tol = tol*exp(lambda)/2;
k = 0;
n = r+k*m;
yic = (lambda.^(n))./(factorial(n));
ytotal = yic;

while sum(sum(abs(yic))) > tol/10
    k = k+1;
    n = r+k*m;
    yic = (lambda.^(n))./(factorial(n));
    ytotal = ytotal+yic;
end

% r
% m
% lambda
% k
% yic
% ytotal
y = ytotal*exp(-lambda);