--- ray/src/gen/genprism.c	1989/03/14 15:08:35	1.2
+++ ray/src/gen/genprism.c	2003/07/21 22:30:18	2.11
@@ -1,42 +1,311 @@
-/* Copyright (c) 1987 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: genprism.c,v 2.11 2003/07/21 22:30:18 schorsch Exp $";
 #endif
-
 /*
  *  genprism.c - generate a prism.
  *		2D vertices in the xy plane are given on the
  *		command line or from a file.  Their order together
  *		with the extrude direction will determine surface
  *		orientation.
- *
- *     8/24/87
  */
 
 #include  <stdio.h>
-
+#include  <string.h>
+#include <stdlib.h>
+#include  <math.h>
 #include  <ctype.h>
 
 #define  MAXVERT	1024		/* maximum # vertices */
 
+#define  FTINY		1e-6
+
 char  *pmtype;		/* material type */
 char  *pname;		/* name */
 
 double  lvect[3] = {0.0, 0.0, 1.0};
+int	lvdir = 1;
+double	llen = 1.0;
 
 double  vert[MAXVERT][2];
 int  nverts = 0;
 
+double	u[MAXVERT][2];		/* edge unit vectors */
+double	a[MAXVERT];		/* corner trim sizes */
+
 int  do_ends = 1;		/* include end caps */
 int  iscomplete = 0;		/* polygon is already completed */
+double  crad = 0.0;		/* radius for corners (sign from lvdir) */
 
+static double  compute_rounding(void);
 
+
+static void
+readverts(fname)		/* read vertices from a file */
+char  *fname;
+{
+	FILE  *fp;
+
+	if (fname == NULL)
+		fp = stdin;
+	else if ((fp = fopen(fname, "r")) == NULL) {
+		fprintf(stderr, "%s: cannot open\n", fname);
+		exit(1);
+	}
+	while (fscanf(fp, "%lf %lf", &vert[nverts][0], &vert[nverts][1]) == 2)
+		nverts++;
+	fclose(fp);
+}
+
+
+static void
+side(n0, n1)			/* print single side */
+register int  n0, n1;
+{
+	printf("\n%s polygon %s.%d\n", pmtype, pname, n0+1);
+	printf("0\n0\n12\n");
+	printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n0][0],
+			vert[n0][1], 0.0);
+	printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n0][0]+lvect[0],
+			vert[n0][1]+lvect[1], lvect[2]);
+	printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n1][0]+lvect[0],
+			vert[n1][1]+lvect[1], lvect[2]);
+	printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n1][0],
+			vert[n1][1], 0.0);
+}
+
+
+static void
+rside(n0, n1)			/* print side with rounded edge */
+register int  n0, n1;
+{
+	double  s, c, t[3];
+
+					/* compute tanget offset vector */
+	s = lvdir*(lvect[1]*u[n1][0] - lvect[0]*u[n1][1])/llen;
+	if (s < -FTINY || s > FTINY) {
+		c = sqrt(1. - s*s);
+		t[0] = (c - 1.)*u[n1][1];
+		t[1] = (1. - c)*u[n1][0];
+		t[2] = s;
+	} else
+		t[0] = t[1] = t[2] = 0.;
+					/* output side */
+	printf("\n%s polygon %s.%d\n", pmtype, pname, n0+1);
+	printf("0\n0\n12\n");
+	printf("\t%18.12g\t%18.12g\t%18.12g\n",
+			vert[n0][0] + crad*(t[0] - a[n0]*u[n1][0]),
+			vert[n0][1] + crad*(t[1] - a[n0]*u[n1][1]),
+			crad*(t[2] + 1.));
+	printf("\t%18.12g\t%18.12g\t%18.12g\n",
+			vert[n0][0] + lvect[0] + crad*(t[0] - a[n0]*u[n1][0]),
+			vert[n0][1] + lvect[1] + crad*(t[1] - a[n0]*u[n1][1]),
+			lvect[2] + crad*(t[2] - 1.));
+	printf("\t%18.12g\t%18.12g\t%18.12g\n",
+			vert[n1][0] + lvect[0] + crad*(t[0] + a[n1]*u[n1][0]),
+			vert[n1][1] + lvect[1] + crad*(t[1] + a[n1]*u[n1][1]),
+			lvect[2] + crad*(t[2] - 1.));
+	printf("\t%18.12g\t%18.12g\t%18.12g\n",
+			vert[n1][0] + crad*(t[0] + a[n1]*u[n1][0]),
+			vert[n1][1] + crad*(t[1] + a[n1]*u[n1][1]),
+			crad*(t[2] + 1.));
+					/* output joining edge */
+	if (lvdir*a[n1] < 0.)
+		return;
+	printf("\n%s cylinder %s.e%d\n", pmtype, pname, n0+1);
+	printf("0\n0\n7\n");
+	printf("\t%18.12g\t%18.12g\t%18.12g\n",
+		vert[n1][0] + crad*(a[n1]*u[n1][0] - u[n1][1]),
+		vert[n1][1] + crad*(a[n1]*u[n1][1] + u[n1][0]),
+		crad);
+	printf("\t%18.12g\t%18.12g\t%18.12g\n",
+		vert[n1][0] + lvect[0] + crad*(a[n1]*u[n1][0] - u[n1][1]),
+		vert[n1][1] + lvect[1] + crad*(a[n1]*u[n1][1] + u[n1][0]),
+		lvect[2] - crad);
+	printf("\t%18.12g\n", lvdir*crad);
+}
+
+
+static double
+compute_rounding()		/* compute vectors for rounding operations */
+{
+	register int  i;
+	register double	*v0, *v1;
+	double  l, asum;
+
+	v0 = vert[nverts-1];
+	for (i = 0; i < nverts; i++) {		/* compute u[*] */
+		v1 = vert[i];
+		u[i][0] = v0[0] - v1[0];
+		u[i][1] = v0[1] - v1[1];
+		l = sqrt(u[i][0]*u[i][0] + u[i][1]*u[i][1]);
+		if (l <= FTINY) {
+			fprintf(stderr, "Degenerate side in prism\n");
+			exit(1);
+		}
+		u[i][0] /= l;
+		u[i][1] /= l;
+		v0 = v1;
+	}
+	asum = 0.;
+	v1 = u[0];
+	for (i = nverts; i--; ) {		/* compute a[*] */
+		v0 = u[i];
+		l = v0[0]*v1[0] + v0[1]*v1[1];
+		if (1+l <= FTINY) {
+			fprintf(stderr, "Overlapping sides in prism\n");
+			exit(1);
+		}
+		if (1-l <= 0.)
+			a[i] = 0.;
+		else {
+			a[i] = sqrt((1-l)/(1+l));
+			asum += l = v1[0]*v0[1]-v1[1]*v0[0];
+			if (l < 0.)
+				a[i] = -a[i];
+		}
+		v1 = v0;
+	}
+	return(asum*.5);
+}
+
+
+static void
+printends()			/* print ends of prism */
+{
+	register int  i;
+						/* bottom face */
+	printf("\n%s polygon %s.b\n", pmtype, pname);
+	printf("0\n0\n%d\n", nverts*3);
+	for (i = 0; i < nverts; i++)
+		printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[i][0],
+				vert[i][1], 0.0);
+						/* top face */
+	printf("\n%s polygon %s.t\n", pmtype, pname);
+	printf("0\n0\n%d\n", nverts*3);
+	for (i = nverts; i--; )
+		printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[i][0]+lvect[0],
+				vert[i][1]+lvect[1], lvect[2]);
+}
+
+
+static void
+printrends()			/* print ends of prism with rounding */
+{
+	register int  i;
+	double	c0[3], c1[3], cl[3];
+						/* bottom face */
+	printf("\n%s polygon %s.b\n", pmtype, pname);
+	printf("0\n0\n%d\n", nverts*3);
+	for (i = 0; i < nverts; i++) {
+		printf("\t%18.12g\t%18.12g\t%18.12g\n",
+			vert[i][0] + crad*(a[i]*u[i][0] - u[i][1]),
+			vert[i][1] + crad*(a[i]*u[i][1] + u[i][0]),
+			0.0);
+	}
+						/* top face */
+	printf("\n%s polygon %s.t\n", pmtype, pname);
+	printf("0\n0\n%d\n", nverts*3);
+	for (i = nverts; i--; ) {
+		printf("\t%18.12g\t%18.12g\t%18.12g\n",
+		vert[i][0] + lvect[0] + crad*(a[i]*u[i][0] - u[i][1]),
+		vert[i][1] + lvect[1] + crad*(a[i]*u[i][1] + u[i][0]),
+		lvect[2]);
+	}
+						/* bottom corners and edges */
+	c0[0] = cl[0] = vert[nverts-1][0] +
+			crad*(a[nverts-1]*u[nverts-1][0] - u[nverts-1][1]);
+	c0[1] = cl[1] = vert[nverts-1][1] +
+			crad*(a[nverts-1]*u[nverts-1][1] + u[nverts-1][0]);
+	c0[2] = cl[2] = crad;
+	for (i = 0; i < nverts; i++) {
+		if (i < nverts-1) {
+			c1[0] = vert[i][0] + crad*(a[i]*u[i][0] - u[i][1]);
+			c1[1] = vert[i][1] + crad*(a[i]*u[i][1] + u[i][0]);
+			c1[2] = crad;
+		} else {
+			c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
+		}
+		if (lvdir*a[i] > 0.) {
+			printf("\n%s sphere %s.bc%d\n", pmtype, pname, i+1);
+			printf("0\n0\n4 %18.12g %18.12g %18.12g %18.12g\n",
+				c1[0], c1[1], c1[2], lvdir*crad);
+		}
+		printf("\n%s cylinder %s.be%d\n", pmtype, pname, i+1);
+		printf("0\n0\n7\n");
+		printf("\t%18.12g\t%18.12g\t%18.12g\n", c0[0], c0[1], c0[2]);
+		printf("\t%18.12g\t%18.12g\t%18.12g\n", c1[0], c1[1], c1[2]);
+		printf("\t%18.12g\n", lvdir*crad);
+		c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
+	}
+						/* top corners and edges */
+	c0[0] = cl[0] = vert[nverts-1][0] + lvect[0] +
+			crad*(a[nverts-1]*u[nverts-1][0] - u[nverts-1][1]);
+	c0[1] = cl[1] = vert[nverts-1][1] + lvect[1] +
+			crad*(a[nverts-1]*u[nverts-1][1] + u[nverts-1][0]);
+	c0[2] = cl[2] = lvect[2] - crad;
+	for (i = 0; i < nverts; i++) {
+		if (i < nverts-1) {
+			c1[0] = vert[i][0] + lvect[0] +
+					crad*(a[i]*u[i][0] - u[i][1]);
+			c1[1] = vert[i][1] + lvect[1] +
+					crad*(a[i]*u[i][1] + u[i][0]);
+			c1[2] = lvect[2] - crad;
+		} else {
+			c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
+		}
+		if (lvdir*a[i] > 0.) {
+			printf("\n%s sphere %s.tc%d\n", pmtype, pname, i+1);
+			printf("0\n0\n4 %18.12g %18.12g %18.12g %18.12g\n",
+				c1[0], c1[1], c1[2], lvdir*crad);
+		}
+		printf("\n%s cylinder %s.te%d\n", pmtype, pname, i+1);
+		printf("0\n0\n7\n");
+		printf("\t%18.12g\t%18.12g\t%18.12g\n", c0[0], c0[1], c0[2]);
+		printf("\t%18.12g\t%18.12g\t%18.12g\n", c1[0], c1[1], c1[2]);
+		printf("\t%18.12g\n", lvdir*crad);
+		c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
+	}
+}
+
+
+static void
+printsides(round)		/* print prism sides */
+int  round;
+{
+	register int  i;
+
+	for (i = 0; i < nverts-1; i++)
+		if (round)
+			rside(i, i+1);
+		else
+			side(i, i+1);
+	if (!iscomplete) {
+		if (round)
+			rside(nverts-1, 0);
+		else
+			side(nverts-1, 0);
+	}
+}
+
+
+static void
+printhead(ac, av)		/* print command header */
+register int  ac;
+register char  **av;
+{
+	putchar('#');
+	while (ac--) {
+		putchar(' ');
+		fputs(*av++, stdout);
+	}
+	putchar('\n');
+}
+
+
 main(argc, argv)
 int  argc;
 char  **argv;
 {
-	double  atof();
 	int  an;
 	
 	if (argc < 4)
@@ -74,7 +343,22 @@ char  **argv;
 			lvect[0] = atof(argv[++an]);
 			lvect[1] = atof(argv[++an]);
 			lvect[2] = atof(argv[++an]);
+			if (lvect[2] < -FTINY)
+				lvdir = -1;
+			else if (lvect[2] > FTINY)
+				lvdir = 1;
+			else {
+				fprintf(stderr,
+					"%s: illegal extrusion vector\n",
+						argv[0]);
+				exit(1);
+			}
+			llen = sqrt(lvect[0]*lvect[0] + lvect[1]*lvect[1] +
+					lvect[2]*lvect[2]);
 			break;
+		case 'r':				/* radius */
+			crad = atof(argv[++an]);
+			break;
 		case 'e':				/* ends */
 			do_ends = !do_ends;
 			break;
@@ -85,105 +369,29 @@ char  **argv;
 			goto userr;
 		}
 	}
-
-	printhead(argc, argv);
-
-	if (do_ends)
-		printends();
-	printsides();
-
-	return(0);
+	if (crad > FTINY) {
+		if (crad > lvdir*lvect[2]) {
+			fprintf(stderr, "%s: rounding greater than height\n",
+					argv[0]);
+			exit(1);
+		}
+		crad *= lvdir;		/* simplifies formulas */
+		compute_rounding();
+		printhead(argc, argv);
+		if (do_ends)
+			printrends();
+		printsides(1);
+	} else {
+		printhead(argc, argv);
+		if (do_ends)
+			printends();
+		printsides(0);
+	}
+	exit(0);
 userr:
 	fprintf(stderr, "Usage: %s material name ", argv[0]);
 	fprintf(stderr, "{ - | vfile | N v1 v2 .. vN } ");
-	fprintf(stderr, "[-l lvect][-c][-e]\n");
+	fprintf(stderr, "[-l lvect][-r radius][-c][-e]\n");
 	exit(1);
 }
 
-
-readverts(fname)		/* read vertices from a file */
-char  *fname;
-{
-	FILE  *fp;
-
-	if (fname == NULL)
-		fp = stdin;
-	else if ((fp = fopen(fname, "r")) == NULL) {
-		fprintf(stderr, "%s: cannot open\n", fname);
-		exit(1);
-	}
-	while (fscanf(fp, "%lf %lf", &vert[nverts][0], &vert[nverts][1]) == 2)
-		nverts++;
-	fclose(fp);
-}
-
-
-printends()			/* print ends of prism */
-{
-	register int  i;
-
-	printf("\n%s polygon %s.b\n", pmtype, pname);
-	printf("0\n0\n%d\n", nverts*3);
-	for (i = 0; i < nverts; i++) {
-		printf("\t%18.12g\t%18.12g\t%18.12g\n",
-				vert[i][0],
-				vert[i][1],
-				0.0);
-	}
-	printf("\n%s polygon %s.e\n", pmtype, pname);
-	printf("0\n0\n%d\n", nverts*3);
-	for (i = nverts-1; i >= 0; i--) {
-		printf("\t%18.12g\t%18.12g\t%18.12g\n",
-				vert[i][0]+lvect[0],
-				vert[i][1]+lvect[1],
-				lvect[2]);
-	}
-}
-
-
-printsides()			/* print prism sides */
-{
-	register int  i;
-
-	for (i = 0; i < nverts-1; i++)
-		side(i, i+1);
-	if (!iscomplete)
-		side(nverts-1, 0);
-}
-
-
-side(n1, n2)			/* print single side */
-register int  n1, n2;
-{
-	printf("\n%s polygon %s.%d\n", pmtype, pname, n1+1);
-	printf("0\n0\n12\n");
-	printf("\t%18.12g\t%18.12g\t%18.12g\n",
-			vert[n1][0],
-			vert[n1][1],
-			0.0);
-	printf("\t%18.12g\t%18.12g\t%18.12g\n",
-			vert[n1][0]+lvect[0],
-			vert[n1][1]+lvect[1],
-			lvect[2]);
-	printf("\t%18.12g\t%18.12g\t%18.12g\n",
-			vert[n2][0]+lvect[0],
-			vert[n2][1]+lvect[1],
-			lvect[2]);
-	printf("\t%18.12g\t%18.12g\t%18.12g\n",
-			vert[n2][0],
-			vert[n2][1],
-			0.0);
-}
-
-
-printhead(ac, av)		/* print command header */
-register int  ac;
-register char  **av;
-{
-	putchar('#');
-	while (ac--) {
-		putchar(' ');
-		fputs(*av++, stdout);
-	}
-	putchar('\n');
-}