@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

The function is already available in OSL. This should not be a surprise since it is so fundamental to rendering and scattering simulations of all varieties! This is from v1.7 of the Open Shading Language 1.7 Language Specification - which is outdated. A new and/or better link to 1.8 or whatever is newest is appreciated!

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()