Bremsstrahlung

In one of my classes (High Energy Astrophysics) I must do a semester project. The assignment was rather open-ended, and for my project i proposed to do a Monte-Carlo Simulation. So far so good. The trouble arises when I consider the starting point of the simulation, which happens to be black body radiation. Black body radiation is emitted according to a Planck distribution (with the physical constants combined into constants a and b):

It's not as though it's that ugly. Now, all I need to do is generate random numbers according to this distribution. For such a task, I turned to Numerical Recipes in C++, where i learned that the easiest way, if the math is tractable, to generate randomly distributed numbers with a new distribution is to apply the inverse of the anti-derivative of that distribution function to a uniform random distribution. It works marvelously for the exponential distribution, which I also happen to need for my project.

This seemed like a method I should certainly try; after all, I know the functional form of my distribution, and it's not all that complicated. First stop: The CRC Standard mathematical tables to look for that integral. Unsettlingly, no integral of such a form appears in the table¹. That means it must be either hard or uncommon. Hopefully it's the latter. Next, ask Maxima what it knows about such a thing:

integrate((x^3)/(exp(x)-1), x);

Uh, oh. What are those lis? Polylogarithms, apparently. Not only do I not know how to compute those, I definitely don't know how to invert them, much less this entire thing. So much for the analytical approach.

The other method recommended by Numerical Recipes is to use a Rejection Method. It's not quite as elegant, but it can get the job done pretty well for tougher distributions. The only tricky point of using the rejection method, is that one has to start with some easily computable distribution, from which to reject points. That starting distribution must have finite area beneath it, lie everywhere above the actually desired distribution, and have an invertible anti-derivative. So far, I think that by best bet for such a bounding distribution is an exponential distribution; it'll be a bit wasteful, but hopefully not too badly so. The trouble is that I haven't yet been able to analytically determine the parameters for such a distribution so as to guarantee that it will always be greater than the planck distribution, adjusting as needed for changes in a and b. Just trying to find the peak of the planck distribution hasn't been pretty, it seems to be at:

That isn't much help anyway, since I need my distribution to pass not through that point, but above it and all other points on the planck distribution.