--- ray/src/rt/srcsupp.c	2007/02/26 21:16:02	2.16
+++ ray/src/rt/srcsupp.c	2024/08/07 18:33:33	2.25
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: srcsupp.c,v 2.16 2007/02/26 21:16:02 greg Exp $";
+static const char	RCSid[] = "$Id: srcsupp.c,v 2.25 2024/08/07 18:33:33 greg Exp $";
 #endif
 /*
  *  Support routines for source objects and materials
@@ -19,7 +19,7 @@ static const char	RCSid[] = "$Id: srcsupp.c,v 2.16 200
 
 #include  "face.h"
 
-#define SRCINC		8		/* realloc increment for array */
+#define SRCINC		32		/* realloc increment for array */
 
 SRCREC  *source = NULL;			/* our list of sources */
 int  nsources = 0;			/* the number of sources */
@@ -28,7 +28,7 @@ SRCFUNC  sfun[NUMOTYPE];		/* source dispatch table */
 
 
 void
-initstypes()			/* initialize source dispatch table */
+initstypes(void)			/* initialize source dispatch table */
 {
 	extern VSMATERIAL  mirror_vs, direct1_vs, direct2_vs;
 	static SOBJECT  fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
@@ -49,7 +49,7 @@ initstypes()			/* initialize source dispatch table */
 
 
 int
-newsource()			/* allocate new source in our array */
+newsource(void)			/* allocate new source in our array */
 {
 	if (nsources == 0)
 		source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
@@ -69,39 +69,36 @@ newsource()			/* allocate new source in our array */
 
 
 void
-setflatss(src)				/* set sampling for a flat source */
-register SRCREC  *src;
+setflatss(				/* set sampling for a flat source */
+	SRCREC  *src
+)
 {
 	double  mult;
-	register int  i;
+	int  i;
 
-	src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
+	getperpendicular(src->ss[SU], src->snorm, rand_samp);
+	mult = .5 * sqrt( src->ss2 );
 	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]);
 }
 
 
 void
-fsetsrc(src, so)			/* set a face as a source */
-register SRCREC  *src;
-OBJREC  *so;
+fsetsrc(			/* set a face as a source */
+	SRCREC  *src,
+	OBJREC  *so
+)
 {
-	register FACE  *f;
-	register int  i, j;
+	FACE  *f;
+	int  i, j;
 	double  d;
 	
 	src->sa.success = 2*AIMREQT-1;		/* bitch on second failure */
 	src->so = so;
 						/* get the face */
 	f = getface(so);
-	if (f->area == 0.0)
+	if (f->area == 0.)
 		objerror(so, USER, "zero source area");
 						/* find the center */
 	for (j = 0; j < 3; j++) {
@@ -124,20 +121,41 @@ OBJREC  *so;
 	}
 	src->srad = sqrt(src->srad);
 						/* compute size vectors */
-	if (f->nv == 4)				/* parallelogram case */
+	if (f->nv == 4) {				/* parallelogram case */
 		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);
+	} else if (f->nv == 3) {			/* triangle case */
+		int	near0 = 2;
+		double	dmin = dist2line(src->sloc, VERTEX(f,2), VERTEX(f,0));
+		for (i = 0; i < 2; i++) {
+			double	d2 = dist2line(src->sloc, VERTEX(f,i), VERTEX(f,i+1));
+			if (d2 >= dmin)
+				continue;
+			near0 = i;
+			dmin = d2;			/* radius = min distance */
+		}
+		if (dmin < .08*f->area)
+			objerror(so, WARNING, "triangular source with poor aspect");
+		i = (near0 + 1) % 3;
+		for (j = 0; j < 3; j++)
+			src->ss[SU][j] = VERTEX(f,i)[j] - VERTEX(f,near0)[j];
+		normalize(src->ss[SU]);
+		dmin = sqrt(dmin);
+		for (j = 0; j < 3; j++)
+			src->ss[SU][j] *= dmin;
+		fcross(src->ss[SV], f->norm, src->ss[SU]);
+	} else
+		setflatss(src);				/* hope for convex! */
 }
 
 
 void
-ssetsrc(src, so)			/* set a source as a source */
-register SRCREC  *src;
-register OBJREC  *so;
+ssetsrc(			/* set a source as a source */
+	SRCREC  *src,
+	OBJREC  *so
+)
 {
 	double  theta;
 	
@@ -145,7 +163,7 @@ register OBJREC  *so;
 	src->so = so;
 	if (so->oargs.nfargs != 4)
 		objerror(so, USER, "bad arguments");
-	src->sflags |= SDISTANT;
+	src->sflags |= (SDISTANT|SCIR);
 	VCOPY(src->sloc, so->oargs.farg);
 	if (normalize(src->sloc) == 0.0)
 		objerror(so, USER, "zero direction");
@@ -162,11 +180,12 @@ register OBJREC  *so;
 
 
 void
-sphsetsrc(src, so)			/* set a sphere as a source */
-register SRCREC  *src;
-register OBJREC  *so;
+sphsetsrc(			/* set a sphere as a source */
+	SRCREC  *src,
+	OBJREC  *so
+)
 {
-	register int  i;
+	int  i;
 
 	src->sa.success = 2*AIMREQT-1;		/* bitch on second failure */
 	src->so = so;
@@ -174,33 +193,36 @@ register OBJREC  *so;
 		objerror(so, USER, "bad # arguments");
 	if (so->oargs.farg[3] <= FTINY)
 		objerror(so, USER, "illegal source radius");
+	src->sflags |= SCIR;
 	VCOPY(src->sloc, so->oargs.farg);
 	src->srad = so->oargs.farg[3];
 	src->ss2 = PI * src->srad * src->srad;
+	memset(src->ss, 0, sizeof(src->ss));
 	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];
+		src->ss[i][i] = 0.7236 * so->oargs.farg[3];
 }
 
 
 void
-rsetsrc(src, so)			/* set a ring (disk) as a source */
-register SRCREC  *src;
-OBJREC  *so;
+rsetsrc(			/* set a ring (disk) as a source */
+	SRCREC  *src,
+	OBJREC  *so
+)
 {
-	register CONE  *co;
+	CONE  *co;
 	
 	src->sa.success = 2*AIMREQT-1;		/* bitch on second failure */
 	src->so = so;
 						/* get the ring */
 	co = getcone(so, 0);
+	if (co == NULL)
+		objerror(so, USER, "illegal source");
 	if (CO_R1(co) <= FTINY)
 		objerror(so, USER, "illegal source radius");
 	VCOPY(src->sloc, CO_P0(co));
 	if (CO_R0(co) > 0.0)
 		objerror(so, USER, "cannot hit source center");
-	src->sflags |= SFLAT;
+	src->sflags |= (SFLAT|SCIR);
 	VCOPY(src->snorm, co->ad);
 	src->srad = CO_R1(co);
 	src->ss2 = PI * src->srad * src->srad;
@@ -209,17 +231,20 @@ OBJREC  *so;
 
 
 void
-cylsetsrc(src, so)			/* set a cylinder as a source */
-register SRCREC  *src;
-OBJREC  *so;
+cylsetsrc(			/* set a cylinder as a source */
+	SRCREC  *src,
+	OBJREC  *so
+)
 {
-	register CONE  *co;
-	register int  i;
+	CONE  *co;
+	int  i;
 	
 	src->sa.success = 4*AIMREQT-1;		/* bitch on fourth failure */
 	src->so = so;
 						/* get the cylinder */
 	co = getcone(so, 0);
+	if (co == NULL)
+		objerror(so, USER, "illegal source");
 	if (CO_R0(co) <= FTINY)
 		objerror(so, USER, "illegal source radius");
 	if (CO_R0(co) > .2*co->al)		/* heuristic constraint */
@@ -232,24 +257,19 @@ OBJREC  *so;
 						/* 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;
+	getperpendicular(src->ss[SW], co->ad, rand_samp);
 	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;
+makespot(			/* make a spotlight */
+	OBJREC  *m
+)
 {
-	register SPOT  *ns;
+	SPOT  *ns;
 
 	if ((ns = (SPOT *)m->os) != NULL)
 		return(ns);
@@ -267,9 +287,10 @@ register OBJREC  *m;
 
 
 int
-spotout(r, s)			/* check if we're outside spot region */
-register RAY  *r;
-register SPOT  *s;
+spotout(			/* check if we're outside spot region */
+	RAY  *r,
+	SPOT  *s
+)
 {
 	double  d;
 	FVECT  vd;
@@ -297,14 +318,15 @@ register SPOT  *s;
 
 
 double
-fgetmaxdisk(ocent, op)		/* get center and squared radius of face */
-FVECT  ocent;
-OBJREC  *op;
+fgetmaxdisk(		/* get center and squared radius of face */
+	FVECT  ocent,
+	OBJREC  *op
+)
 {
 	double  maxrad2;
 	double  d;
-	register int  i, j;
-	register FACE  *f;
+	int  i, j;
+	FACE  *f;
 	
 	f = getface(op);
 	if (f->area == 0.)
@@ -329,24 +351,28 @@ OBJREC  *op;
 
 
 double
-rgetmaxdisk(ocent, op)		/* get center and squared radius of ring */
-FVECT  ocent;
-OBJREC  *op;
+rgetmaxdisk(		/* get center and squared radius of ring */
+	FVECT  ocent,
+	OBJREC  *op
+)
 {
-	register CONE  *co;
+	CONE  *co;
 	
 	co = getcone(op, 0);
+	if (co == NULL)
+		return(0.);
 	VCOPY(ocent, CO_P0(co));
 	return(CO_R1(co)*CO_R1(co));
 }
 
 
 double
-fgetplaneq(nvec, op)			/* get plane equation for face */
-FVECT  nvec;
-OBJREC  *op;
+fgetplaneq(			/* get plane equation for face */
+	FVECT  nvec,
+	OBJREC  *op
+)
 {
-	register FACE  *fo;
+	FACE  *fo;
 
 	fo = getface(op);
 	VCOPY(nvec, fo->norm);
@@ -355,22 +381,29 @@ OBJREC  *op;
 
 
 double
-rgetplaneq(nvec, op)			/* get plane equation for ring */
-FVECT  nvec;
-OBJREC  *op;
+rgetplaneq(			/* get plane equation for ring */
+	FVECT  nvec,
+	OBJREC  *op
+)
 {
-	register CONE  *co;
+	CONE  *co;
 
 	co = getcone(op, 0);
+	if (co == NULL) {
+		memset(nvec, 0, sizeof(FVECT));
+		return(0.);
+	}
 	VCOPY(nvec, co->ad);
 	return(DOT(nvec, CO_P0(co)));
 }
 
 
 int
-commonspot(sp1, sp2, org)	/* set sp1 to intersection of sp1 and sp2 */
-register SPOT  *sp1, *sp2;
-FVECT  org;
+commonspot(		/* set sp1 to intersection of sp1 and sp2 */
+	SPOT  *sp1,
+	SPOT  *sp2,
+	FVECT  org
+)
 {
 	FVECT  cent;
 	double  rad2, cos1, cos2;
@@ -392,20 +425,19 @@ FVECT  org;
 
 
 int
-commonbeam(sp1, sp2, dir)	/* set sp1 to intersection of sp1 and sp2 */
-register SPOT  *sp1, *sp2;
-FVECT  dir;
+commonbeam(		/* set sp1 to intersection of sp1 and sp2 */
+	SPOT  *sp1,
+	SPOT  *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];
+	VSUM(c1, sp1->aim, dir, -d);
 	d = DOT(sp2->aim, dir);
-	for (i = 0; i < 3; i++)
-		c2[i] = sp2->aim[i] - d*dir[i];
+	VSUM(c2, sp2->aim, dir, -d);
 					/* compute overlap */
 	rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
 	if (rad2 <= FTINY)
@@ -417,9 +449,10 @@ FVECT  dir;
 
 
 int
-checkspot(sp, nrm)		/* check spotlight for behind source */
-register SPOT  *sp;	/* spotlight */
-FVECT  nrm;		/* source surface normal */
+checkspot(			/* check spotlight for behind source */
+	SPOT  *sp,	/* spotlight */
+	FVECT  nrm	/* source surface normal */
+)
 {
 	double  d, d1;
 
@@ -433,15 +466,15 @@ FVECT  nrm;		/* source surface normal */
 
 
 double
-spotdisk(oc, op, sp, pos)	/* intersect spot with object op */
-FVECT  oc;
-OBJREC  *op;
-register SPOT  *sp;
-FVECT  pos;
+spotdisk(		/* intersect spot with object op */
+	FVECT  oc,
+	OBJREC  *op,
+	SPOT  *sp,
+	FVECT  pos
+)
 {
 	FVECT  onorm;
 	double  offs, d, dist;
-	register int  i;
 
 	offs = getplaneq(onorm, op);
 	d = -DOT(onorm, sp->aim);
@@ -450,46 +483,45 @@ FVECT  pos;
 	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];
+	VSUM(oc, pos, sp->aim, dist);
 	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;
+beamdisk(		/* intersect beam with object op */
+	FVECT  oc,
+	OBJREC  *op,
+	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];
+	VSUM(oc, sp->aim, dir, dist);
 	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 */
+intercircle(				/* intersect two circles */
+	FVECT  cc,		/* midpoint (return value) */
+	FVECT  c1,		/* circle centers */
+	FVECT  c2,
+	double  r1s,		/* radii squared */
+	double	r2s
+)
 {
 	double  a2, d2, l;
 	FVECT  disp;
-	register int  i;
 
-	for (i = 0; i < 3; i++)
-		disp[i] = c2[i] - c1[i];
+	VSUB(disp, c2, c1);
 	d2 = DOT(disp,disp);
 					/* circle within overlap? */
 	if (r1s < r2s) {
@@ -509,7 +541,6 @@ double  r1s, r2s;		/* radii squared */
 		return(0.);
 					/* overlap, compute center */
 	l = sqrt((r1s - a2)/d2);
-	for (i = 0; i < 3; i++)
-		cc[i] = c1[i] + l*disp[i];
+	VSUM(cc, c1, disp, l);
 	return(a2);
 }