             volHits++;
         }
 
+        // If ray was absorbed by the volume, terminate
+        if (absorb) break;
         if (volHits == maxVolHits) hit = false;
         vec3 woW = -rayDir;
 

            rayLength = length(pW_next - pW);
        }

        /*  Volume rendering logic:
            ray state: pW -> pW_next,  hitMaterial(pW_next), inVolume

            while (inVolume || hitMaterial(pW_next)==VOLUME)

                - if inVolume
                    - if !inDielectric, check for scatter on segment pW -> pW_next
                    - if no scatter, 
                            - set pW = pW_next
                            - if hitMaterial is surface, break and handle it (note, still inside volume)
                            - if hitMaterial is volume boundary, 
                                - displace pW into volume exterior 
                                - inVolume = false
                            - (or if !hit, that's a numerical glitch as we should be inside a closed volume boundary, so just break)
                    - else if scatter,
                            - set pW = pScatter
                            - update rayDir via phase function sampling      

                - else
                    - advance to volume boundary:  pW = pW_next
                    - displace pW into volume interior
                    - inVolume = true
                    
                - re-trace from new pW, giving new pW_next, hit, hitMaterial(pW_next)
                    
            - if !hit (ray escaped to infinity) add env light contibution, and terminate
            - otherwise, hit a surface. Add contribution, fire bounce ray pW -> pW_next, continue
        */


        // Deal with all volume interactions until the ray either hits a surface or escapes
        int volHits = 0;
        while ((inVolume || hitMaterial==MAT_VOLUM) && volHits<maxVolHits)
        {
            // If ray lies in the volume interior

            {
                float segmentLength = length(pW_next - pW);

                // Do Woodcock tracking to check for a possible volume absorption/scattering event along the in-volume segment pW -> pW_next
                // (except if the current ray lies inside a dielectric, in which case there is no volumetric scattering, i.e. the dielectric displaces the volume)
                vec3 pScatter;
                float sigma_t;
                bool interacted = false;
                if (!inDielectric)
                {
                    float distanceMarched = -log(rand(rnd)) * inv_sigma_t_max;
                    int steps = 0; 
                    while (distanceMarched<segmentLength-eps && steps<maxWoodcockSteps)
                    {
                        pScatter = pW + distanceMarched*rayDir;
                        sigma_t = VOLUME_EXTINCTION(pScatter);
                        interacted = bool(rand(rnd) < sigma_t * inv_sigma_t_max);
                        if (interacted) break;
                        distanceMarched += -log(rand(rnd)) * inv_sigma_t_max;
                        steps++;
                    }
                }
            
                // Absorption/scattering event
                if (interacted)
                {
                    vec3 wiW = samplePhaseFunction(rayDir, rnd); // sample scattered dir (NB, PF is the angle PDF so the MC PDF denom. cancels it in the L estimator)
                    float emission = abs(dot(RGB, VOLUME_EMISSION(pScatter)));
                    vec3 XYZ_albedo = rgbToXyz(VOLUME_ALBEDO(pScatter));
                    vec3 XYZ_spec = xyz_to_spectrum(XYZ_albedo);
                    float albedo = dot(XYZ, XYZ_spec);
                    absorb = (rand(rnd) > albedo);
                    if (absorb)
                    {
                        L += throughput * emission; // add volume emission term
                        break;
                    }

                    L += throughput * directVolumeLighting(pScatter, -rayDir, XYZ, rnd, inVolume); // add contribution due to direct lighting at vertex
                    throughput *= albedo; // update throughput due to scattering
                    if (throughput < THROUGHPUT_EPSILON) break;
                    pW = pScatter; // continue to next vertex
                    rayDir = wiW;
