All the foregoing can be generalized from Rayleigh scattering to Mie scattering cases. Since we seldom are dealing with particles of a single size, we must consider scattering by a distribution of particle sizes. Here is J code to compute the scattering cross-sections and elements of the scattering matrix for power-law distributions of particle sizes: Power_law.ijs. The equations for the distributions considered can be found in the appendix of Harrington, Lame, White and Borkowski (1997) AJ 113, 2147. (Use Polar_plot.ijs to see the resulting scattering phase functions.)

As an example, we can take the old Rouleau & Martin [1991: ApJ 377,526] optical constants for carbon (benzene burned in air, "BE1"). We will adopt the MRN power law distribution: N(a) proportional to the grain radius "a" to the -3.5 power, for "a" between 0.005 and 0.25 microns. We consider a wavelength of 5000A. Here is the polar plot of the resulting phase function:BE1-MRN. In this figure, the blue line is for incoming radiation polarized perpendicular to the scattering plane, while the red line is for incoming radiation polarized in the scattering plane. We see the distribution is strongly forward scattering (by a ratio of 12.6). The green line is the fractional polarization of the scattered radiation, which never exceeds 60 percent -- the maximum polarization occurs at an angle slightly greater than 90 degrees.

Here is the code for a beam into a slab, like "Pol_beam.ijs" given above, modified to use the scattering matrix from our power-law distribution:PB_rj.ijs. The equations are given in sections 4 and 5 of notes3.pdf. Since this code uses the "rejection method", we encounter the problem mentioned in section 5 of the notes: too many trials of angles greater than 90 degrees will be rejected because the scattering is strongly forward scattering. This code thus splits the tests between forward and back scattering to diminish this problem. (In addition, we include the Stokes parameter V, since some circular polarization can arise depending on the scattering matrix.)

As an example of a "PB_rj" run, let's consider a slab of total optical thickness tau=2, with the opacity due to the distribution of carbon particles illustrated above. We illuminate this slab with a beam of unpolarized light entering with angle of about 45 deg (mu= cos theta= 0.7). The albedo of this distribution is 0.3076, so most of the incident light will be absorbed. We consider many scatterings, though with this albedo, only the first few will be significant. We find that 4.39% of the radiation is scattered back out the top of the slab, exp(-2/0.7) = 5.74% passes through without absorption, and 3.70% is scattered out the bottom of the slab. Thus 13.8% escaped the slab while 86.2% is absorbed.
Here is the intensity distribution of the scattered radiation as a function of azimuth (y-axis) and zenith distance (x-axis - 0 deg is perpendicular to the surface the slab, 90 deg tangent). The beam enters at azimuth = 0 and zenith distance = 45. Scattered from the top face: C-int-up. From the lower face: C-int-down. The spike is where radiation which passes through without interaction emerges at azimuth = 180 and z.d. = 45. (This spike should be much higher, but is suppressed to show the scattered radiation.) Note how the downward scattered radiation is concentrated around that spike due to the strong forward scattering of this particle distribution. The upward scattered radiation is strongest for z.d. near 90 deg and azimuths near 180: this is where the angle of scattering is greatest.
Here is the fractional polarization of the scattered radiation from the upper face C-pol-up and from the lower: C-pol-down. Finally, we have the angle of polarization for the two faces: C-ang-up ; C-ang-down.

Here we present variation of "PB_rj.ijs" which uses inverse interpolation to find the "theta" scattering angles: PB_mixed.ijs. This code should be good even if the scattering matrix has a strongly peaked phase function. (Still a bug in this one: after running for a day or so, it may halt on a NaN :-( ... otherwise, seems OK). Further, here is a version of the code for a point source in a slab, PPaS.ijs, which uses the inverse interpolation method for a general scattering matrix: PPaSmix.ijs. In addition, here is code for a point source on the axis on a cylinder: Cylmix.ijs. And here is the code for an external beam scattered by an infinite slab (generalized from Pol_beam.ijs): PB_Mie.ijs. (Note that these codes also do Rayleigh or electron scattering: just feed them the pre-computed Rayleigh scattering matrix generated by Rayleigh_table.ijs.)