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
hence
where
and we can find 
. To integrate the function
over the sphere we let
where
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.
