--- ray/src/util/glarendx.c	1992/11/19 09:59:08	2.3
+++ ray/src/util/glarendx.c	1994/03/24 11:06:59	2.4
@@ -7,6 +7,10 @@ static char SCCSid[] = "$SunId$ LBL";
 /*
  * Compute Glare Index given by program name or -t option:
  *
+ *	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)
@@ -15,23 +19,27 @@ static char SCCSid[] = "$SunId$ LBL";
  *	vert_ill -	Total vertical illuminance
  *
  *		12 April 1991	Greg Ward	EPFL
+ *		19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
  */
-
+ 
 #include "standard.h"
 #include "view.h"
-
+ 
 extern double	erfc();
-
+ 
 double	posindex();
-
+ 
 double	direct(), guth_dgr(), guth_vcp(), cie_cgi(),
-	indirect(), total();
-
+	indirect(), total(), dgi(), brs_gi(), ugr();
+ 
 struct named_func {
 	char	*name;
 	double	(*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"},
@@ -40,30 +48,30 @@ struct named_func {
 	{"vert_ind", indirect, "Indirect Vertical Illuminance"},
 	{NULL}
 };
-
+ 
 struct glare_src {
 	FVECT	dir;		/* source direction */
 	double	dom;		/* solid angle */
 	double	lum;		/* average luminance */
 	struct glare_src	*next;
 } *all_srcs = NULL;
-
+ 
 struct glare_dir {
 	double	ang;		/* angle (in radians) */
 	double	indirect;	/* indirect illuminance */
 	struct glare_dir	*next;
 } *all_dirs = NULL;
-
+ 
 #define newp(type)	(type *)malloc(sizeof(type))
-
+ 
 char	*progname;
 int	print_header = 1;
-
+ 
 VIEW	midview = STDVIEW;
-
+ 
 int	wrongformat = 0;
-
-
+ 
+ 
 main(argc, argv)
 int	argc;
 char	*argv[];
@@ -118,13 +126,13 @@ userr:
 		fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
 	exit(1);
 }
-
-
+ 
+ 
 headline(s)			/* get line from header */
 char	*s;
 {
 	char	fmt[32];
-
+ 
 	if (print_header)		/* copy to output */
 		fputs(s, stdout);
 	if (isview(s))
@@ -134,8 +142,8 @@ char	*s;
 		wrongformat = strcmp(fmt, "ascii");
 	}
 }
-
-
+ 
+ 
 init()				/* initialize calculation */
 {
 					/* read header */
@@ -151,8 +159,8 @@ init()				/* initialize calculation */
 		exit(1);
 	}
 }
-
-
+ 
+ 
 read_input()			/* read glare sources from stdin */
 {
 #define	S_SEARCH	0
@@ -162,7 +170,7 @@ read_input()			/* read glare sources from stdin */
 	char	buf[128];
 	register struct glare_src	*gs;
 	register struct glare_dir	*gd;
-
+ 
 	while (fgets(buf, sizeof(buf), stdin) != NULL)
 		switch (state) {
 		case S_SEARCH:
@@ -212,18 +220,18 @@ readerr:
 #undef S_SOURCE
 #undef S_DIREC
 }
-
-
+ 
+ 
 print_values(funp)		/* print out calculations */
 double	(*funp)();
 {
 	register 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));
 }
-
-
+ 
+ 
 double
 direct(gd)			/* compute direct vertical illuminance */
 struct glare_dir	*gd;
@@ -231,7 +239,7 @@ struct glare_dir	*gd;
 	FVECT	mydir;
 	double	d, dval;
 	register struct glare_src	*gs;
-
+ 
 	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
 	dval = 0.0;
 	for (gs = all_srcs; gs != NULL; gs = gs->next) {
@@ -241,24 +249,24 @@ struct glare_dir	*gd;
 	}
 	return(dval);
 }
-
-
+ 
+ 
 double
 indirect(gd)			/* return indirect vertical illuminance */
 struct glare_dir	*gd;
 {
 	return(gd->indirect);
 }
-
-
+ 
+ 
 double
 total(gd)			/* return total vertical illuminance */
 struct glare_dir	*gd;
 {
 	return(direct(gd)+gd->indirect);
 }
-
-
+ 
+ 
 /*
  * posindex -	compute glare position index from:
  *
@@ -272,34 +280,103 @@ struct glare_dir	*gd;
  * 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;
 {
 	double	sigma, tau;
 	double	d;
-
+ 
 	d = DOT(sd,vd);
 	if (d <= 0.0)
 		return(-1.0);
 	if (d >= 1.0)
 		return(1.0);
 	sigma = acos(d) * (180./PI);
-	d = DOT(sd,vu)/sqrt(1.0-d*d);
+	d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
 	if (d >= 1.0)
 		tau = 0.0;
-	else if (d <= -1.0)
-		tau = 180.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
+dgi(gd)		/* compute Daylight Glare Index */
+struct glare_dir	*gd;
+{
+	register struct glare_src	*gs;
+	FVECT	mydir,testdir[7],vhor;
+	double	r,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]);
+			p[i] = pow(posindex(testdir[i],mydir,midview.vup),-2.0);
+			if (p[i] <= FTINY) p[i] = 0.0;
+		}
+		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) );
+}
+ 
+
 double
+brs_gi(gd)		/* compute BRS Glare Index */
+struct glare_dir	*gd;
+{
+	register 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));
+}
+
+
+double
 guth_dgr(gd)		/* compute Guth discomfort glare rating */
 struct glare_dir	*gd;
 {
@@ -310,7 +387,7 @@ struct glare_dir	*gd;
 	double	sum;
 	double	wtot, brsum;
 	int	n;
-
+ 
 	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
 	sum = wtot = brsum = 0.0; n = 0;
 	for (gs = all_srcs; gs != NULL; gs = gs->next) {
@@ -324,35 +401,34 @@ struct glare_dir	*gd;
 	}
 	if (n == 0)
 		return(0.0);
-
 	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
-
+ 
 #define norm_integral(z)	(1.-.5*erfc((z)/M_SQRT2))
-
-
+ 
+ 
 double
 guth_vcp(gd)		/* compute Guth visual comfort probability */
 struct glare_dir	*gd;
 {
 	double	dgr;
-
+ 
 	dgr = guth_dgr(gd);
 	if (dgr <= FTINY)
 		return(100.0);
 	return(100.*norm_integral(6.374-1.3227*log(dgr)));
 }
-
-
+ 
+ 
 double
 cie_cgi(gd)		/* compute CIE Glare Index */
 struct glare_dir	*gd;
@@ -362,7 +438,7 @@ struct glare_dir	*gd;
 	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) {
@@ -375,4 +451,27 @@ struct glare_dir	*gd;
 		return(0.0);
 	dillum = direct(gd);
 	return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
+}
+
+
+double
+ugr(gd)		/* compute Unified Glare Rating */
+struct glare_dir	*gd;
+{
+	register 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));
 }