Random numbers from a distribution and the Lambert W function

When I was making Physibounce, one of the things that took me a while to figure out is something that is unfortunately barely noticeable in gameplay. When the speed of light becomes small, I wanted the speed distribution of particles to be a 2D Maxwell-Juttner distribution rather than the non-relativistic Maxwell speed distribution. The probability density function for particles allowed to be close to the speed of light is:

f(v) = \frac{\gamma(v)^4 m^2 c^2 v}{k_B^2 T^2(1+\frac{mc^2}{k_B T})} \exp\left(-\frac{mc^2}{k_B T}\left(\sqrt{1+\gamma(v)^2\frac{v^2}{c^2}}-1\right)\right).

Where v is the velocity of the particle, \gamma(v) = (1-\frac{v^2}{c^2})^{-1/2} is the Lorentz factor, k_B is the Boltzmann constant, m is the mass of the mass of the particle, T is the temperature of the system, and c is the speed of light. This distribution will look like a Maxwell distribution in the c\to\infty limit, and is normalized so that \int_0^c f(v) dv = 1.

The simplest way to generate random numbers following a distribution is the acceptance-rejection method. This is what I do in Physibounce when the speed of light is not changed. However, I was not able to find an expression for a suitable (and efficient) bound (or maximum value) for the distribution given above. Because of this, I decided to use inverse sampling instead. To do this, I first find the cumulative distribution function as

 F = \int_0^v f(v') dv'.

Amazingly, Mathematica was able to do this analytically. By the definition, it is clear F takes on values in the interval [0,1]. Thus, if I find a relation that expressed v in terms of F (equivalent to inverting the function), I have an expression that takes uniform values on the interval [0,1] and maps them to a velocity in the interval [0,c] that follows the desired distribution. Mathematica was also able to do the inversion, but the result was in terms of the Lambert W function.

The Lambert W function is the function that satisfies

 z = W(z) e^{W(z)}.

So it seems reasonable something like this would come up when inverting the cumulative distribution function. An important thing about this function is that it is double valued! There are two branches. When I asked Mathematica to solve the problem, it chose a branch and the branch it chose turned out to be the wrong one. I wondered for a while why I was not getting reasonable velocity values when I plugged in random numbers, and after a bit I realized I needed to choose the other branch.

The next problem I ran into was that I did know of a way to actually evaluate the Lambert W function in actionscript. So I turned to the asymptotic expansions. It turned out I was in a regime where the argument of the Lambert W tended to be large. I took the first few terms in the expansion about \infty for the Lambert W. I verified that over the parameter ranges I could expect in Physibounce that the difference between the Lambert W and the expansion was less than 1%. Since it's pretty hard to look at the speeds of a bunch of objects and determine the distribution of their speeds, I thought this was more than good enough, though I did include a check to make sure no particles were created with speed greater than c. In the end the code ended up only being a few lines, but quite a bit of thought went into each of the lines!

See the ActionScript code here.