--- ray/src/rt/raytrace.c	2005/04/14 18:04:12	2.47
+++ ray/src/rt/raytrace.c	2005/06/23 09:11:38	2.53
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: raytrace.c,v 2.47 2005/04/14 18:04:12 greg Exp $";
+static const char RCSid[] = "$Id: raytrace.c,v 2.53 2005/06/23 09:11:38 greg Exp $";
 #endif
 /*
  *  raytrace.c - routines for tracing and shading rays.
@@ -13,6 +13,7 @@ static const char RCSid[] = "$Id: raytrace.c,v 2.47 20
 #include  "source.h"
 #include  "otypes.h"
 #include  "otspecial.h"
+#include  "random.h"
 
 #define  MAXCSET	((MAXSET+1)*2-1)	/* maximum check set size */
 
@@ -36,14 +37,22 @@ static void checkset(OBJECT  *os, OBJECT  *cs);
 
 extern int
 rayorigin(		/* start new ray from old one */
-	register RAY  *r,
-	register RAY  *ro,
+	RAY  *r,
 	int  rt,
-	double  rw
+	const RAY  *ro,
+	const COLOR rc
 )
 {
-	double	re;
-
+	double	rw, re;
+						/* assign coefficient/weight */
+	if (rc == NULL) {
+		rw = 1.0;
+		setcolor(r->rcoef, 1., 1., 1.);
+	} else {
+		rw = intens(rc);
+		if (rc != r->rcoef)
+			copycolor(r->rcoef, rc);
+	}
 	if ((r->parent = ro) == NULL) {		/* primary ray */
 		r->rlvl = 0;
 		r->rweight = rw;
@@ -55,10 +64,11 @@ rayorigin(		/* start new ray from old one */
 		copycolor(r->albedo, salbedo);
 		r->gecc = seccg;
 		r->slights = NULL;
-	} else if (ro->rot >= FHUGE) {		/* illegal continuation */
-		memset(r, 0, sizeof(RAY));
-		return(-1);
 	} else {				/* spawned ray */
+		if (ro->rot >= FHUGE) {
+			memset(r, 0, sizeof(RAY));
+			return(-1);		/* illegal continuation */
+		}
 		r->rlvl = ro->rlvl;
 		if (rt & RAYREFL) {
 			r->rlvl++;
@@ -78,14 +88,36 @@ rayorigin(		/* start new ray from old one */
 		r->crtype = ro->crtype | (r->rtype = rt);
 		VCOPY(r->rorg, ro->rop);
 		r->rweight = ro->rweight * rw;
-						/* estimate absorption */
+						/* estimate extinction */
 		re = colval(ro->cext,RED) < colval(ro->cext,GRN) ?
 				colval(ro->cext,RED) : colval(ro->cext,GRN);
 		if (colval(ro->cext,BLU) < re) re = colval(ro->cext,BLU);
-		if (re > 0.)
-			r->rweight *= exp(-re*ro->rot);
+		re *= ro->rot;
+		if (re > 0.1)
+			if (re > 92.)
+				r->rweight = 0.0;
+			else
+				r->rweight *= exp(-re);
 	}
 	rayclear(r);
+	if (r->rweight <= 0.0)			/* check for expiration */
+		return(-1);
+	if (r->crtype & SHADOW)			/* shadow commitment */
+		return(0);
+	if (maxdepth <= 0 && rc != NULL) {	/* Russian roulette */
+		if (minweight <= 0.0)
+			error(USER, "zero ray weight in Russian roulette");
+		if (maxdepth < 0 && r->rlvl > -maxdepth)
+			return(-1);		/* upper reflection limit */
+		if (r->rweight >= minweight)
+			return(0);
+		if (frandom() < r->rweight/minweight)
+			return(-1);
+		rw = minweight/r->rweight;	/* promote survivor */
+		scalecolor(r->rcoef, rw);
+		r->rweight = minweight;
+		return(0);
+	}
 	return(r->rlvl <= maxdepth && r->rweight >= minweight ? 0 : -1);
 }
 
@@ -122,10 +154,10 @@ raytrace(			/* trace a ray and compute its value */
 	} else if (sourcehit(r))
 		rayshade(r, r->ro->omod);	/* distant source */
 
-	rayparticipate(r);		/* for participating medium */
-
 	if (trace != NULL)
 		(*trace)(r);		/* trace execution */
+
+	rayparticipate(r);		/* for participating medium */
 }
 
 
@@ -147,7 +179,7 @@ raytrans(			/* transmit ray as is */
 {
 	RAY  tr;
 
-	if (rayorigin(&tr, r, TRANS, 1.0) == 0) {
+	if (rayorigin(&tr, TRANS, r, NULL) == 0) {
 		VCOPY(tr.rdir, r->rdir);
 		rayvalue(&tr);
 		copycolor(r->rcol, tr.rcol);
@@ -209,10 +241,10 @@ rayparticipate(			/* compute ray medium participation 
 		ge *= 1. - colval(r->albedo,GRN);
 		be *= 1. - colval(r->albedo,BLU);
 	}
-	setcolor(ce,	re<=0. ? 1. : re>92. ? 0. : exp(-re),
-			ge<=0. ? 1. : ge>92. ? 0. : exp(-ge),
-			be<=0. ? 1. : be>92. ? 0. : exp(-be));
-	multcolor(r->rcol, ce);			/* path absorption */
+	setcolor(ce,	re<=FTINY ? 1. : re>92. ? 0. : exp(-re),
+			ge<=FTINY ? 1. : ge>92. ? 0. : exp(-ge),
+			be<=FTINY ? 1. : be>92. ? 0. : exp(-be));
+	multcolor(r->rcol, ce);			/* path extinction */
 	if (r->crtype & SHADOW || intens(r->albedo) <= FTINY)
 		return;				/* no scattering */
 	setcolor(ca,
@@ -305,7 +337,7 @@ raymixture(		/* mix modifiers */
 
 extern double
 raydist(		/* compute (cumulative) ray distance */
-	register RAY  *r,
+	register const RAY  *r,
 	register int  flags
 )
 {
@@ -319,6 +351,32 @@ raydist(		/* compute (cumulative) ray distance */
 }
 
 
+extern void
+raycontrib(		/* compute (cumulative) ray contribution */
+	double  rc[3],
+	const RAY  *r,
+	int  flags
+)
+{
+	double	eext[3];
+	int	i;
+
+	eext[0] = eext[1] = eext[2] = 0.;
+	rc[0] = rc[1] = rc[2] = 1.;
+
+	while (r != NULL && r->crtype&flags) {
+		for (i = 3; i--; ) {
+			rc[i] *= colval(r->rcoef,i);
+			eext[i] += r->rot * colval(r->cext,i);
+		}
+		r = r->parent;
+	}
+	for (i = 3; i--; )
+		rc[i] *= (eext[i] <= FTINY) ? 1. :
+				(eext[i] > 92.) ? 0. : exp(-eext[i]);
+}
+
+
 extern double
 raynormal(		/* compute perturbed normal for ray */
 	FVECT  norm,
@@ -365,7 +423,7 @@ newrayxf(			/* get new tranformation matrix for ray */
 		FULLXF  xf;
 	}  xfseed = { &xfseed }, *xflast = &xfseed;
 	register struct xfn  *xp;
-	register RAY  *rp;
+	register const RAY  *rp;
 
 	/*
 	 * Search for transform in circular list that