How to simply find the code/script behind the volume scatter node so I can see the math?

Question

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?

---

Answer 1

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