How to simply find the code/script behind the volume scatter node so I can see the math?
@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&amp;data_type=PDF_HIGH&amp;whole_paper=YES&amp;type=PRINTER&amp;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&amp;defaultprint=YES&amp;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) &lt; 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()