         throughput *= fOverPdf * abs(dot(wiW, nW));
         if (throughput < THROUGHPUT_EPSILON) break;
 
+        // Add direct lighting term at current surface vertex
+        float lightPdf = 0.0;
+        if (!inDielectric)
+        {
+            L += throughput * directSurfaceLighting(pW, basis, woW, hitMaterial, wavelength_nm, XYZ, rnd, inVolume, lightPdf);
+        }
+
         // Prepare for tracing the bounce ray
-        pW += nW * sign(dot(wiW, nW)) * eps; // perturb vertex into half-space of scattered ray
-        float lightPdf = pdfHemisphere(wiL); // compute MIS weight for bounce ray
-        float misWeight = powerHeuristic(bsdfPdf, lightPdf);
+        pW += nW * sign(dot(wiW, nW)) * eps;                 // perturb vertex into half-space of scattered ray
+        float misWeight = powerHeuristic(bsdfPdf, lightPdf); // compute MIS weight for bounce ray
     }
 
     // Compute tristimulus contribution from estimated radiance

            volHits++;
        }

        if (volHits == maxVolHits) hit = false;
        vec3 woW = -rayDir;


float sampleLightInVolume(in vec3 XYZ, inout vec4 rnd, inout vec3 wiW, inout float lightPdf)
{
    // Light sampling (choose either sun or sky)
    vec3 wiW_sun = sampleSunDir(rnd);
    float sunPdf = 1.0; // convenient to decouple total sun power from its angular size
    float sunWeight = sunPower * sunPdf;
    
    float skyPdf;
    vec3 wiW_sky = sampleSphere(rnd, skyPdf);
    float skyWeight = skyPower * skyPdf;

    float sunProb = clamp(max(sunWeight, DENOM_TOLERANCE) / max(skyWeight+sunWeight, DENOM_TOLERANCE), 0.0, 1.0);
    bool chooseSun = (rand(rnd) <= sunProb);
    if (chooseSun)
    {
        lightPdf = sunPdf * sunProb; 
        wiW = wiW_sun;
        return sunRadiance(wiW_sun, XYZ);
    }
    lightPdf = skyPdf * max(PDF_EPSILON, 1.0-sunProb);
    wiW = wiW_sky; 
    return environmentRadiance(wiW, XYZ);
}

// Estimate direct radiance at the given volumetric vertex
