Simulation of a pressure field developed by a song Greetings!

What I am doing:

In short, I want to write a program that is fed a song as a wav file, and spits out an animated contour plot indicating the relative pressure/loudness/quantity indicating amplitude at various points in a plane. Each speaker is modeled as a point source a set distance from each other hanging in an infinite void. My weapon of choice is Matlab.

How I am doing it:

Poorly. Song is read and split into channels, then cut into little slices. Each slice is FFT'd and the amplitudes and phases are split. For each point in the simulation space every frequency is taken, wavelength calculated and used in combination with the phase to determine the contribution of that frequency to the dimensionless loudness/pressure/sound intensity indicating quantity in that cell (the exact unit isn't that important - explained in next section). Do this for each slice of song, draw contour plots for each slice, and animate! 

What my problem is:

The above seems kinda clumsy, and my preliminary results tell me I'm likely neglecting something huge or making some kind of fundamental error. Please note the idea here isn't to make a simulation for speaker setup or any kind of real physical use, I want to draw pretty pictures of what songs do with some basis in reality (specific application: psytrance). This is why the units are not that important, whether I use the normalised voltages in the wav file or convert it to dB or pressure... apart from looking a bit different in lin or log space the relative difference should be the same, at least for my pretty-picture-drawing-purposes. Once the program works as planned some number-fudging will occur to make it look better. 

My Questions

1) Can someone suggest a better way of doing this? I feel like I'm doing this in a too complex manner, but I don't even know where to begin looking for a different method to accomplish my goals.

2) Am I right in doing phase corrections with the values given in the FFT? Gut says the error is here - Fourier coefficients are returned as complex numbers a+bi - for the use I have described here, how would I use this to yield "amplitude" and phase angle? It's a simple question... but humour me, self-kicking might follow your answer.

Edit: On second thought this may not be the most appropriate place to post this (at least for part 1 of the question), if someone can suggest a better place, please do.
For a perfect isotropic point source the pressure or air displacement at a point is just the pressure/displacement stimulating signal delayed by the time it took the disturbance to travel to that point. 

That means you get spherical shells of constant displacement, and the displacement is the signal at the time $t-d/c$ where $d$ is the distance from the source and $c$ is the propagation speed.

Addendum:

If you do not require the full resolution of the pressure field that would allow you to actually reconstruct the sound from the field then you may want to average the pressure over spatial volumes of constant size. That is not equivalent to a time average of the power, because the influence of spatial coarse graining on the pressure history depends on where you are. Only in the far field where the radius of the constant pressure sphere is much larger than the typical coarse graining scale, time averaging ans spatial averaging become equivalent.

So it depends on what you want exactly and what you need it for. The physics involved is very simple for a perfect point source, but gets very complicated as soon as you have participating objects. Maybe you can give a more detailed description of what you are trying to achieve with this.