% flatbeam.m
%   Steven L. Jacques, Dec. 1998. Oregon Medical Laser Center.

% A MATLAB file which calculates the radial positions for
% launching photons to yield a uniform-irradiance
% circular beam of diameter w (radius a).
%   (1) Creates histogram of Irradiance E(r) [mm^-2]
%   (2) Creates analytic expression for E(r) = 1/(pi*a^2)


w = 1;                 % arbitrary choice of beam diameter. 
a = w/2;               % radius
rmax = w;
r = (0:rmax/30:rmax);  % domain of r values

N = 30000;             % Number of photons launched
RND = rand(1,N);       % Array of random numbers generated
y = 2*r/a^2;           % Array of y = p(r) for flat circular beam,
                       %   such that integrate[p(r), {r,0, Infinity}] = 1.
E = 1/(pi*a^2);        % Constant value of irradiance E(r) [mm^-1].
r1 = a*sqrt(RND);      % Choosing launch positions r1 based on RND

nb = 30;               % number of bins
[n,x] = hist(r1,nb);   % n(x) = histogram(r1) using nb bins
dx = x(3)-x(2);        % incremental bin size
n = n/N/dx;            % Normalize by Ndx.
sum(n*dx)              % Verify that this sum equals unity.
nr = n./(2*pi*x);      % Normalize by circumference of ring at r=x.
                       %   This histogram nr is the same as irradiance E(r).

figure
plot(x, nr, 'yo')      % Plot histogram
xlabel('Radial position, r [mm]')
ylabel('Irradiance, E(r) [mm^-2]  =  1/(Ľa^2)')
hold on
plot([0, rmax], [E, E], 'r-')       % Plot analytic E(r)
title(['flat circular beam, dia = ' int2str(w) ' [mm]'])
axis([0 rmax 0 2*E]);
hold off