--- ray/src/gen/genprism.c	1996/03/19 16:30:56	2.6
+++ ray/src/gen/genprism.c	2003/07/12 09:41:41	2.10
@@ -1,9 +1,6 @@
-/* Copyright (c) 1996 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: genprism.c,v 2.10 2003/07/12 09:41:41 schorsch Exp $";
 #endif
-
 /*
  *  genprism.c - generate a prism.
  *		2D vertices in the xy plane are given on the
@@ -13,23 +10,20 @@ static char SCCSid[] = "$SunId$ LBL";
  */
 
 #include  <stdio.h>
-
+#include  <string.h>
+#include <stdlib.h>
 #include  <math.h>
-
 #include  <ctype.h>
 
 #define  MAXVERT	1024		/* maximum # vertices */
 
 #define  FTINY		1e-6
 
-#ifdef  DCL_ATOF
-extern double  atof();
-#endif
-
 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];
@@ -40,100 +34,12 @@ 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 */
+double  crad = 0.0;		/* radius for corners (sign from lvdir) */
 
-#define rounding	(crad > FTINY)
+static double  compute_rounding(void);
 
 
-main(argc, argv)
-int  argc;
-char  **argv;
-{
-	int  an;
-	
-	if (argc < 4)
-		goto userr;
-
-	pmtype = argv[1];
-	pname = argv[2];
-
-	if (!strcmp(argv[3], "-")) {
-		readverts(NULL);
-		an = 4;
-	} else if (isdigit(argv[3][0])) {
-		nverts = atoi(argv[3]);
-		if (argc-3 < 2*nverts)
-			goto userr;
-		for (an = 0; an < nverts; an++) {
-			vert[an][0] = atof(argv[2*an+4]);
-			vert[an][1] = atof(argv[2*an+5]);
-		}
-		an = 2*nverts+4;
-	} else {
-		readverts(argv[3]);
-		an = 4;
-	}
-	if (nverts < 3) {
-		fprintf(stderr, "%s: not enough vertices\n", argv[0]);
-		exit(1);
-	}
-
-	for ( ; an < argc; an++) {
-		if (argv[an][0] != '-')
-			goto userr;
-		switch (argv[an][1]) {
-		case 'l':				/* length vector */
-			lvect[0] = atof(argv[++an]);
-			lvect[1] = atof(argv[++an]);
-			lvect[2] = atof(argv[++an]);
-			llen = sqrt(lvect[0]*lvect[0] + lvect[1]*lvect[1] +
-					lvect[2]*lvect[2]);
-			if (llen <= FTINY) {
-				fprintf(stderr, "%s: zero extrusion vector\n",
-						argv[0]);
-				exit(1);
-			}
-			break;
-		case 'r':				/* radius */
-			crad = atof(argv[++an]);
-			break;
-		case 'e':				/* ends */
-			do_ends = !do_ends;
-			break;
-		case 'c':				/* complete */
-			iscomplete = !iscomplete;
-			break;
-		default:
-			goto userr;
-		}
-	}
-	if (rounding) {
-		if (crad > fabs(lvect[2])) {
-			fprintf(stderr, "%s: rounding greater than height\n",
-					argv[0]);
-			exit(1);
-		}
-		compute_rounding();
-	}
-	printhead(argc, argv);
-
-	if (do_ends)
-		if (rounding)
-			printrends();
-		else
-			printends();
-
-	printsides(rounding);
-
-	exit(0);
-userr:
-	fprintf(stderr, "Usage: %s material name ", argv[0]);
-	fprintf(stderr, "{ - | vfile | N v1 v2 .. vN } ");
-	fprintf(stderr, "[-l lvect][-r radius][-c][-e]\n");
-	exit(1);
-}
-
-
+static void
 readverts(fname)		/* read vertices from a file */
 char  *fname;
 {
@@ -151,11 +57,80 @@ char  *fname;
 }
 
 
+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;
+	double  l, asum;
 
 	v0 = vert[nverts-1];
 	for (i = 0; i < nverts; i++) {		/* compute u[*] */
@@ -171,6 +146,7 @@ compute_rounding()		/* compute vectors for rounding op
 		u[i][1] /= l;
 		v0 = v1;
 	}
+	asum = 0.;
 	v1 = u[0];
 	for (i = nverts; i--; ) {		/* compute a[*] */
 		v0 = u[i];
@@ -183,14 +159,17 @@ compute_rounding()		/* compute vectors for rounding op
 			a[i] = 0.;
 		else {
 			a[i] = sqrt((1-l)/(1+l));
-			if ((v1[0]*v0[1]-v1[1]*v0[0] > 0.) != (lvect[2] > 0.))
+			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;
@@ -209,6 +188,7 @@ printends()			/* print ends of prism */
 }
 
 
+static void
 printrends()			/* print ends of prism with rounding */
 {
 	register int  i;
@@ -222,6 +202,15 @@ printrends()			/* print ends of prism with rounding */
 			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]);
@@ -236,27 +225,18 @@ printrends()			/* print ends of prism with rounding */
 		} else {
 			c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
 		}
-		if (a[i] > 0.) {
+		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], crad);
+				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", crad);
+		printf("\t%18.12g\n", lvdir*crad);
 		c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
 	}
-						/* 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]);
-	}
 						/* 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]);
@@ -273,21 +253,22 @@ printrends()			/* print ends of prism with rounding */
 		} else {
 			c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
 		}
-		if (a[i] > 0.) {
+		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], crad);
+				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", crad);
+		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;
 {
@@ -306,72 +287,7 @@ int  round;
 }
 
 
-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);
-}
-
-
-rside(n0, n1)			/* print side with rounded edge */
-register int  n0, n1;
-{
-	double  s, c, t[3];
-
-					/* compute tanget offset vector */
-	s = (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 (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", crad);
-}
-
-
+static void
 printhead(ac, av)		/* print command header */
 register int  ac;
 register char  **av;
@@ -383,3 +299,98 @@ register char  **av;
 	}
 	putchar('\n');
 }
+
+
+main(argc, argv)
+int  argc;
+char  **argv;
+{
+	int  an;
+	
+	if (argc < 4)
+		goto userr;
+
+	pmtype = argv[1];
+	pname = argv[2];
+
+	if (!strcmp(argv[3], "-")) {
+		readverts(NULL);
+		an = 4;
+	} else if (isdigit(argv[3][0])) {
+		nverts = atoi(argv[3]);
+		if (argc-3 < 2*nverts)
+			goto userr;
+		for (an = 0; an < nverts; an++) {
+			vert[an][0] = atof(argv[2*an+4]);
+			vert[an][1] = atof(argv[2*an+5]);
+		}
+		an = 2*nverts+4;
+	} else {
+		readverts(argv[3]);
+		an = 4;
+	}
+	if (nverts < 3) {
+		fprintf(stderr, "%s: not enough vertices\n", argv[0]);
+		exit(1);
+	}
+
+	for ( ; an < argc; an++) {
+		if (argv[an][0] != '-')
+			goto userr;
+		switch (argv[an][1]) {
+		case 'l':				/* length vector */
+			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;
+		case 'c':				/* complete */
+			iscomplete = !iscomplete;
+			break;
+		default:
+			goto userr;
+		}
+	}
+	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][-r radius][-c][-e]\n");
+	exit(1);
+}
+