#define HALF_PI 1.57079632679

surface test (

color SurfaceColor=color(0.8,0.74,0.68); 
float diffuseCoeff=.8; 
float Ka2=0.5; 
color sheenColor2=color(1,1,1); 
float eta2=0.714; 
float thickness2=0.5; 
color SurfaceColor2=color(.8,.5,.5); 
color incandescence=color(0,0,0); 
float translucence=0; 
color specularColor=color(0.5,0.5,0.5); 
float eccentricity=.3; 
float specularRollOff=.7; ) 
{
/** Surface main-code start **/

/* Evaluate the Henyey-Greenstein phase function for two vectors with
an asymmetry value g. v1 and v2 should be normalized and g should 
be in the range (-1, 1). Negative values of g correspond to more
back-scattering and positive values correspond to more forward scattering.
*/
float phase(vector v1, v2; float g) {
float costheta = -v1 . v2;
return (1. - g*g) / pow(1. + g*g - 2.*g*costheta, 1.5);
}

/* Compute a the single-scattering approximation to scattering from
a one-dimensional volumetric surface. Given incident and outgoing
directions wi and wo, surface normal n, asymmetry value g (see above),
scattering albedo (between 0 and 1 for physically-valid volumes),
and the thickness of the volume, use the closed-form single-scattering
equation to approximate overall scattering.
*/
float singleScatter(vector wi, wo; normal n; float g, albedo, thickness) {
float win = abs(wi . n);
float won = abs(wo . n);

return albedo * phase(wo, wi, g) / (win + won) *
(1. - exp(-(1/win + 1/won) * thickness));
}

vector efresnel(vector II; normal NN; float eta; output float Kr, Kt;) {
vector R, T;
fresnel(II, NN, eta, Kr, Kt, R, T);
Kr = smoothstep(0., .5, Kr);
Kt = 1. - Kr;
return normalize(T);
}


/* Implements overall skin subsurface shading model. Takes viewing and
surface normal information, the base color of the skin, a
color for an oily surface sheen, the ratio of the indices of 
refraction of the incoming ray (typically 
1 for air) to the index
of refraction for the transmitted ray (say something like 1.4 for
skin), and the overall thickness of the skin layer. Then loops
over light sources with illuminance() and computes the reflected
skin color.
*/
color subsurfaceSkin(vector Vf; normal Nn; color skinColor, sheenColor;
float eta, thickness) {
extern point P;
float Kr, Kt, Kr2, Kt2;
color C = 0;

vector T = efresnel(-Vf, Nn, eta, Kr, Kt);

illuminance(P, Nn, PI/2) {
vector Ln = normalize(L);

vector H = normalize(Ln + Vf);
if (H . Nn > 0)
C += Kr * sheenColor * Cl * (Ln . Nn) * pow(H . Nn, 4.);
C += Kr * sheenColor * Cl * (Ln . Nn) * .2;

vector T2 = efresnel(-Ln, Nn, eta, Kr2, Kt2);
C += skinColor * Cl * (Ln . Nn) * Kt * Kt2 * 
(singleScatter(T, T2, Nn, .8, .8, thickness) +
singleScatter(T, T2, Nn, .3, .5, thickness) +
singleScatter(T, T2, Nn, 0., .4, thickness));
}
return C;
}

normal Nn2 = faceforward(normalize(N), I);
vector Vf2 = -normalize(I);
/* #define HALF_PI 1.57079632679 */

normal Nf3;
vector H3, Ln3, V3;
color Ia3, Id3, Itr3, Is3;
float NH3, NH23, NHSQ3, Dd3, Gg3, VN3, VH3, LN3, Ff3, tmp3;
/* uniform float E3 = 1 / (eccentricity * eccentricity - 1); */
float E3 = 1 / (eccentricity * eccentricity - 1);

Nf3 = faceforward(normalize(N), I);
Ia3 = ambient() + subsurfaceSkin(Vf2, Nn2, SurfaceColor2, sheenColor2, eta2, thickness2);
Id3 = diffuseCoeff * diffuse(Nf3);

Itr3 = 0;
if (translucence != 0) {
illuminance(P, Nf3, PI)
Itr3 += Cl;
Itr3 *= translucence;
}

Is3 = 0;
V3 = normalize(-I);
VN3 = V3.Nf3;
illuminance(P, Nf3, HALF_PI) {
Ln3 = normalize(L);
H3 = normalize(Ln3+V3);
NH3 = Nf3.H3;
NHSQ3 = NH3*NH3;
NH23 = NH3 * 2;
Dd3 = (E3+1) / (NHSQ3 + E3);
Dd3 *= Dd3;
VH3 = V3.H3;
LN3 = Ln3.Nf3;
if (VN3 < LN3) {
if (VN3 * NH23 < VH3)
Gg3 = NH23 / VH3;
else
Gg3 = 1 / VN3;
}
else {
if (LN3 * NH23 < VH3)
Gg3 = (LN3 * NH23) / (VH3 * VN3);
else
Gg3 = 1 / VN3;
}
tmp3 = pow((1 - VH3), 3);
Ff3 = tmp3 + (1 - tmp3) * specularRollOff;
Is3 += Cl * Dd3 * Gg3 * Ff3;
}
Is3 *= specularColor;

/** Surface main-code end **/

Ci = SurfaceColor * (Ia3 + Id3 + Itr3 + incandescence + Is3);

}