% beam.m

w = 1;               % arbitrary choice of beam width. 

% CREATE SET OF RANDOM NUMBERS THAT SPECIFY SET OF r1 VALUES
N = 10000;           % choose a set of 10,000 random numbers
RND = rand(1,N);     % random number generator yields RND
r1 = w*sqrt(-log(RND)); 
                     % choice of r1 based on RND (from Monte Carlo expression)

% CREATE HISTOGRAM AND NORMALIZE TO MATCH p(r)
nb = 30;             % choose 30 bins for histogram
[n,x] = hist(r1,nb); % n(x) = histo(x=r1) using nb bins
dx = x(3)-x(2);      % dx is the bin size
p = n/N/dx;          % normalize n by N and dx to yield prob. density fcn. p
sum(p*dx)            % verify that integral of p over all x equals unity. YES.
pr = p./(2*pi*x);    % normalize by circumference of ring at r=x. 
                     % This is pr(r) generated by random number generator.

% CREATE ANALYTIC EXPRESSION FOR GAUSSIAN, pr(r)
r = (0:.1*w:3*w);    % domain of r values
y = exp(-r.^2/w^2)/(pi*w^2);  
                     % pure Gaussian, such that 
                     % Integral[ y (2 pi r) dr, for r = 0 to Infinity] = 1

% DRAW GRAPH COMPARING HISTOGRAM (pr) AND ANALYTIC (y).
figure
plot(x,(pr),'yo')
xlabel('r')
ylabel('p(r)/(2 pi r)')
hold on
plot(r,(y),'r-')
title(['w = ' int2str(w)])
hold off