% generate a rician distribution for 4 K values

sig=3;
sig2=sig^2;
A=5;
r=[0:1500]/100;
figure;
% Create Rician distribution
ex=exp(-(r.^2+A^2)/(2*sig2));
I=besseli(0,A*r/sig2);
p=r.*ex.*I/sig2;
plot(r,p,'k')
hold on
xlabel('r')
ylabel('p(r)')
axis([0 15 0 .2])
x=[r(500) r(500:1500) r(1500) r(500)];
y=[0 p(500:1500) 0 0];
fill(x,y,[.9 .9 .9])
xlabel('r (mV)')
ylabel('p(r)')