@ DuarteFarrajotaRamos 's comment turns out to have been a good suggestion - ask someone. I found the IRC channel on one of the linked pages, asked, and got the following link straight away. Exactly what I needed!

https://developer.blender.org/diffusion/B/browse/master/intern/cycles/kernel/closure/volume.h

note: it's the .h file that has the goodies I'm after, not the .c. Line 35 has the math (several lines shown here):

/* HENYEY-GREENSTEIN CLOSURE */ /* Given cosine between rays, return probability density that a photon bounces * to that direction. The g parameter controls how different it is from the * uniform sphere. g=0 uniform diffuse-like, g=1 close to sharp single ray. */ ccl_device float single_peaked_henyey_greenstein(float cos_theta, float g) { return ((1.0f - g * g) / safe_powf(1.0f + g * g - 2.0f * g * cos_theta, 1.5f)) * (M_1_PI_F * 0.25f); }; 

Below are some plots of the shape as a function of g, along with the python script to make the plots them. The Henyey-Greenstein dates back to 1941 in a paper analyzing astronomical scattering of dust. The distribution is not profoundly magic, but it has some roots in science and works nicely, with one parameter g from -1 through 0 to +1 reproducing backwards scattering to isotropic to forward scattering.

If you are interested, here are a number of links with a wide variety of depth and focus. Plots are below that.

http://www.oceanopticsbook.info/view/scattering/the_henyeygreenstein_phase_function

http://adsabs.harvard.edu/abs/1941ApJ....93...70H

http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1941ApJ....93...70H&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf

https://www.astro.umd.edu/~jph/HG_note.pdf

http://omlc.org/education/ece532/class3/hg.html

https://en.wikipedia.org/wiki/Monte_Carlo_method_for_photon_transport

http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1985A%26A...146...67H&defaultprint=YES&filetype=.pdf

https://www.cs.dartmouth.edu/~wjarosz/publications/dissertation/chapter4.pdf

http://web.cs.wpi.edu/~emmanuel/courses/cs563/S07/talks/Paulo_volumeScattering_wk10_p1.pdf

def pu(cos_theta, g): top = oneoverfourpi * (1. - g**2) bot = (1 + g**2 - 2*g*cos_theta)**1.5 return top / bot import numpy as np import matplotlib.pyplot as plt halfpi, pi, twopi, fourpi = [f*np.pi for f in [0.5, 1.0, 2.0, 4.0]] oneoverfourpi = 1. / fourpi degs, rads = 180/pi, pi/180 theta = np.linspace(0, pi, 200) cos_theta, sin_theta = np.cos(theta), np.sin(theta) dtheta = theta[1] - theta[0] g = np.arange(0.9, -1, -0.15)[:, None] Pu = pu(cos_theta, g) plt.figure() plt.subplot(3, 1, 1) for thing in Pu: plt.plot(degs*theta, thing) plt.subplot(3, 1, 2) for gval, thing in zip(g, Pu): if np.abs(gval) < 1E-05: plt.plot(degs*theta, thing, '--k') else: plt.plot(degs*theta, thing) plt.yscale('log') plt.subplot(3, 1, 3) for thing in Pu: plt.plot(degs*theta, twopi * sin_theta * dtheta * thing) plt.show()