Here is a short description of how the computations were performed. As described above, the initial conditions were formed by developing a set of coefficients for the cosine functions in the direction and the Bessel functions in the direction. Then at each time step, and for each point in the grid, the sum in equation (23) was computed. The programming was done in the Interactive Data Language from Research Systems Incorporated. The script that I used can be found here.
The computations are considerable. To save time, I precomputed the and terms and stored them in an array for easy access. Even so, it must be remembered that a computation of terms must occur for each frame, and there were about 100 frames. Total running time was about three hours on a Pentium II 450 MHz Linux PC. The final data array was saved and then rendered in a second step.