what is the probability distribution for the angle of an approximate laserbeam

Question

I'm trying to simulate the light distribution characteristics from a Gaussian laser beam, but having difficulty with the angular distribution.

I need to generate a large number of points on an x/y plane, along with pointing vectors down the z-axis, such that their aggregate approximates the power distribution of a laser beam. The distribution on the x/y plane is a Gaussian, while the pointing vector down the z-axis approximates the laser divergence.

For example, a beam with a 1 mm beam waist with a 1.5 mrad divergence has the following beam irradiance, E(r):

$E(r) = exp(-r^2/b^2)/(\pi b^2)$

where $r$ is the radius from the beam center, $b$ is the beam waist (1/e radius). Therefore, I can sample from this distribution by the following equation:

$r = b\sqrt(-ln(1-\mathbb{R}))$

where $\mathbb{R}$ is a random number uniformly distribution on [0,1], to get the radius from the origin on the beams starting point on the x/y plane.

OK, so the question I'm having trouble with, is how do I randomly choose the polar angle $\theta$ of divergence such that it approximates the laser divergence of 1.5 mradians? Do I just choose the polar angle to be uniform on $[0,1.5e^{-3}]$ radians?

I guess I'm getting confused by whether the distribution should be uniform over polar angle, or uniform irradiance over the solid angle, and how to sample from that. Below is an illustration to help sort things out. I'm trying to determine how to distribute the polar angle in order to approximate the laser divergence.

Answer 1

This process is treated in this paper: http://www.springerlink.com/content/2ww1gtp5cvbrerhm/fulltext.pdf

The relevant equations are:

$x_s=\frac{w(-d)}{\sqrt{2}}\text{Erf}^{-1}(2r_1-1)$

$y_s=\frac{w(-d)}{\sqrt{2}}\text{Erf}^{-1}(2r_2-1)$

$x_f=\frac{w_0}{\sqrt{2}}\text{Erf}^{-1}(2r_3-1)$

$y_f=\frac{w_0}{\sqrt{2}}\text{Erf}^{-1}(2r_4-1)$

These give the x and y positions of rays at the beam waist ($x_f$) and at a point far away ($x_s$), with d much greater than the Rayleigh range. $w(z)$ is the standard gaussian beam waist formula, and the $r_i$ are uniform random variates on (0,1).

EDIT: Realized I didn't really address your divergence issue. The above is valid for a collimated diffraction limited beam, so to add excess divergence you could apply the ABCD matrix for a thin lens to the rays in order to create a diverging beam.