Radiative Transfer

by Created: 01 Jun 1996 Updated: 24 May 2014

Radiative transfer is the process by which radiation propagates through space. Radiation is absorbed and emitted by interstellar media. At NASA, I implemented a radiative transfer Monte Carlo simulation – effectively a a volumetric backwards ray tracer. We used this to predict results from Hubble Space Telescope observations.

RT, as we called it, introduced me to volumetric ray tracing and Monte Carlo fluid dynamics simulation, and also the importance of high-quality random number generation. I developed RT in C++, which had just recently added templates to the language. Susan had also developed a protostellar density model in Fortran. RT called the density model with points (in a different coordinate system), and the model returned a density value.

Other astronomers at IPAC were stunned that we ran these simulations on a single Mac desktop computer. At the time, astronomy software mainly ran on Sun Solaris workstations, and simulations typically ran on supercomputers or very beefy workstations (or multiple machines). It was the very early days of commodity PCs beginning to be competitive with workstation-class machines.

project/rt/rt_process.png

In computer graphics, ray tracing normally works backwards from the pixels you see using Kajiya’s rendering equation. In physics, ray tracing usually starts with the light source and then propagates forward to the observer. This “wastes” most of the work performed because most rays go off in a direction other than what your camera is looking at. However, we were able to exploit rotational symmetry in the model to be able to recapture most of the photons traced.

project/rt/adaptive_grid.png

Other astronomy simulations often made gross simplifications to raycasting – one recent paper at the time limited rays to 27 directions chosen to match grid orientations. RT handled arbitrary directions (uniformly random on a sphere, or biased to match scattering physics). RT used an adaptive grid implemented with a BSP tree, and I developed a novel technique for efficiently iterating through ray intersections points with the grid. RT used this grid to model a astronomical-sized volume of space with sufficient resolution near areas of interest and still fit into available memory. It modeled photon absorption and scattering, ionization of interstellar gases, and time-delayed recombination to ground state (emitting new photos). Ionization and recombination make the simulation stateful and time-dependent.

Monte Carlo simulations require good random number generators. I initially implemented the generalized feedback shift register (GFSR) random number generator GFSR(250,103), also known as R250, with a period of almost 2250. I later combined this with GFSR(521,168) to overcome statistical defects, and later still, replaced it with the newly discovered Mersenne Twister GFSR(624,397) which was slower than R250 but faster than R250+521, and has a period of length 219937+1. Today, (SF)MT is still popular, and the WELL generator is state-of-the-art.

RT simulated the Strömgren sphere, giving some confidence in the accuracy of the simulation. We used RT to simulate approximately what the Hubble Space Telescope would see if it observed a protostellar system that matched Susan’s density model. Susan used these theoretical results to get HST observing time, leading to the published papers below.

One of these days, I hope to dust off the source code and update it for modern times. This problem is especially amenable to GPGPU computation and parallelization in a distributed computing environment.

Publications

Additional Resources