--- ray/src/rt/srcsamp.c	1991/10/21 14:27:36	1.2
+++ ray/src/rt/srcsamp.c	2009/06/06 02:11:44	2.16
@@ -1,47 +1,47 @@
-/* Copyright (c) 1991 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: srcsamp.c,v 2.16 2009/06/06 02:11:44 greg Exp $";
 #endif
-
 /*
  * Source sampling routines
+ *
+ *  External symbols declared in source.h
  */
 
-#include  "standard.h"
+#include "copyright.h"
 
-#include  "object.h"
+#include  "ray.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];
+	SRCREC  *srcp;
 	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;
 	}
@@ -52,13 +52,14 @@ tryagain:
 	if (!skipparts(cent, size, parr, si->spt))
 		error(CONSISTENCY, "bad source partition in nextssamp");
 					/* compute sample */
+	srcp = source + si->sn;
 	if (dstrsrc > FTINY) {			/* jitter sample */
 		dimlist[ndims] = si->sn + 8831;
 		dimlist[ndims+1] = si->sp + 3109;
 		d = urand(ilhash(dimlist,ndims+2)+samplendx);
-		if (source[si->sn].sflags & SFLAT) {
+		if (srcp->sflags & SFLAT) {
 			multisamp(vpos, 2, d);
-			vpos[2] = 0.5;
+			vpos[SW] = 0.5;
 		} else
 			multisamp(vpos, 3, d);
 		for (i = 0; i < 3; i++)
@@ -69,43 +70,69 @@ tryagain:
 
 	for (i = 0; i < 3; i++)
 		vpos[i] += (double)cent[i]/MAXSPART;
+					/* avoid circular aiming failures */
+	if ((srcp->sflags & SCIR) && (si->np > 1 || dstrsrc > 0.7)) {
+		FVECT	trim;
+		if (srcp->sflags & (SFLAT|SDISTANT)) {
+			d = 1.12837917;		/* correct setflatss() */
+			trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]);
+			trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]);
+			trim[SW] = 0.0;
+		} else {
+			trim[SW] = trim[SU] = vpos[SU]*vpos[SU];
+			d = vpos[SV]*vpos[SV];
+			if (d > trim[SW]) trim[SW] = d;
+			trim[SU] += d;
+			d = vpos[SW]*vpos[SW];
+			if (d > trim[SW]) trim[SW] = d;
+			trim[SU] += d;
+			if (trim[SU] > FTINY*FTINY) {
+				d = 1.0/0.7236;	/* correct sphsetsrc() */
+				trim[SW] = trim[SV] = trim[SU] =
+						d*sqrt(trim[SW]/trim[SU]);
+			} else
+				trim[SW] = trim[SV] = trim[SU] = 0.0;
+		}
+		for (i = 0; i < 3; i++)
+			vpos[i] *= trim[i];
+	}
 					/* compute direction */
 	for (i = 0; i < 3; i++)
-		dir[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];
+		r->rdir[i] = srcp->sloc[i] +
+				vpos[SU]*srcp->ss[SU][i] +
+				vpos[SV]*srcp->ss[SV][i] +
+				vpos[SW]*srcp->ss[SW][i];
 
-	if (!(source[si->sn].sflags & SDISTANT))
+	if (!(srcp->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);
-	} else if (source[si->sn].sflags & SCYL) {
-		si->dom *= scylform(si->sn, dir);
-		si->dom *= (double)size[SU]/MAXSPART;
+	if (srcp->sflags & SFLAT) {
+		si->dom = sflatform(si->sn, r->rdir);
+		si->dom *= size[SU]*size[SV]/(MAXSPART*(double)MAXSPART);
+	} else if (srcp->sflags & SCYL) {
+		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 (srcp->sflags & SDISTANT) {
+		si->dom *= srcp->ss2;
 		return(FHUGE);
-	si->dom /= d*d;
+	}
+	if (si->dom <= 1e-4)
+		goto nextsample;		/* behind source? */
+	si->dom *= srcp->ss2/(d*d);
 	return(d);		/* sample OK, return distance */
 }
 
 
+int
 skipparts(ct, sz, pp, pt)		/* skip to requested partition */
 int  ct[3], sz[3];		/* center and size of partition (returned) */
 register int  pp[2];		/* current index, number to skip (modified) */
@@ -115,12 +142,13 @@ unsigned char  *pt;		/* partition array */
 					/* check this partition */
 	p = spart(pt, pp[0]);
 	pp[0]++;
-	if (p == S0)			/* leaf partition */
+	if (p == S0) {			/* leaf partition */
 		if (pp[1]) {
 			pp[1]--;
 			return(0);	/* not there yet */
 		} else
 			return(1);	/* we've arrived */
+	}
 				/* else check lower */
 	sz[p] >>= 1;
 	ct[p] -= sz[p];
@@ -137,9 +165,10 @@ unsigned char  *pt;		/* partition array */
 }
 
 
-nopart(si, ro)			/* single source partition */
+void
+nopart(si, r)			/* single source partition */
 register SRCINDEX  *si;
-FVECT  ro;
+RAY  *r;
 {
 	clrpart(si->spt);
 	setpart(si->spt, 0, S0);
@@ -147,9 +176,10 @@ FVECT  ro;
 }
 
 
-cylpart(si, ro)			/* partition a cylinder */
+void
+cylpart(si, r)			/* partition a cylinder */
 SRCINDEX  *si;
-FVECT  ro;
+register RAY  *r;
 {
 	double  dist2, safedist2, dist2cent, rad2;
 	FVECT  v;
@@ -157,11 +187,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 +201,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);
 }
 
@@ -223,21 +253,35 @@ double  d2;
 }
 
 
-flatpart(si, ro)			/* partition a flat source */
+void
+flatpart(si, r)				/* partition a flat source */
 register SRCINDEX  *si;
-FVECT  ro;
+register RAY  *r;
 {
-	register double  *vp;
+	register RREAL  *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) <= 0.) {		/* 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);
 }
 
@@ -298,7 +342,7 @@ scylform(sn, dir)		/* compute cosine for cylinder's pr
 int  sn;
 register FVECT  dir;		/* assume normalized */
 {
-	register double  *dv;
+	register RREAL  *dv;
 	double  d;
 
 	dv = source[sn].ss[SU];