How to simply find the code/script behind the volume scatter node so I can see the math? tl;dr I need help finding the code or OSL script for the volume scatter node. Not asking how it works, just a specific, noob-friendly instruction how to find it, or a link to it. Blender is big, and this is my first time, so...



The volume scatter node has an Anisotropy input to control the relative strength of scattering in different directions. When changed from +1 to -1 the scattering changes from more forward to more backward scattering.

It turns out there is a published excerpt from the Blender Cycles Encyclopedia that shows this - I've included two screenshots below.

The (unnormalized) distribution could be as simple as:

I(theta) ~ 1 + A * cos(theta)


where A is the Anisotropy factor, or it could be completely different.

I could build a science experiment in blender to "measure" the intensity distribution as a function of angle. I don't know if there is a "light meter" node yet, but I could do it with cameras, saving rendered images, and analyzing them with NumPy.

But I am wondering if I can just read the math somewhere? Are nodes like this written in OSL, or in C? (I learned about OSL in this helpful answer.) How would I go about finding the code for this node to find how the shape of the distribution is calculated?
@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()