--- ray/src/rt/srcsamp.c	1991/10/21 12:57:13	1.1
+++ ray/src/rt/srcsamp.c	1994/07/20 12:48:44	2.6
@@ -8,40 +8,38 @@ static char SCCSid[] = "$SunId$ LBL";
  * Source sampling routines
  */
 
-#include  "standard.h"
+#include  "ray.h"
 
-#include  "object.h"
-
 #include  "source.h"
 
 #include  "random.h"
 
 
-extern int  dimlist[];		/* dimension list for distribution */
-extern int  ndims;		/* number of dimensions so far */
-extern int  samplendx;		/* index for this sample */
+static int  cyl_partit(), flt_partit();
 
 
 double
-nextssamp(org, dir, si)		/* compute sample for source, rtn. distance */
-FVECT  org, dir;		/* origin is read only, direction is set */
+nextssamp(r, si)		/* compute sample for source, rtn. distance */
+register RAY  *r;		/* origin is read, direction is set */
 register SRCINDEX  *si;		/* source index (modified to current) */
 {
 	int  cent[3], size[3], parr[2];
 	FVECT  vpos;
 	double  d;
 	register int  i;
-tryagain:
+nextsample:
 	while (++si->sp >= si->np) {	/* get next sample */
 		if (++si->sn >= nsources)
 			return(0.0);	/* no more */
-		if (srcsizerat <= FTINY)
-			nopart(si, org);
+		if (source[si->sn].sflags & SSKIP)
+			si->np = 0;
+		else if (srcsizerat <= FTINY)
+			nopart(si, r);
 		else {
 			for (i = si->sn; source[i].sflags & SVIRTUAL;
 					i = source[i].sa.sv.sn)
 				;		/* partition source */
-			(*sfun[source[i].so->otype].of->partit)(si, org);
+			(*sfun[source[i].so->otype].of->partit)(si, r);
 		}
 		si->sp = -1;
 	}
@@ -71,37 +69,36 @@ tryagain:
 		vpos[i] += (double)cent[i]/MAXSPART;
 					/* compute direction */
 	for (i = 0; i < 3; i++)
-		dir[i] = source[si->sn].sloc[i] +
+		r->rdir[i] = source[si->sn].sloc[i] +
 				vpos[SU]*source[si->sn].ss[SU][i] +
 				vpos[SV]*source[si->sn].ss[SV][i] +
 				vpos[SW]*source[si->sn].ss[SW][i];
 
 	if (!(source[si->sn].sflags & SDISTANT))
 		for (i = 0; i < 3; i++)
-			dir[i] -= org[i];
+			r->rdir[i] -= r->rorg[i];
 					/* compute distance */
-	if ((d = normalize(dir)) == 0.0)
-		goto tryagain;	/* at source! */
+	if ((d = normalize(r->rdir)) == 0.0)
+		goto nextsample;		/* at source! */
 
 					/* compute sample size */
-	si->dom  = source[si->sn].ss2;
 	if (source[si->sn].sflags & SFLAT) {
-		si->dom *= sflatform(si->sn, dir);
-		if (si->dom <= FTINY) {		/* behind source */
-			si->sp = si->np;
-			goto tryagain;
-		}
-		si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART);
+		si->dom = sflatform(si->sn, r->rdir);
+		si->dom *= size[SU]*size[SV]/(MAXSPART*(double)MAXSPART);
 	} else if (source[si->sn].sflags & SCYL) {
-		si->dom *= scylform(si->sn, dir);
-		si->dom *= (double)size[SU]/MAXSPART;
+		si->dom = scylform(si->sn, r->rdir);
+		si->dom *= size[SU]/(double)MAXSPART;
 	} else {
-		si->dom *= (double)(size[SU]*size[SV]*size[SW]) /
-				(MAXSPART*MAXSPART*MAXSPART) ;
+		si->dom = size[SU]*size[SV]*(double)size[SW] /
+				(MAXSPART*MAXSPART*(double)MAXSPART) ;
 	}
-	if (source[si->sn].sflags & SDISTANT)
+	if (source[si->sn].sflags & SDISTANT) {
+		si->dom *= source[si->sn].ss2;
 		return(FHUGE);
-	si->dom /= d*d;
+	}
+	if (si->dom <= 1e-4)
+		goto nextsample;		/* behind source? */
+	si->dom *= source[si->sn].ss2/(d*d);
 	return(d);		/* sample OK, return distance */
 }
 
@@ -137,9 +134,9 @@ unsigned char  *pt;		/* partition array */
 }
 
 
-nopart(si, ro)			/* single source partition */
+nopart(si, r)			/* single source partition */
 register SRCINDEX  *si;
-FVECT  ro;
+RAY  *r;
 {
 	clrpart(si->spt);
 	setpart(si->spt, 0, S0);
@@ -147,9 +144,9 @@ FVECT  ro;
 }
 
 
-cylpart(si, ro)			/* partition a cylinder */
+cylpart(si, r)			/* partition a cylinder */
 SRCINDEX  *si;
-FVECT  ro;
+register RAY  *r;
 {
 	double  dist2, safedist2, dist2cent, rad2;
 	FVECT  v;
@@ -157,11 +154,11 @@ FVECT  ro;
 	int  pi;
 					/* first check point location */
 	clrpart(si->spt);
-	sp = &source[si->sn];
-	rad2 = 1.273 * DOT(sp->ss[SV],sp->ss[SV]);
-	v[0] = ro[0] - sp->sloc[0];
-	v[1] = ro[1] - sp->sloc[1];
-	v[2] = ro[2] - sp->sloc[2];
+	sp = source + si->sn;
+	rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
+	v[0] = r->rorg[0] - sp->sloc[0];
+	v[1] = r->rorg[1] - sp->sloc[1];
+	v[2] = r->rorg[2] - sp->sloc[2];
 	dist2 = DOT(v,sp->ss[SU]);
 	safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
 	dist2 *= dist2 / safedist2;
@@ -171,15 +168,15 @@ FVECT  ro;
 		si->np = 0;
 		return;
 	}
-	safedist2 *= 4./(srcsizerat*srcsizerat);
-	if (dist2 <= 4.*rad2 ||		/* point too close to subdivide? */
-			dist2cent >= safedist2) {
+	safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
+	if (dist2 <= 4.*rad2 ||		/* point too close to subdivide */
+			dist2cent >= safedist2) {	/* or too far */
 		setpart(si->spt, 0, S0);
 		si->np = 1;
 		return;
 	}
 	pi = 0;
-	si->np = cyl_partit(ro, si->spt, &pi, MAXSPART,
+	si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
 			sp->sloc, sp->ss[SU], safedist2);
 }
 
@@ -212,32 +209,45 @@ double  d2;
 	newct[0] = cent[0] - newax[0];
 	newct[1] = cent[1] - newax[1];
 	newct[2] = cent[2] - newax[2];
-	npl = cyl_partit(ro, pt, pi, mp*3/4, newct, newax, d2);
+	npl = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
 					/* upper half */
 	newct[0] = cent[0] + newax[0];
 	newct[1] = cent[1] + newax[1];
 	newct[2] = cent[2] + newax[2];
-	npu = cyl_partit(ro, pt, pi, mp-npl, newct, newax, d2);
+	npu = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
 					/* return total */
 	return(npl + npu);
 }
 
 
-flatpart(si, ro)			/* partition a flat source */
+flatpart(si, r)				/* partition a flat source */
 register SRCINDEX  *si;
-FVECT  ro;
+register RAY  *r;
 {
-	register double  *vp;
+	register FLOAT  *vp;
+	FVECT  v;
 	double  du2, dv2;
 	int  pi;
 
+	clrpart(si->spt);
+	vp = source[si->sn].sloc;
+	v[0] = r->rorg[0] - vp[0];
+	v[1] = r->rorg[1] - vp[1];
+	v[2] = r->rorg[2] - vp[2];
+	vp = source[si->sn].snorm;
+	if (DOT(v,vp) <= FTINY) {	/* behind source */
+		si->np = 0;
+		return;
+	}
+	dv2 = 2.*r->rweight/srcsizerat;
+	dv2 *= dv2;
 	vp = source[si->sn].ss[SU];
-	du2 = 4./(srcsizerat*srcsizerat) * DOT(vp,vp);
+	du2 = dv2 * DOT(vp,vp);
 	vp = source[si->sn].ss[SV];
-	dv2 = 4./(srcsizerat*srcsizerat) * DOT(vp,vp);
-	clrpart(si->spt);
+	dv2 *= DOT(vp,vp);
 	pi = 0;
-	si->np = flt_partit(ro, si->spt, &pi, MAXSPART, source[si->sn].sloc,
+	si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
+		source[si->sn].sloc,
 		source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
 }
 
@@ -282,12 +292,12 @@ double  du2, dv2;
 	newct[0] = cent[0] - newax[0];
 	newct[1] = cent[1] - newax[1];
 	newct[2] = cent[2] - newax[2];
-	npl = flt_partit(ro, pt, pi, mp*3/4, newct, u, v, du2, dv2);
+	npl = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
 					/* upper half */
 	newct[0] = cent[0] + newax[0];
 	newct[1] = cent[1] + newax[1];
 	newct[2] = cent[2] + newax[2];
-	npu = flt_partit(ro, pt, pi, mp-npl, newct, u, v, du2, dv2);
+	npu = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
 				/* return total */
 	return(npl + npu);
 }
@@ -298,7 +308,7 @@ scylform(sn, dir)		/* compute cosine for cylinder's pr
 int  sn;
 register FVECT  dir;		/* assume normalized */
 {
-	register double  *dv;
+	register FLOAT  *dv;
 	double  d;
 
 	dv = source[sn].ss[SU];