function vonrand = vonMisesrand(nrand,mu,sigma)
% Usage:   vonrand = vonMisesrand(nrand,mu,sigma)
%  Returns nrand samples from circular normal (vonMises) distribution,
%  centered at mu with standard deviation sigma.
%
% inverse cumulative method, executed by table lookup with
% linear interpolation.  Written by Paul Schrater 2001
%
% C. Olman: 2006 - imbedded vonMisespdf so only have to keep track of one
% file.

if sigma > 0.25,
	% build sampled cdf
	x = (-pi:2*pi/(2e3):pi);
	pofx = vonMisespdf(x,0,sigma);
	cofx = cumsum(pofx/sum(pofx));
	%
	u = rand(1,nrand);
	vonrand = interp1(cofx,x,u)+mu;
	vonrand = min(max(vonrand,-pi+eps),pi-eps);
else
	vonrand = randn(1,nrand)*sigma +mu;
end	


function pofx = vonMisespdf(x,mu,sigma)
% For -pi <= x <= pi

% force x-mu within -pi to pi
y = angle(exp(i*(x-mu)));
kappa = 1/(sigma)^2;

pofx = exp(kappa*cos(y))/(2*pi*besseli(0,kappa));

