/gaussianBeam.m

 1% beam.m
 2
 3% A MATLAB program
 4%   Steven L. Jacques, Dec. 1998.
 5%
 6% Generates the probability density function p(r)
 7% and the relative irradiance E(r) [mm^-2]
 8% for a 2D Gaussian beam profile associated with
 9% a laser beam. Histograms based on 10,000 random numbers
10% are created, and compared with the analytic expressions.
11%
12% Two figures are generated: p(r) and E(r).
13%
14
15close
16close
17clear
18
19b = 1;                       % arbitrary choice of 1/e beam radius. 
20
21N = 10000;                   % number of random numbers used
22rnd1 = rand(1,N);            % create array of random numbers
23r1 = b*sqrt(-log(1 - rnd1)); % array of selected values r1
24                             % Note: log() is base e
25
26nb = 30;                     % number of bins
27[n,x] = hist(r1,nb);         % n(x) = histo(x=r1) using nb bins
28dx = x(3)-x(2);              % bin size
29n = n/N/dx;                  % normalize by N and dx
30sum(n*dx)                    % verify that integral equals unity
31nr = n./(2*pi*x);            % normalize by circumference of ring at r=x
32
33r = (0:.1*b:3*b);            % create an array of r values
34E = exp(-r.^2/b^2)/(pi*b^2); % calculate an array of E values,
35                             %   using a pure 2D Gaussian.
36
37% Create figure of p(r)
38figure
39plot(x, n,'yo')
40hold on
41plot(r, E*2*pi.*r,'r-')
42xlabel('r')
43ylabel('p(r)')
44xlabel('r [mm]')
45title(['b = ' int2str(b) ' [mm]'])
46hold off
47
48% Create figure of E(r)
49figure
50plot(x,nr,'o')
51xlabel('r')
52ylabel('p(r)')
53hold on
54plot(r,E,'r-')
55title(['b = ' int2str(b) ' [mm]'])
56ylabel('E(r) [mm^-2]')
57xlabel('r [mm]')
58hold off