 float sampleLightInVolume(in vec3 XYZ, inout vec4 rnd, inout vec3 wiW, inout float lightPdf)
 {
     // Light sampling (choose either sun or sky)
-    float relSunPower = sunPower/max(skyPower+sunPower, DENOM_TOLERANCE);
-    bool chooseSun = (rand(rnd) < relSunPower);
-    float Li;
-    if (chooseSun) // Sample sun
-    {
-        wiW = sampleSunDir(rnd);
-        lightPdf = relSunPower;
-        Li = sunRadiance(wiW, XYZ);
-    }
-    else // Sample sky
-    {
-        wiW = sampleSphere(rnd, lightPdf);
-        lightPdf *= max(PDF_EPSILON, 1.0-relSunPower);
-        Li = environmentRadiance(wiW, XYZ);
-    }
-    return Li;
+    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


vec3 sampleSunDir(inout vec4 rnd)
{
    float theta = sunAngularSize * M_PI/180.0 * sqrt(rand(rnd));
    float phi = 2.0 * M_PI * rand(rnd);
    float costheta = cos(theta);
    float sintheta = sin(theta);
    float cosphi = cos(phi);
    float sinphi = sin(phi);
    float x = sintheta * cosphi;

        rayDir = wiW; // Update ray direction

        // Detect dielectric transmission
        if (hitMaterial==MAT_DIELE && dot(woW, nW)*dot(wiW, nW) < 0.0) 
        {
            inDielectric = !inDielectric;
        }
