function P=clusterdistrib(r,VAR,D,R)

SD=(2*(pi^(D/2)))/GAMMA(D/2);

if VAR==0
   P=zeros(1,size(r,2));
   P(find(r==0))=inf;
else
   P=SD*(1/((2*pi*VAR)^(D/2)))*(r.^(D-1)).*exp(-((r.^2)/(2*VAR)));
   P(find(r==0))=1;
end

%P(find(r>R))=0;