function [Ksd] = HoopStats_KsDensity(Theta,P,K)
% "HoopStats_KsDensity"
% R02.00
% By Sam Berens (sam.berens@york.ac.uk).
%
% Returns kernel density estimates for a sample of circularly distributed
% data in "P", evaluated at the angles in "Theta". "K" is a kappa value for
% the von Mises distribution. Here, it defines the kernel bandwidth of the
%  that spreads the density function around each point in "P".
%
% [MLE] = HoopStats_KsDensity(Theta,P,K)

%% Use default if "K" is not provided:
if nargin < 3
    K = 2;
end

%% Compute Ksd:
T = Theta(:);
P = P(:)';
D = exp(1i.*T)*exp(-1i.*P); % Distances;
A = angle(D); % Angles;
C = 1/(2*pi*besseli(0,K)); % A constant on K;
Pd = C * exp(K*cos(A));
Ksd = mean(Pd,2);
Ksd = reshape(Ksd,size(Theta));

return