--- ray/src/rt/pmap.c	2018/06/07 19:26:04	2.15
+++ ray/src/rt/pmap.c	2018/11/21 19:33:30	2.16
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: pmap.c,v 2.15 2018/06/07 19:26:04 rschregle Exp $";
+static const char RCSid[] = "$Id: pmap.c,v 2.16 2018/11/21 19:33:30 rschregle Exp $";
 #endif
 
 
@@ -13,7 +13,7 @@ static const char RCSid[] = "$Id: pmap.c,v 2.15 2018/0
        supported by the Swiss National Science Foundation (SNSF, #147053)
    ======================================================================
    
-   $Id: pmap.c,v 2.15 2018/06/07 19:26:04 rschregle Exp $
+   $Id: pmap.c,v 2.16 2018/11/21 19:33:30 rschregle Exp $
 */
 
 
@@ -52,7 +52,8 @@ static int photonParticipate (RAY *ray)
    int i;
    RREAL cosTheta, cosPhi, du, dv;
    const float cext = colorAvg(ray -> cext),
-               albedo = colorAvg(ray -> albedo);
+               albedo = colorAvg(ray -> albedo),
+               gecc2 = ray -> gecc * ray -> gecc;
    FVECT u, v;
    COLOR cvext;
 
@@ -60,6 +61,16 @@ static int photonParticipate (RAY *ray)
    ray -> rmax = -log(pmapRandom(mediumState)) / cext;
    
    while (!localhit(ray, &thescene)) {
+      if (!incube(&thescene, ray -> rop)) {
+         /* Terminate photon if it has leaked from the scene */
+#ifdef DEBUG_PMAP
+         fprintf(stderr, 
+                 "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", 
+                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
+#endif                 
+         return 0;
+      }
+         
       setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]),
                       exp(-ray -> rmax * ray -> cext [1]),
                       exp(-ray -> rmax * ray -> cext [2]));
@@ -75,7 +86,7 @@ static int photonParticipate (RAY *ray)
       /* Store volume photons unconditionally in mist to also account for
          direct inscattering from sources */
       if (albedo > FTINY)
-#endif
+#endif      
          /* Add to volume photon map */
          newPhoton(volumePmap, ray);
          
@@ -88,12 +99,11 @@ static int photonParticipate (RAY *ray)
       scalecolor(ray -> rcol, 1 / albedo);    
       
       /* Scatter photon */
-      cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1
-                                      : 1 / (2 * ray -> gecc) * 
-                                            (1 + ray -> gecc * ray -> gecc - 
-                                               (1 - ray -> gecc * ray -> gecc) / 
-                                               (1 - ray -> gecc + 2 * ray -> gecc *
-                                                  pmapRandom(scatterState)));
+      cosTheta = ray -> gecc <= FTINY 
+                    ? 2 * pmapRandom(scatterState) - 1
+                    : 0.5 * (1 + gecc2 - 
+                         (1 - gecc2) / (1 - ray -> gecc + 2 * ray -> gecc * 
+                         pmapRandom(scatterState))) / ray -> gecc;
                                                   
       cosPhi = cos(2 * PI * pmapRandom(scatterState));
       du = dv = sqrt(1 - cosTheta * cosTheta);   /* sin(theta) */
@@ -115,16 +125,16 @@ static int photonParticipate (RAY *ray)
       ray -> rlvl++;
       ray -> rmax = -log(pmapRandom(mediumState)) / cext;
    }  
-    
+   
+   /* Passed through medium until intersecting local object */  
    setcolor(cvext, exp(-ray -> rot * ray -> cext [0]),
                    exp(-ray -> rot * ray -> cext [1]),
                    exp(-ray -> rot * ray -> cext [2]));
                    
    /* Modify ray color and normalise */
    multcolor(ray -> rcol, cvext);
-   colorNorm(ray -> rcol);
-   
-   /* Passed through medium */  
+   colorNorm(ray -> rcol);   
+
    return 1;
 }