We wish to solve for in equation 1. is a direction vector and can be described by spherical polar coordinates i.e. . Hence we can view 1 as a five dimensional integro-differential equation.
We discretize equation 1 spatially on a 3D grid, and angularly in either 6 (axial) or 26 (axial+diagonal) directions, using a forward-differencing to approximate the directional derivative.
To evaluate the integral we initially tried to approximate as a sum of spherical harmonics as defined above. Unfortunately this turned out to be a pathological case for interpolating with due to our fixed set of directions falling directly on the lobes of the basis functions causing ill-conditioning, evidenced by zeros in the singular value decomposition of the of the interpolation matrix.
To obtain a smooth interpolation function (see figure 1) we chose a set of independant basis functions composed of products of sine and cosine functions (see appendix A). We want to interpolate a function defined on a sphere -
where k ranges over 0..5 or 0..25. Given our discrete array of function values
and we can find . To integrate the function over the sphere we let
Figure 1: Intensity distribution.
Note that the vector can be precalculated.
While maintaining suitable boundary conditions equation 2 can be iterated until it converges, i.e. the transport of light reaches a steady state.