--- ray/src/util/glarendx.c	1991/04/16 16:40:22	1.1
+++ ray/src/util/glarendx.c	2018/08/02 18:33:50	2.13
@@ -1,38 +1,30 @@
-/* Copyright (c) 1991 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: glarendx.c,v 2.13 2018/08/02 18:33:50 greg Exp $";
 #endif
-
 /*
- * Compute Glare Index given by program name:
+ * Compute Glare Index given by program name or -t option:
  *
- *	gnuth_dgr -	Gnuth discomfort glare rating
- *	gnuth_vcp -	Gnuth visual comfort probability
+ *	dgi -		Daylight Glare Index
+ *	brs_gi -	Building Research Station Glare Index (Petherbridge
+ *							       & Hopkinson)
+ *	ugr -		Unified Glare Rating System (Fischer)
+ *	guth_dgr -	Guth discomfort glare rating
+ *	guth_vcp -	Guth visual comfort probability
+ *	cie_cgi -	CIE Glare Index (1983, due to Einhorn)
+ *	vert_dir -	Direct vertical illuminance
+ *	vert_ind -	Indirect vertical illuminance (from input)
+ *	vert_ill -	Total vertical illuminance
  *
  *		12 April 1991	Greg Ward	EPFL
+ *		19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
  */
 
+#include <string.h>
+
 #include "standard.h"
 #include "view.h"
 
-extern double	erfc();
 
-double	posindex();
-int	headline();
-
-double	direct(), gnuth_dgr(), gnuth_vcp();
-
-struct named_func {
-	char	*name;
-	double	(*func)();
-} all_funcs[] = {
-	{"direct", direct},
-	{"gnuth_dgr", gnuth_dgr},
-	{"gnuth_vcp", gnuth_vcp},
-	{NULL}
-};
-
 struct glare_src {
 	FVECT	dir;		/* source direction */
 	double	dom;		/* solid angle */
@@ -46,6 +38,42 @@ struct glare_dir {
 	struct glare_dir	*next;
 } *all_dirs = NULL;
 
+typedef double gdfun(struct glare_dir *gd);
+
+static gdfun dgi;
+static gdfun brs_gi;
+static gdfun ugr;
+static gdfun guth_dgr;
+static gdfun guth_vcp;
+static gdfun cie_cgi;
+static gdfun direct;
+static gdfun indirect;
+static gdfun total;
+
+static gethfunc headline;
+static void init(void);
+static void read_input(void);
+static void print_values(gdfun *func);
+static double posindex(FVECT sd, FVECT vd, FVECT vu);
+
+
+struct named_func {
+	char	*name;
+	gdfun	*func;
+	char	*descrip;
+} all_funcs[] = {
+	{"dgi", dgi, "Daylight Glare Index"},
+	{"brs_gi", brs_gi, "BRS Glare Index"},
+	{"ugr", ugr, "Unified Glare Rating"},
+	{"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
+	{"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
+	{"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
+	{"vert_dir", direct, "Direct Vertical Illuminance"},
+	{"vert_ill", total, "Total Vertical Illuminance"},
+	{"vert_ind", indirect, "Indirect Vertical Illuminance"},
+	{NULL}
+};
+
 #define newp(type)	(type *)malloc(sizeof(type))
 
 char	*progname;
@@ -53,18 +81,21 @@ int	print_header = 1;
 
 VIEW	midview = STDVIEW;
 
+int	wrongformat = 0;
 
-main(argc, argv)
-int	argc;
-char	*argv[];
+
+int
+main(
+	int	argc,
+	char	*argv[]
+)
 {
-	extern char	*rindex();
 	struct named_func	*funp;
 	char	*progtail;
 	int	i;
 					/* get program name */
 	progname = argv[0];
-	progtail = rindex(progname, '/');	/* final component */
+	progtail = strrchr(progname, '/');	/* final component */
 	if (progtail == NULL)
 		progtail = progname;
 	else
@@ -88,53 +119,76 @@ char	*argv[];
 			perror(argv[i]);
 			exit(1);
 		}
-					/* read header */
-	getheader(stdin, headline);
-	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");
-		exit(1);
-	}
-					/* get findglare data */
-	read_input();
-					/* find calculation */
+					/* find and run calculation */
 	for (funp = all_funcs; funp->name != NULL; funp++)
 		if (!strcmp(funp->name, progtail)) {
+			init();
+			read_input();
+			if (print_header) {
+				printargs(i, argv, stdout);
+				putchar('\n');
+			}
 			print_values(funp->func);
 			exit(0);		/* we're done */
 		}
 					/* invalid function */
-	fprintf(stderr, "%s: unknown function!\n", progtail);
-	exit(1);
 userr:
-	fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
+	fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
+	fprintf(stderr, "\twhere 'type' is one of the following:\n");
+	for (funp = all_funcs; funp->name != NULL; funp++)
+		fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
 	exit(1);
 }
 
 
-headline(s)			/* get line from header */
-char	*s;
+static int
+headline(			/* get line from header */
+	char	*s,
+	void	*p
+)
 {
+	char	fmt[MAXFMTLEN];
+
 	if (print_header)		/* copy to output */
 		fputs(s, stdout);
-	if (!strncmp(s, VIEWSTR, VIEWSTRL))
-		sscanview(&midview, s+VIEWSTRL);
+	if (isview(s))
+		sscanview(&midview, s);
+	else if (isformat(s)) {
+		formatval(fmt, s);
+		wrongformat = strcmp(fmt, "ascii");
+	}
+	return(0);
 }
 
 
-read_input()			/* read glare sources from stdin */
+static void
+init(void)				/* initialize calculation */
 {
+					/* read header */
+	getheader(stdin, headline, NULL);
+	if (wrongformat) {
+		fprintf(stderr, "%s: bad input format\n", progname);
+		exit(1);
+	}
+					/* set view */
+	if (setview(&midview) != NULL) {
+		fprintf(stderr, "%s: bad view information in input\n",
+				progname);
+		exit(1);
+	}
+}
+
+
+static void
+read_input(void)			/* read glare sources from stdin */
+{
 #define	S_SEARCH	0
 #define S_SOURCE	1
 #define S_DIREC		2
 	int	state = S_SEARCH;
 	char	buf[128];
-	register struct glare_src	*gs;
-	register struct glare_dir	*gd;
+	struct glare_src	*gs;
+	struct glare_dir	*gd;
 
 	while (fgets(buf, sizeof(buf), stdin) != NULL)
 		switch (state) {
@@ -151,10 +205,12 @@ read_input()			/* read glare sources from stdin */
 			}
 			if ((gs = newp(struct glare_src)) == NULL)
 				goto memerr;
-			if (sscanf(buf, "%lf %lf %lf %lf %lf",
-					&gs->dir[0], &gs->dir[1], &gs->dir[2],
-					&gs->dom, &gs->lum) != 5)
+			if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
+						&gs->dir[2]) != 3 ||
+					sscanf(sskip2(buf, 3), "%lf %lf",
+						&gs->dom, &gs->lum) != 2)
 				goto readerr;
+			normalize(gs->dir);
 			gs->next = all_srcs;
 			all_srcs = gs;
 			break;
@@ -186,23 +242,26 @@ readerr:
 }
 
 
-print_values(funp)		/* print out calculations */
-double	(*funp)();
+static void
+print_values(		/* print out calculations */
+	gdfun *func
+)
 {
-	register struct glare_dir	*gd;
+	struct glare_dir	*gd;
 
 	for (gd = all_dirs; gd != NULL; gd = gd->next)
-		printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
+		printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
 }
 
 
-double
-direct(gd)			/* compute direct illuminance */
-struct glare_dir	*gd;
+static double
+direct(			/* compute direct vertical illuminance */
+	struct glare_dir	*gd
+)
 {
 	FVECT	mydir;
 	double	d, dval;
-	register struct glare_src	*gs;
+	struct glare_src	*gs;
 
 	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
 	dval = 0.0;
@@ -215,6 +274,24 @@ struct glare_dir	*gd;
 }
 
 
+static double
+indirect(			/* return indirect vertical illuminance */
+	struct glare_dir	*gd
+)
+{
+	return(gd->indirect);
+}
+
+
+static double
+total(			/* return total vertical illuminance */
+	struct glare_dir	*gd
+)
+{
+	return(direct(gd)+gd->indirect);
+}
+
+
 /*
  * posindex -	compute glare position index from:
  *
@@ -225,12 +302,16 @@ 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.
  */
 
-double
-posindex(sd, vd, vu)			/* compute position index */
-FVECT	sd, vd, vu;
+static double
+posindex(			/* compute position index */
+	FVECT	sd,
+	FVECT	vd,
+	FVECT	vu
+)
 {
 	double	sigma, tau;
 	double	d;
@@ -238,44 +319,128 @@ 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 = fabs(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 )
 		) );
 }
 
 
-double
-gnuth_dgr(gd)		/* compute Gnuth discomfort glare rating */
-struct glare_dir	*gd;
+static double
+dgi(		/* compute Daylight Glare Index */
+	struct glare_dir	*gd
+)
 {
+	struct glare_src	*gs;
+	FVECT	mydir,testdir[7],vhor;
+	double	r,posn,omega,p[7],sum;
+	int	i,n;
+
+	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
+	sum = 0.0; n = 0;
+	for (gs = all_srcs; gs != NULL; gs = gs->next) {
+
+		/* compute 1/p^2 weighted solid angle of the source */
+		r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
+		fcross(vhor,gs->dir,midview.vup);
+		normalize(vhor);
+		VCOPY(testdir[0],gs->dir);
+		fvsum(testdir[1],gs->dir,vhor,r);
+		fvsum(testdir[2],gs->dir,vhor,0.5*r);
+		fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
+		fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
+		fvsum(testdir[3],gs->dir,vhor,-r);
+		fvsum(testdir[4],gs->dir,vhor,-0.5*r);
+		fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
+		fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
+		for (i = 0; i < 7; i++) {
+			normalize(testdir[i]);
+			posn = posindex(testdir[i],mydir,midview.vup);
+			if (posn <= FTINY)
+				p[i] = 0.0;
+			else
+				p[i] = 1./(posn*posn);
+		}
+		r = 1-gs->dom/2./PI;
+		omega = gs->dom*p[0];
+		omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
+		omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
+			  +0.3334*(p[2]+p[4]+p[5]+p[6]));
+
+		sum += pow(gs->lum,1.6) * pow(omega,0.8) /
+		       (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
+		n++;
+	}
+	if (n == 0)
+		return(0.0);
+	return( 10*log10(0.478*sum) );
+}
+
+
+static double
+brs_gi(		/* compute BRS Glare Index */
+	struct glare_dir	*gd
+)
+{
+	struct glare_src	*gs;
+	FVECT	mydir;
+	double	p;
+	double	sum;
+
+	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
+	sum = 0.0;
+	for (gs = all_srcs; gs != NULL; gs = gs->next) {
+		p = posindex(gs->dir, mydir, midview.vup);
+		if (p <= FTINY)
+			continue;
+		sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
+	}
+	if (sum <= FTINY)
+		return(0.0);
+	sum /= gd->indirect/PI;
+	return(10*log10(0.478*sum));
+}
+
+
+static double
+guth_dgr(		/* 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;
+	struct glare_src	*gs;
+	FVECT	mydir;
 	double	p;
 	double	sum;
+	double	wtot, brsum;
 	int	n;
 
-	sum = 0.0; n = 0;
+	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
+	sum = wtot = brsum = 0.0; n = 0;
 	for (gs = all_srcs; gs != NULL; gs = gs->next) {
-		p = posindex(gs->dir, midview.vdir, midview.vup);
+		p = posindex(gs->dir, mydir, 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),
+	return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
 			pow((double)n, -.0914) ) );
 #undef q
 }
 
 
-extern double	erf(), erfc();
-
 #ifndef M_SQRT2
 #define	M_SQRT2	1.41421356237309504880
 #endif
@@ -283,9 +448,66 @@ extern double	erf(), erfc();
 #define norm_integral(z)	(1.-.5*erfc((z)/M_SQRT2))
 
 
-double
-gnuth_vcp(gd)		/* compute Gnuth visual comfort probability */
-struct glare_dir	*gd;
+static double
+guth_vcp(		/* compute Guth visual comfort probability */
+	struct glare_dir	*gd
+)
 {
-	return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd))));
+	extern double	erfc();
+	double	dgr;
+
+	dgr = guth_dgr(gd);
+	if (dgr <= FTINY)
+		return(100.0);
+	return(100.*norm_integral(6.374-1.3227*log(dgr)));
+}
+
+
+static double
+cie_cgi(		/* compute CIE Glare Index */
+	struct glare_dir	*gd
+)
+{
+	struct glare_src	*gs;
+	FVECT	mydir;
+	double	dillum;
+	double	p;
+	double	sum;
+
+	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
+	sum = 0.0;
+	for (gs = all_srcs; gs != NULL; gs = gs->next) {
+		p = posindex(gs->dir, mydir, midview.vup);
+		if (p <= FTINY)
+			continue;
+		sum += gs->lum*gs->lum * gs->dom / (p*p);
+	}
+	if (sum <= FTINY)
+		return(0.0);
+	dillum = direct(gd);
+	return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
+}
+
+
+static double
+ugr(		/* compute Unified Glare Rating */
+	struct glare_dir	*gd
+)
+{
+	struct glare_src	*gs;
+	FVECT	mydir;
+	double	p;
+	double	sum;
+
+	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
+	sum = 0.0;
+	for (gs = all_srcs; gs != NULL; gs = gs->next) {
+		p = posindex(gs->dir, mydir, midview.vup);
+		if (p <= FTINY)
+			continue;
+		sum += gs->lum*gs->lum * gs->dom / (p*p);
+	}
+	if (sum <= FTINY)
+		return(0.0);
+	return(8.*log10(0.25*sum*PI/gd->indirect));
 }