--- ray/src/rt/srcsupp.c	1991/06/21 13:23:05	1.3
+++ ray/src/rt/srcsupp.c	1991/11/08 16:34:24	1.16
@@ -18,6 +18,7 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #include  "face.h"
 
+#define SRCINC		4		/* realloc increment for array */
 
 SRCREC  *source = NULL;			/* our list of sources */
 int  nsources = 0;			/* the number of sources */
@@ -27,19 +28,24 @@ SRCFUNC  sfun[NUMOTYPE];		/* source dispatch table */
 
 initstypes()			/* initialize source dispatch table */
 {
-	extern VSMATERIAL  mirror_vs;
-	extern int  fsetsrc(), ssetsrc(), sphsetsrc(), rsetsrc();
+	extern VSMATERIAL  mirror_vs, direct1_vs, direct2_vs;
+	extern int  fsetsrc(), ssetsrc(), sphsetsrc(), cylsetsrc(), rsetsrc();
+	extern int  nopart(), flatpart(), cylpart();
 	extern double  fgetplaneq(), rgetplaneq();
 	extern double  fgetmaxdisk(), rgetmaxdisk();
-	static SOBJECT  fsobj = {fsetsrc, fgetplaneq, fgetmaxdisk};
-	static SOBJECT  ssobj = {ssetsrc};
-	static SOBJECT  sphsobj = {sphsetsrc};
-	static SOBJECT  rsobj = {rsetsrc, rgetplaneq, rgetmaxdisk};
+	static SOBJECT  fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
+	static SOBJECT  ssobj = {ssetsrc, nopart};
+	static SOBJECT  sphsobj = {sphsetsrc, nopart};
+	static SOBJECT  cylsobj = {cylsetsrc, cylpart};
+	static SOBJECT  rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
 
 	sfun[MAT_MIRROR].mf = &mirror_vs;
+	sfun[MAT_DIRECT1].mf = &direct1_vs;
+	sfun[MAT_DIRECT2].mf = &direct2_vs;
 	sfun[OBJ_FACE].of = &fsobj;
 	sfun[OBJ_SOURCE].of = &ssobj;
 	sfun[OBJ_SPHERE].of = &sphsobj;
+	sfun[OBJ_CYLINDER].of = &cylsobj;
 	sfun[OBJ_RING].of = &rsobj;
 }
 
@@ -48,10 +54,10 @@ int
 newsource()			/* allocate new source in our array */
 {
 	if (nsources == 0)
-		source = (SRCREC *)malloc(sizeof(SRCREC));
-	else
+		source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
+	else if (nsources%SRCINC == 0)
 		source = (SRCREC *)realloc((char *)source,
-				(unsigned)(nsources+1)*sizeof(SRCREC));
+				(unsigned)(nsources+SRCINC)*sizeof(SRCREC));
 	if (source == NULL)
 		return(-1);
 	source[nsources].sflags = 0;
@@ -61,12 +67,32 @@ newsource()			/* allocate new source in our array */
 }
 
 
+setflatss(src)				/* set sampling for a flat source */
+register SRCREC  *src;
+{
+	double  mult;
+	register int  i;
+
+	src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
+	for (i = 0; i < 3; i++)
+		if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
+			break;
+	src->ss[SV][i] = 1.0;
+	fcross(src->ss[SU], src->ss[SV], src->snorm);
+	mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
+	for (i = 0; i < 3; i++)
+		src->ss[SU][i] *= mult;
+	fcross(src->ss[SV], src->snorm, src->ss[SU]);
+}
+
+
 fsetsrc(src, so)			/* set a face as a source */
 register SRCREC  *src;
 OBJREC  *so;
 {
 	register FACE  *f;
 	register int  i, j;
+	double  d;
 	
 	src->sa.success = 2*AIMREQT-1;		/* bitch on second failure */
 	src->so = so;
@@ -83,8 +109,24 @@ OBJREC  *so;
 		objerror(so, USER, "cannot hit center");
 	src->sflags |= SFLAT;
 	VCOPY(src->snorm, f->norm);
-	src->ss = sqrt(f->area / PI);
 	src->ss2 = f->area;
+						/* find maximum radius */
+	src->srad = 0.;
+	for (i = 0; i < f->nv; i++) {
+		d = dist2(VERTEX(f,i), src->sloc);
+		if (d > src->srad)
+			src->srad = d;
+	}
+	src->srad = sqrt(src->srad);
+						/* compute size vectors */
+	if (f->nv == 4 || (f->nv == 5 &&	/* parallelogram case */
+			dist2(VERTEX(f,0),VERTEX(f,4)) <= FTINY*FTINY))
+		for (j = 0; j < 3; j++) {
+			src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
+			src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
+		}
+	else
+		setflatss(src);
 }
 
 
@@ -105,8 +147,12 @@ register OBJREC  *so;
 	theta = PI/180.0/2.0 * so->oargs.farg[3];
 	if (theta <= FTINY)
 		objerror(so, USER, "zero size");
-	src->ss = theta >= PI/4 ? 1.0 : tan(theta);
 	src->ss2 = 2.0*PI * (1.0 - cos(theta));
+					/* the following is approximate */
+	src->srad = sqrt(src->ss2/PI);
+	VCOPY(src->snorm, src->sloc);
+	setflatss(src);			/* hey, whatever works */
+	src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
 }
 
 
@@ -114,6 +160,8 @@ sphsetsrc(src, so)			/* set a sphere as a source */
 register SRCREC  *src;
 register OBJREC  *so;
 {
+	register int  i;
+
 	src->sa.success = 2*AIMREQT-1;		/* bitch on second failure */
 	src->so = so;
 	if (so->oargs.nfargs != 4)
@@ -121,8 +169,12 @@ register OBJREC  *so;
 	if (so->oargs.farg[3] <= FTINY)
 		objerror(so, USER, "illegal radius");
 	VCOPY(src->sloc, so->oargs.farg);
-	src->ss = so->oargs.farg[3];
-	src->ss2 = PI * src->ss * src->ss;
+	src->srad = so->oargs.farg[3];
+	src->ss2 = PI * src->srad * src->srad;
+	for (i = 0; i < 3; i++)
+		src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
+	for (i = 0; i < 3; i++)
+		src->ss[i][i] = .7236 * so->oargs.farg[3];
 }
 
 
@@ -141,16 +193,50 @@ OBJREC  *so;
 		objerror(so, USER, "cannot hit center");
 	src->sflags |= SFLAT;
 	VCOPY(src->snorm, co->ad);
-	src->ss = CO_R1(co);
-	src->ss2 = PI * src->ss * src->ss;
+	src->srad = CO_R1(co);
+	src->ss2 = PI * src->srad * src->srad;
+	setflatss(src);
 }
 
 
+cylsetsrc(src, so)			/* set a cylinder as a source */
+register SRCREC  *src;
+OBJREC  *so;
+{
+	register CONE  *co;
+	register int  i;
+	
+	src->sa.success = 4*AIMREQT-1;		/* bitch on fourth failure */
+	src->so = so;
+						/* get the cylinder */
+	co = getcone(so, 0);
+	if (CO_R0(co) > .2*co->al)		/* heuristic constraint */
+		objerror(so, WARNING, "source aspect too small");
+	src->sflags |= SCYL;
+	for (i = 0; i < 3; i++)
+		src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
+	src->srad = .5*co->al;
+	src->ss2 = 2.*CO_R0(co)*co->al;
+						/* set sampling vectors */
+	for (i = 0; i < 3; i++)
+		src->ss[SU][i] = .5 * co->al * co->ad[i];
+	src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
+	for (i = 0; i < 3; i++)
+		if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
+			break;
+	src->ss[SV][i] = 1.0;
+	fcross(src->ss[SW], src->ss[SV], co->ad);
+	normalize(src->ss[SW]);
+	for (i = 0; i < 3; i++)
+		src->ss[SW][i] *= .8559 * CO_R0(co);
+	fcross(src->ss[SV], src->ss[SW], co->ad);
+}
+
+
 SPOT *
 makespot(m)			/* make a spotlight */
 register OBJREC  *m;
 {
-	extern double  cos();
 	register SPOT  *ns;
 
 	if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
@@ -169,24 +255,27 @@ FVECT  ocent;
 OBJREC  *op;
 {
 	double  maxrad2;
-	double  d2;
+	double  d;
 	register int  i, j;
 	register FACE  *f;
 	
 	f = getface(op);
+	if (f->area == 0.)
+		return(0.);
 	for (i = 0; i < 3; i++) {
 		ocent[i] = 0.;
 		for (j = 0; j < f->nv; j++)
 			ocent[i] += VERTEX(f,j)[i];
 		ocent[i] /= (double)f->nv;
 	}
-	if (f->area == 0.)
-		return(0.);
+	d = DOT(ocent,f->norm);
+	for (i = 0; i < 3; i++)
+		ocent[i] += (f->offset - d)*f->norm[i];
 	maxrad2 = 0.;
 	for (j = 0; j < f->nv; j++) {
-		d2 = dist2(VERTEX(f,j), ocent);
-		if (d2 > maxrad2)
-			maxrad2 = d2;
+		d = dist2(VERTEX(f,j), ocent);
+		if (d > maxrad2)
+			maxrad2 = d;
 	}
 	return(maxrad2);
 }
@@ -231,6 +320,151 @@ OBJREC  *op;
 }
 
 
+commonspot(sp1, sp2, org)	/* set sp1 to intersection of sp1 and sp2 */
+register SPOT  *sp1, *sp2;
+FVECT  org;
+{
+	FVECT  cent;
+	double  rad2, cos1, cos2;
+
+	cos1 = 1. - sp1->siz/(2.*PI);
+	cos2 = 1. - sp2->siz/(2.*PI);
+	if (sp2->siz >= 2.*PI-FTINY)		/* BIG, just check overlap */
+		return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
+					sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
+				/* compute and check disks */
+	rad2 = intercircle(cent, sp1->aim, sp2->aim,
+			1./(cos1*cos1) - 1.,  1./(cos2*cos2) - 1.);
+	if (rad2 <= FTINY || normalize(cent) == 0.)
+		return(0);
+	VCOPY(sp1->aim, cent);
+	sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
+	return(1);
+}
+
+
+commonbeam(sp1, sp2, dir)	/* set sp1 to intersection of sp1 and sp2 */
+register SPOT  *sp1, *sp2;
+FVECT  dir;
+{
+	FVECT  cent, c1, c2;
+	double  rad2, d;
+	register int  i;
+					/* move centers to common plane */
+	d = DOT(sp1->aim, dir);
+	for (i = 0; i < 3; i++)
+		c1[i] = sp1->aim[i] - d*dir[i];
+	d = DOT(sp2->aim, dir);
+	for (i = 0; i < 3; i++)
+		c2[i] = sp2->aim[i] - d*dir[i];
+					/* compute overlap */
+	rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
+	if (rad2 <= FTINY)
+		return(0);
+	VCOPY(sp1->aim, cent);
+	sp1->siz = PI*rad2;
+	return(1);
+}
+
+
+checkspot(sp, nrm)		/* check spotlight for behind source */
+register SPOT  *sp;	/* spotlight */
+FVECT  nrm;		/* source surface normal */
+{
+	double  d, d1;
+
+	d = DOT(sp->aim, nrm);
+	if (d > FTINY)			/* center in front? */
+		return(1);
+					/* else check horizon */
+	d1 = 1. - sp->siz/(2.*PI);
+	return(1.-FTINY-d*d < d1*d1);
+}
+
+
+double
+spotdisk(oc, op, sp, pos)	/* intersect spot with object op */
+FVECT  oc;
+OBJREC  *op;
+register SPOT  *sp;
+FVECT  pos;
+{
+	FVECT  onorm;
+	double  offs, d, dist;
+	register int  i;
+
+	offs = getplaneq(onorm, op);
+	d = -DOT(onorm, sp->aim);
+	if (d >= -FTINY && d <= FTINY)
+		return(0.);
+	dist = (DOT(pos, onorm) - offs)/d;
+	if (dist < 0.)
+		return(0.);
+	for (i = 0; i < 3; i++)
+		oc[i] = pos[i] + dist*sp->aim[i];
+	return(sp->siz*dist*dist/PI/(d*d));
+}
+
+
+double
+beamdisk(oc, op, sp, dir)	/* intersect beam with object op */
+FVECT  oc;
+OBJREC  *op;
+register SPOT  *sp;
+FVECT  dir;
+{
+	FVECT  onorm;
+	double  offs, d, dist;
+	register int  i;
+
+	offs = getplaneq(onorm, op);
+	d = -DOT(onorm, dir);
+	if (d >= -FTINY && d <= FTINY)
+		return(0.);
+	dist = (DOT(sp->aim, onorm) - offs)/d;
+	for (i = 0; i < 3; i++)
+		oc[i] = sp->aim[i] + dist*dir[i];
+	return(sp->siz/PI/(d*d));
+}
+
+
+double
+intercircle(cc, c1, c2, r1s, r2s)	/* intersect two circles */
+FVECT  cc;			/* midpoint (return value) */
+FVECT  c1, c2;			/* circle centers */
+double  r1s, r2s;		/* radii squared */
+{
+	double  a2, d2, l;
+	FVECT  disp;
+	register int  i;
+
+	for (i = 0; i < 3; i++)
+		disp[i] = c2[i] - c1[i];
+	d2 = DOT(disp,disp);
+					/* circle within overlap? */
+	if (r1s < r2s) {
+		if (r2s >= r1s + d2) {
+			VCOPY(cc, c1);
+			return(r1s);
+		}
+	} else {
+		if (r1s >= r2s + d2) {
+			VCOPY(cc, c2);
+			return(r2s);
+		}
+	}
+	a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
+					/* no overlap? */
+	if (a2 <= 0.)
+		return(0.);
+					/* overlap, compute center */
+	l = sqrt((r1s - a2)/d2);
+	for (i = 0; i < 3; i++)
+		cc[i] = c1[i] + l*disp[i];
+	return(a2);
+}
+
+
 sourcehit(r)			/* check to see if ray hit distant source */
 register RAY  *r;
 {
@@ -243,7 +477,7 @@ register RAY  *r;
 		first = 0; last = nsources-1;
 	}
 	for (i = first; i <= last; i++)
-		if (source[i].sflags & SDISTANT)
+		if ((source[i].sflags & (SDISTANT|SVIRTUAL)) == SDISTANT)
 			/*
 			 * Check to see if ray is within
 			 * solid angle of source.
@@ -266,17 +500,71 @@ register RAY  *r;
 }
 
 
-#define  wrongsource(m, r)	(m->otype!=MAT_ILLUM && \
-				r->rsrc>=0 && \
-				source[r->rsrc].so!=r->ro)
+/****************************************************************
+ * The following macros were separated from the m_light() routine
+ * because they are very nasty and difficult to understand.
+ */
 
+/* wrongillum *
+ *
+ * We cannot allow an illum to pass to another illum, because that
+ * would almost certainly constitute overcounting.
+ * However, we do allow an illum to pass to another illum
+ * that is actually going to relay to a virtual light source.
+ */
+
+#define  wrongillum(m, r)	(!(source[r->rsrc].sflags&SVIRTUAL) && \
+			objptr(source[r->rsrc].so->omod)->otype==MAT_ILLUM)
+
+/* wrongsource *
+ *
+ * This source is the wrong source (ie. overcounted) if we are
+ * aimed to a different source than the one we hit and the one
+ * we hit is not an illum which should be passed.
+ */
+
+#define  wrongsource(m, r)	(r->rsrc>=0 && source[r->rsrc].so!=r->ro && \
+				(m->otype!=MAT_ILLUM || wrongillum(m,r)))
+
+/* distglow *
+ *
+ * A distant glow is an object that sometimes acts as a light source,
+ * but is too far away from the test point to be one in this case.
+ */
+
+#define  distglow(m, r)		(m->otype==MAT_GLOW && \
+				r->rot > m->oargs.farg[3])
+
+/* badambient *
+ *
+ * We must avoid including counting light sources in the ambient calculation,
+ * since the direct component is handled separately.  Therefore, any
+ * ambient ray which hits an active light source must be discarded.
+ */
+
 #define  badambient(m, r)	((r->crtype&(AMBIENT|SHADOW))==AMBIENT && \
-				!(m->otype==MAT_GLOW&&r->rot>m->oargs.farg[3]))
+				!distglow(m, r))
 
+/* passillum *
+ *
+ * An illum passes to another material type when we didn't hit it
+ * on purpose (as part of a direct calculation), or it is relaying
+ * a virtual light source.
+ */
+
 #define  passillum(m, r)	(m->otype==MAT_ILLUM && \
-				!(r->rsrc>=0&&source[r->rsrc].so==r->ro))
+				(r->rsrc<0 || source[r->rsrc].so!=r->ro || \
+				source[r->rsrc].sflags&SVIRTUAL))
 
+/* srcignore *
+ *
+ * The -di flag renders light sources invisible, and here is the test.
+ */
 
+#define  srcignore(m, r)	(directinvis && !(r->crtype&SHADOW) && \
+				!distglow(m, r))
+
+
 m_light(m, r)			/* ray hit a light source */
 register OBJREC  *m;
 register RAY  *r;
@@ -286,24 +574,25 @@ register RAY  *r;
 		return;
 						/* check for passed illum */
 	if (passillum(m, r)) {
-
 		if (m->oargs.nsargs < 1 || !strcmp(m->oargs.sarg[0], VOIDID))
 			raytrans(r);
 		else
 			rayshade(r, modifier(m->oargs.sarg[0]));
-
-						/* otherwise treat as source */
-	} else {
+		return;
+	}
+					/* otherwise treat as source */
 						/* check for behind */
-		if (r->rod < 0.0)
-			return;
+	if (r->rod < 0.0)
+		return;
+						/* check for invisibility */
+	if (srcignore(m, r))
+		return;
 						/* get distribution pattern */
-		raytexture(r, m->omod);
+	raytexture(r, m->omod);
 						/* get source color */
-		setcolor(r->rcol, m->oargs.farg[0],
-				  m->oargs.farg[1],
-				  m->oargs.farg[2]);
+	setcolor(r->rcol, m->oargs.farg[0],
+			  m->oargs.farg[1],
+			  m->oargs.farg[2]);
 						/* modify value */
-		multcolor(r->rcol, r->pcol);
-	}
+	multcolor(r->rcol, r->pcol);
 }