--- ray/src/util/glarendx.c	1991/04/16 16:40:22	1.1
+++ ray/src/util/glarendx.c	1991/04/22 10:11:21	1.6
@@ -7,8 +7,8 @@ static char SCCSid[] = "$SunId$ LBL";
 /*
  * Compute Glare Index given by program name:
  *
- *	gnuth_dgr -	Gnuth discomfort glare rating
- *	gnuth_vcp -	Gnuth visual comfort probability
+ *	guth_dgr -	Guth discomfort glare rating
+ *	guth_vcp -	Guth visual comfort probability
  *
  *		12 April 1991	Greg Ward	EPFL
  */
@@ -21,15 +21,15 @@ extern double	erfc();
 double	posindex();
 int	headline();
 
-double	direct(), gnuth_dgr(), gnuth_vcp();
+double	direct(), guth_dgr(), guth_vcp();
 
 struct named_func {
 	char	*name;
 	double	(*func)();
 } all_funcs[] = {
 	{"direct", direct},
-	{"gnuth_dgr", gnuth_dgr},
-	{"gnuth_vcp", gnuth_vcp},
+	{"guth_dgr", guth_dgr},
+	{"guth_vcp", guth_vcp},
 	{NULL}
 };
 
@@ -53,7 +53,9 @@ int	print_header = 1;
 
 VIEW	midview = STDVIEW;
 
+int	wrongformat = 0;
 
+
 main(argc, argv)
 int	argc;
 char	*argv[];
@@ -89,14 +91,19 @@ char	*argv[];
 			exit(1);
 		}
 					/* read header */
-	getheader(stdin, headline);
+	getheader(stdin, headline, NULL);
+	if (wrongformat) {
+		fprintf(stderr, "%s: bad input format\n", progname);
+		exit(1);
+	}
 	if (print_header) {		/* add to header */
 		printargs(i, argv, stdout);
 		putchar('\n');
 	}
 					/* set view */
 	if (setview(&midview) != NULL) {
-		fprintf(stderr, "%s: bad view information in input\n");
+		fprintf(stderr, "%s: bad view information in input\n",
+				progname);
 		exit(1);
 	}
 					/* get findglare data */
@@ -119,10 +126,16 @@ userr:
 headline(s)			/* get line from header */
 char	*s;
 {
+	char	fmt[32];
+
 	if (print_header)		/* copy to output */
 		fputs(s, stdout);
 	if (!strncmp(s, VIEWSTR, VIEWSTRL))
 		sscanview(&midview, s+VIEWSTRL);
+	else if (isformat(s)) {
+		formatval(fmt, s);
+		wrongformat = strcmp(fmt, "ascii");
+	}
 }
 
 
@@ -155,6 +168,7 @@ read_input()			/* read glare sources from stdin */
 					&gs->dir[0], &gs->dir[1], &gs->dir[2],
 					&gs->dom, &gs->lum) != 5)
 				goto readerr;
+			normalize(gs->dir);
 			gs->next = all_srcs;
 			all_srcs = gs;
 			break;
@@ -225,6 +239,7 @@ struct glare_dir	*gd;
  * All vectors are assumed to be normalized.
  * This function is an implementation of the method proposed by
  * Robert Levin in his 1975 JIES article.
+ * This calculation presumes the view direction and up vectors perpendicular.
  * We return a value less than zero for improper positions.
  */
 
@@ -238,8 +253,14 @@ FVECT	sd, vd, vu;
 	d = DOT(sd,vd);
 	if (d <= 0.0)
 		return(-1.0);
+	if (d >= 1.0)
+		return(1.0);
 	sigma = acos(d) * (180./PI);
-	tau = acos(DOT(sd,vu)/sqrt(1.0-d*d)) * (180./PI);
+	d = DOT(sd,vu)/sqrt(1.0-d*d);
+	if (d >= 1.0)
+		tau = 0.0;
+	else
+		tau = acos(d) * (180./PI);
 	return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
 			+ sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
 		) );
@@ -247,28 +268,31 @@ FVECT	sd, vd, vu;
 
 
 double
-gnuth_dgr(gd)		/* compute Gnuth discomfort glare rating */
+guth_dgr(gd)		/* compute Guth discomfort glare rating */
 struct glare_dir	*gd;
 {
 #define q(w)	(20.4*w+1.52*pow(w,.2)-.075)
 	register struct glare_src	*gs;
 	double	p;
 	double	sum;
+	double	wtot, brsum;
 	int	n;
 
-	sum = 0.0; n = 0;
+	sum = wtot = brsum = 0.0; n = 0;
 	for (gs = all_srcs; gs != NULL; gs = gs->next) {
 		p = posindex(gs->dir, midview.vdir, midview.vup);
 		if (p <= FTINY)
 			continue;
 		sum += gs->lum * q(gs->dom) / p;
+		brsum += gs->lum * gs->dom;
+		wtot += gs->dom;
 		n++;
 	}
 	if (n == 0)
 		return(0.0);
 	else
 		return( pow(
-			.5*sum/pow(direct(gd)+gd->indirect,.44),
+			.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
 			pow((double)n, -.0914) ) );
 #undef q
 }
@@ -284,8 +308,8 @@ extern double	erf(), erfc();
 
 
 double
-gnuth_vcp(gd)		/* compute Gnuth visual comfort probability */
+guth_vcp(gd)		/* compute Guth visual comfort probability */
 struct glare_dir	*gd;
 {
-	return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd))));
+	return(100.*norm_integral(6.374-1.3227*log(guth_dgr(gd))));
 }