Writing a Monte Carlo radiative transfer simulation for computing the observed blackbody radiance due to the environment

Question

Background

I'm attempting to determine the power observed at some antenna immersed in some perhaps highly anisotropic scattering and absorbing medium due to the blackbody radiance of the environment. To do so, I've begun constructing a simple Monte Carlo simulation which (when finished) will trace rays isotropically distributed in spherical coordinates $\theta$ and $\phi$ from the antenna out to some boundary drawn around the antenna (whose size is whatever is deemed sufficient to capture all non-negligible contributions). I choose to trace rays backwards from the antenna (as opposed to from their "actual" starting point somewhere in the environment) since, in general, the refractive index of the medium will be a function of position, and often there are large regions of the medium from which light will never reach the antenna. By tracing rays backwards, I expect to achieve significant improvement in computational efficiency for many media.

Here's what my "algorithm" looks like so far:

1. I generate some large number of random directions in $\theta$ and $\phi$.
2. I break the antenna passband $B$ into some large number $M$ of discrete "slices", and associate with each slice its center frequency. For each of the random directions, I trace the stochastic path of some large number $N$ rays for each "slice" in frequency $M$ (i.e. a total of $M\times N$ rays for each random direction).
3. Tracing the stochastic path consists of: a) taking a small Runge-Kutta raytracing step (due to refraction), b) computing the probability of absorbing or scattering (a function of the path length traversed since the last interaction), c) either absorbing/scattering in some new direction, or repeating (or quitting, if we have reached the bounds of the simulation)

I track the cumulative loss in intensity for each ray by simply accumulating the loss at each interaction (absorption/scatter).

Question

I'm not quite sure where to go from here. I'd like to take into account the antenna's frequency response and directivity, and it isn't clear to me how I should go about turning the cumulative intensity loss computed for each ray to something I can multiply by the directivity function and integrate over. Moreover, it isn't clear to me how to best compute the actual intensity due to the environmental blackbody radiance. Should I simply compute the blackbody radiance power at each stochastic step, multiply it by the cumulative loss at that point, and sum these up?

Where do I go from here? Is this an appropriate way to treat the problem?