--- ray/src/rt/source.c	1995/12/08 18:22:28	2.19
+++ ray/src/rt/source.c	1997/08/15 10:24:24	2.27
@@ -22,6 +22,10 @@ static char SCCSid[] = "$SunId$ LBL";
 
 extern double  ssampdist;		/* scatter sampling distance */
 
+#ifndef MAXSSAMP
+#define MAXSSAMP	16		/* maximum samples per ray */
+#endif
+
 /*
  * Structures used by direct()
  */
@@ -223,6 +227,7 @@ register RAY  *r;
 			}
 
 	if (r->ro != NULL) {
+		r->robj = objndx(r->ro);
 		for (i = 0; i < 3; i++)
 			r->ron[i] = -r->rdir[i];
 		r->rod = 1.0;
@@ -333,8 +338,7 @@ char  *p;			/* data for f */
 				( sr.ro != source[scp->sno].so ||
 				source[scp->sno].sflags & SFOLLOW )) {
 						/* follow entire path */
-			if (!raycont(&sr))
-				objerror(sr.ro, USER, "material not found");
+			raycont(&sr);
 			rayparticipate(&sr);
 			if (trace != NULL)
 				(*trace)(&sr);	/* trace execution */
@@ -374,70 +378,73 @@ char  *p;			/* data for f */
 srcscatter(r)			/* compute source scattering into ray */
 register RAY  *r;
 {
+	int  oldsampndx;
 	int  nsamps;
 	RAY  sr;
 	SRCINDEX  si;
-	double  t, lastt, d;
-	COLOR  cumval, ctmp;
+	double  t, d;
+	double  re, ge, be;
+	COLOR  cvext;
 	int  i, j;
 
-	if (r->slights == NULL || r->slights[0] == 0 || r->gecc >= 1.-FTINY)
+	if (r->slights == NULL || r->slights[0] == 0
+			|| r->gecc >= 1.-FTINY || r->rot >= FHUGE)
 		return;
 	if (ssampdist <= FTINY || (nsamps = r->rot/ssampdist + .5) < 1)
 		nsamps = 1;
-	initsrcindex(&si);
+#if MAXSSAMP
+	else if (nsamps > MAXSSAMP)
+		nsamps = MAXSSAMP;
+#endif
+	oldsampndx = samplendx;
+	samplendx = random()&0x7fff;		/* randomize */
 	for (i = r->slights[0]; i > 0; i--) {	/* for each source */
-		setcolor(cumval, 0., 0., 0.);
-		lastt = r->rot;
-		for (j = nsamps; j-- > 0; ) {	/* for each sample position */
+		for (j = 0; j < nsamps; j++) {	/* for each sample position */
 			samplendx++;
 			t = r->rot * (j+frandom())/nsamps;
+							/* extinction */
+			re = t*colval(r->cext,RED);
+			ge = t*colval(r->cext,GRN);
+			be = t*colval(r->cext,BLU);
+			setcolor(cvext,	re > 92. ? 0. : exp(-re),
+					ge > 92. ? 0. : exp(-ge),
+					be > 92. ? 0. : exp(-be));
+			if (intens(cvext) <= FTINY)
+				break;			/* too far away */
 			sr.rorg[0] = r->rorg[0] + r->rdir[0]*t;
 			sr.rorg[1] = r->rorg[1] + r->rdir[1]*t;
 			sr.rorg[2] = r->rorg[2] + r->rdir[2]*t;
 			sr.rmax = 0.;
-						/* sample ray to this source */
-			if (si.sp >= si.np-1 || !srcray(&sr, NULL, &si) ||
-					sr.rsrc != r->slights[i]) {
-				si.sn = r->slights[i]-1;	/* reset */
-				si.np = 0;
-				if (!srcray(&sr, NULL, &si) ||
-						sr.rsrc != r->slights[i])
-					continue;		/* no path */
-			}
+			initsrcindex(&si);	/* sample ray to this source */
+			si.sn = r->slights[i];
+			nopart(&si, &sr);
+			if (!srcray(&sr, NULL, &si) ||
+					sr.rsrc != r->slights[i])
+				continue;		/* no path */
 			copycolor(sr.cext, r->cext);
-			sr.albedo = r->albedo;
+			copycolor(sr.albedo, r->albedo);
 			sr.gecc = r->gecc;
+			sr.slights = r->slights;
 			rayvalue(&sr);			/* eval. source ray */
 			if (bright(sr.rcol) <= FTINY)
 				continue;
-							/* compute fall-off */
-			d = lastt - t;
-			setcolor(ctmp,	1.-d*colval(r->cext,RED),
-					1.-d*colval(r->cext,GRN),
-					1.-d*colval(r->cext,BLU));
-			multcolor(cumval, ctmp);
-			lastt = t;
 			if (r->gecc <= FTINY)		/* compute P(theta) */
 				d = 1.;
 			else {
 				d = DOT(r->rdir, sr.rdir);
-				d = sqrt(1. + r->gecc*r->gecc - 2.*r->gecc*d);
-				d = (1. - r->gecc*r->gecc) / (d*d*d);
+				d = 1. + r->gecc*r->gecc - 2.*r->gecc*d;
+				d = (1. - r->gecc*r->gecc) / (d*sqrt(d));
 			}
 							/* other factors */
-			d *= si.dom * r->albedo * r->rot / (4.*PI*nsamps);
+			d *= si.dom * r->rot / (4.*PI*nsamps);
 			multcolor(sr.rcol, r->cext);
+			multcolor(sr.rcol, r->albedo);
 			scalecolor(sr.rcol, d);
-			addcolor(cumval, sr.rcol);
+			multcolor(sr.rcol, cvext);
+			addcolor(r->rcol, sr.rcol);	/* add it in */
 		}
-						/* final fall-off */
-		setcolor(ctmp,	1.-lastt*colval(r->cext,RED),
-				1.-lastt*colval(r->cext,GRN),
-				1.-lastt*colval(r->cext,BLU));
-		multcolor(cumval, ctmp);
-		addcolor(r->rcol, cumval);	/* sum into ray result */
 	}
+	samplendx = oldsampndx;
 }