src/util/glarendx.c

#ifndef lint
static const char       RCSid[] = "$Id: glarendx.c,v 2.8 2003/06/30 14:59:13 schorsch Exp $";
#endif
/*
 * 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)
 *      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"
 
 
double  posindex();
 
double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
        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"},
        {"vert_dir", direct, "Direct Vertical Illuminance"},
        {"vert_ill", total, "Total Vertical Illuminance"},
        {"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[];
{
        struct named_func       *funp;
        char    *progtail;
        int     i;
                                        /* get program name */
        progname = argv[0];
        progtail = strrchr(progname, '/');      /* final component */
        if (progtail == NULL)
                progtail = progname;
        else
                progtail++;
                                        /* get options */
        for (i = 1; i < argc && argv[i][0] == '-'; i++)
                switch (argv[i][1]) {
                case 't':
                        progtail = argv[++i];
                        break;
                case 'h':
                        print_header = 0;
                        break;
                default:
                        goto userr;
                }
        if (i < argc-1)
                goto userr;
        if (i == argc-1)                /* open file */
                if (freopen(argv[i], "r", stdin) == NULL) {
                        perror(argv[i]);
                        exit(1);
                }
                                        /* 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 */
userr:
        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);
}
 
 
int
headline(s)                     /* get line from header */
char    *s;
{
        char    fmt[32];
 
        if (print_header)               /* copy to output */
                fputs(s, stdout);
        if (isview(s))
                sscanview(&midview, s);
        else if (isformat(s)) {
                formatval(fmt, s);
                wrongformat = strcmp(fmt, "ascii");
        }
        return(0);
}
 
 
init()                          /* 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);
        }
}
 
 
read_input()                    /* 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;
 
        while (fgets(buf, sizeof(buf), stdin) != NULL)
                switch (state) {
                case S_SEARCH:
                        if (!strcmp(buf, "BEGIN glare source\n"))
                                state = S_SOURCE;
                        else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
                                state = S_DIREC;
                        break;
                case S_SOURCE:
                        if (!strncmp(buf, "END", 3)) {
                                state = S_SEARCH;
                                break;
                        }
                        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)
                                goto readerr;
                        normalize(gs->dir);
                        gs->next = all_srcs;
                        all_srcs = gs;
                        break;
                case S_DIREC:
                        if (!strncmp(buf, "END", 3)) {
                                state = S_SEARCH;
                                break;
                        }
                        if ((gd = newp(struct glare_dir)) == NULL)
                                goto memerr;
                        if (sscanf(buf, "%lf %lf",
                                        &gd->ang, &gd->indirect) != 2)
                                goto readerr;
                        gd->ang *= PI/180.0;    /* convert to radians */
                        gd->next = all_dirs;
                        all_dirs = gd;
                        break;
                }
        return;
memerr:
        perror(progname);
        exit(1);
readerr:
        fprintf(stderr, "%s: read error on input\n", progname);
        exit(1);
#undef S_SEARCH
#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;
{
        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) {
                d = DOT(mydir,gs->dir);
                if (d > FTINY)
                        dval += d * gs->dom * gs->lum;
        }
        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:
 *
 *      Source Direction
 *      View Direction
 *      View Up Direction
 *
 * 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;
{
        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 = 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
dgi(gd)         /* compute Daylight Glare Index */
struct glare_dir        *gd;
{
        register 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) );
}
 

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;
{
#define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
        register struct glare_src       *gs;
        FVECT   mydir;
        double  p;
        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) {
                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);
        return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
                        pow((double)n, -.0914) ) );
#undef q
}
 
 
#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;
{
        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)));
}
 
 
double
cie_cgi(gd)             /* compute CIE Glare Index */
struct glare_dir        *gd;
{
        register 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)));
}


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));
}
Revision:	2.9
Committed:	Wed Nov 19 16:21:28 2003 UTC (20 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.8:	+7 -4 lines
Log Message:	Fixed bug in glarendx pointed out by Phillip Greenup
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.9	static const char RCSid[] = "$Id: glarendx.c,v 2.8 2003/06/30 14:59:13 schorsch Exp $";
3	greg	1.1	#endif
4			/*
5	greg	1.9	* Compute Glare Index given by program name or -t option:
6	greg	1.1	*
7	greg	2.4	* dgi - Daylight Glare Index
8			* brs_gi - Building Research Station Glare Index (Petherbridge
9			* & Hopkinson)
10			* ugr - Unified Glare Rating System (Fischer)
11	greg	1.2	* guth_dgr - Guth discomfort glare rating
12			* guth_vcp - Guth visual comfort probability
13	greg	1.9	* cie_cgi - CIE Glare Index (1983, due to Einhorn)
14			* vert_dir - Direct vertical illuminance
15			* vert_ind - Indirect vertical illuminance (from input)
16			* vert_ill - Total vertical illuminance
17	greg	1.1	*
18			* 12 April 1991 Greg Ward EPFL
19	greg	2.4	* 19 April 1993 R. Compagnon EPFL (added dgi, brs_gi, ugr)
20	greg	1.1	*/
21	greg	2.4
22	schorsch	2.8	#include <string.h>
23
24	greg	1.1	#include "standard.h"
25			#include "view.h"
26	greg	2.4
27
28	greg	1.1	double posindex();
29	greg	2.4
30	greg	1.9	double direct(), guth_dgr(), guth_vcp(), cie_cgi(),
31	greg	2.4	indirect(), total(), dgi(), brs_gi(), ugr();
32
33	greg	1.1	struct named_func {
34			char *name;
35			double (*func)();
36	greg	1.9	char *descrip;
37	greg	1.1	} all_funcs[] = {
38	greg	2.4	{"dgi", dgi, "Daylight Glare Index"},
39			{"brs_gi", brs_gi, "BRS Glare Index"},
40			{"ugr", ugr, "Unified Glare Rating"},
41	greg	1.9	{"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
42			{"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
43			{"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
44			{"vert_dir", direct, "Direct Vertical Illuminance"},
45			{"vert_ill", total, "Total Vertical Illuminance"},
46			{"vert_ind", indirect, "Indirect Vertical Illuminance"},
47	greg	1.1	{NULL}
48			};
49	greg	2.4
50	greg	1.1	struct glare_src {
51			FVECT dir; /* source direction */
52			double dom; /* solid angle */
53			double lum; /* average luminance */
54			struct glare_src *next;
55			} *all_srcs = NULL;
56	greg	2.4
57	greg	1.1	struct glare_dir {
58			double ang; /* angle (in radians) */
59			double indirect; /* indirect illuminance */
60			struct glare_dir *next;
61			} *all_dirs = NULL;
62	greg	2.4
63	greg	1.1	#define newp(type) (type *)malloc(sizeof(type))
64	greg	2.4
65	greg	1.1	char *progname;
66			int print_header = 1;
67	greg	2.4
68	greg	1.1	VIEW midview = STDVIEW;
69	greg	2.4
70	greg	1.4	int wrongformat = 0;
71	greg	2.4
72
73	greg	1.1	main(argc, argv)
74			int argc;
75			char *argv[];
76			{
77			struct named_func *funp;
78			char *progtail;
79			int i;
80			/* get program name */
81			progname = argv[0];
82	schorsch	2.8	progtail = strrchr(progname, '/'); /* final component */
83	greg	1.1	if (progtail == NULL)
84			progtail = progname;
85			else
86			progtail++;
87			/* get options */
88			for (i = 1; i < argc && argv[i][0] == '-'; i++)
89			switch (argv[i][1]) {
90			case 't':
91			progtail = argv[++i];
92			break;
93			case 'h':
94			print_header = 0;
95			break;
96			default:
97			goto userr;
98			}
99			if (i < argc-1)
100			goto userr;
101			if (i == argc-1) /* open file */
102			if (freopen(argv[i], "r", stdin) == NULL) {
103			perror(argv[i]);
104			exit(1);
105			}
106	greg	1.9	/* find and run calculation */
107	greg	1.1	for (funp = all_funcs; funp->name != NULL; funp++)
108			if (!strcmp(funp->name, progtail)) {
109	greg	1.9	init();
110			read_input();
111			if (print_header) {
112			printargs(i, argv, stdout);
113			putchar('\n');
114			}
115	greg	1.1	print_values(funp->func);
116			exit(0); /* we're done */
117			}
118			/* invalid function */
119			userr:
120	greg	1.9	fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
121			fprintf(stderr, "\twhere 'type' is one of the following:\n");
122			for (funp = all_funcs; funp->name != NULL; funp++)
123			fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
124	greg	1.1	exit(1);
125			}
126	greg	2.4
127
128	gwlarson	2.6	int
129	greg	1.1	headline(s) /* get line from header */
130			char *s;
131			{
132	greg	1.4	char fmt[32];
133	greg	2.4
134	greg	1.1	if (print_header) /* copy to output */
135			fputs(s, stdout);
136	greg	2.2	if (isview(s))
137			sscanview(&midview, s);
138	greg	1.4	else if (isformat(s)) {
139			formatval(fmt, s);
140	greg	1.6	wrongformat = strcmp(fmt, "ascii");
141	greg	1.4	}
142	gwlarson	2.6	return(0);
143	greg	1.1	}
144	greg	2.4
145
146	greg	1.9	init() /* initialize calculation */
147			{
148			/* read header */
149			getheader(stdin, headline, NULL);
150			if (wrongformat) {
151			fprintf(stderr, "%s: bad input format\n", progname);
152			exit(1);
153			}
154			/* set view */
155			if (setview(&midview) != NULL) {
156			fprintf(stderr, "%s: bad view information in input\n",
157			progname);
158			exit(1);
159			}
160			}
161	greg	2.4
162
163	greg	1.1	read_input() /* read glare sources from stdin */
164			{
165			#define S_SEARCH 0
166			#define S_SOURCE 1
167			#define S_DIREC 2
168			int state = S_SEARCH;
169			char buf[128];
170			register struct glare_src *gs;
171			register struct glare_dir *gd;
172	greg	2.4
173	greg	1.1	while (fgets(buf, sizeof(buf), stdin) != NULL)
174			switch (state) {
175			case S_SEARCH:
176			if (!strcmp(buf, "BEGIN glare source\n"))
177			state = S_SOURCE;
178			else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
179			state = S_DIREC;
180			break;
181			case S_SOURCE:
182			if (!strncmp(buf, "END", 3)) {
183			state = S_SEARCH;
184			break;
185			}
186			if ((gs = newp(struct glare_src)) == NULL)
187			goto memerr;
188			if (sscanf(buf, "%lf %lf %lf %lf %lf",
189			&gs->dir[0], &gs->dir[1], &gs->dir[2],
190			&gs->dom, &gs->lum) != 5)
191			goto readerr;
192	greg	1.2	normalize(gs->dir);
193	greg	1.1	gs->next = all_srcs;
194			all_srcs = gs;
195			break;
196			case S_DIREC:
197			if (!strncmp(buf, "END", 3)) {
198			state = S_SEARCH;
199			break;
200			}
201			if ((gd = newp(struct glare_dir)) == NULL)
202			goto memerr;
203			if (sscanf(buf, "%lf %lf",
204			&gd->ang, &gd->indirect) != 2)
205			goto readerr;
206			gd->ang = PI/180.0; / convert to radians */
207			gd->next = all_dirs;
208			all_dirs = gd;
209			break;
210			}
211			return;
212			memerr:
213			perror(progname);
214			exit(1);
215			readerr:
216			fprintf(stderr, "%s: read error on input\n", progname);
217			exit(1);
218			#undef S_SEARCH
219			#undef S_SOURCE
220			#undef S_DIREC
221			}
222	greg	2.4
223
224	greg	1.1	print_values(funp) /* print out calculations */
225			double (*funp)();
226			{
227			register struct glare_dir *gd;
228	greg	2.4
229	greg	1.1	for (gd = all_dirs; gd != NULL; gd = gd->next)
230			printf("%f\t%f\n", gd->ang(180.0/PI), (funp)(gd));
231			}
232	greg	2.4
233
234	greg	1.1	double
235	greg	1.9	direct(gd) /* compute direct vertical illuminance */
236	greg	1.1	struct glare_dir *gd;
237			{
238			FVECT mydir;
239			double d, dval;
240			register struct glare_src *gs;
241	greg	2.4
242	greg	1.1	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
243			dval = 0.0;
244			for (gs = all_srcs; gs != NULL; gs = gs->next) {
245			d = DOT(mydir,gs->dir);
246			if (d > FTINY)
247			dval += d * gs->dom * gs->lum;
248			}
249			return(dval);
250	greg	1.9	}
251	greg	2.4
252
253	greg	1.9	double
254			indirect(gd) /* return indirect vertical illuminance */
255			struct glare_dir *gd;
256			{
257			return(gd->indirect);
258			}
259	greg	2.4
260
261	greg	1.9	double
262			total(gd) /* return total vertical illuminance */
263			struct glare_dir *gd;
264			{
265			return(direct(gd)+gd->indirect);
266	greg	1.1	}
267	greg	2.4
268
269	greg	1.1	/*
270			* posindex - compute glare position index from:
271			*
272			* Source Direction
273			* View Direction
274			* View Up Direction
275			*
276			* All vectors are assumed to be normalized.
277			* This function is an implementation of the method proposed by
278			* Robert Levin in his 1975 JIES article.
279	greg	1.5	* This calculation presumes the view direction and up vectors perpendicular.
280	greg	1.1	* We return a value less than zero for improper positions.
281			*/
282	greg	2.4
283	greg	1.1	double
284			posindex(sd, vd, vu) /* compute position index */
285			FVECT sd, vd, vu;
286			{
287			double sigma, tau;
288			double d;
289	greg	2.4
290	greg	1.1	d = DOT(sd,vd);
291			if (d <= 0.0)
292			return(-1.0);
293	greg	1.5	if (d >= 1.0)
294			return(1.0);
295	greg	1.1	sigma = acos(d) * (180./PI);
296	greg	2.4	d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
297	greg	1.5	if (d >= 1.0)
298			tau = 0.0;
299			else
300			tau = acos(d) * (180./PI);
301	greg	1.1	return( exp( sigma( (35.2 - tau.31889 - 1.22exp(-.22222tau))*1e-3
302			+ sigma(21. + tau(.26667 + tau-.002963))1e-5 )
303			) );
304			}
305	greg	2.4
306
307			double
308			dgi(gd) /* compute Daylight Glare Index */
309			struct glare_dir *gd;
310			{
311			register struct glare_src *gs;
312			FVECT mydir,testdir[7],vhor;
313	greg	2.9	double r,posn,omega,p[7],sum;
314	greg	2.4	int i,n;
315
316			spinvector(mydir, midview.vdir, midview.vup, gd->ang);
317			sum = 0.0; n = 0;
318			for (gs = all_srcs; gs != NULL; gs = gs->next) {
319	greg	1.1
320	greg	2.4	/* compute 1/p^2 weighted solid angle of the source */
321			r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
322			fcross(vhor,gs->dir,midview.vup);
323			normalize(vhor);
324			VCOPY(testdir[0],gs->dir);
325			fvsum(testdir[1],gs->dir,vhor,r);
326			fvsum(testdir[2],gs->dir,vhor,0.5*r);
327			fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
328			fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
329			fvsum(testdir[3],gs->dir,vhor,-r);
330			fvsum(testdir[4],gs->dir,vhor,-0.5*r);
331			fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
332			fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
333			for (i = 0; i < 7; i++) {
334			normalize(testdir[i]);
335	greg	2.9	posn = posindex(testdir[i],mydir,midview.vup);
336			if (posn <= FTINY)
337			p[i] = 0.0;
338			else
339			p[i] = 1./(posn*posn);
340	greg	2.4	}
341			r = 1-gs->dom/2./PI;
342			omega = gs->dom*p[0];
343			omega += (rPI(1+1/r/r)-2PI)(-p[0]+(p[1]+p[2])*0.5);
344			omega += (2PI-rPI(1+1/r/r))(-p[0]-0.1667*(p[1]+p[3])
345			+0.3334*(p[2]+p[4]+p[5]+p[6]));
346	greg	1.1
347	greg	2.4	sum += pow(gs->lum,1.6) * pow(omega,0.8) /
348			(gd->indirect/PI + 0.07sqrt(gs->dom)gs->lum);
349			n++;
350			}
351			if (n == 0)
352			return(0.0);
353			return( 10log10(0.478sum) );
354			}
355
356
357	greg	1.1	double
358	greg	2.4	brs_gi(gd) /* compute BRS Glare Index */
359			struct glare_dir *gd;
360			{
361			register struct glare_src *gs;
362			FVECT mydir;
363			double p;
364			double sum;
365
366			spinvector(mydir, midview.vdir, midview.vup, gd->ang);
367			sum = 0.0;
368			for (gs = all_srcs; gs != NULL; gs = gs->next) {
369			p = posindex(gs->dir, mydir, midview.vup);
370			if (p <= FTINY)
371			continue;
372			sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
373			}
374			if (sum <= FTINY)
375			return(0.0);
376			sum /= gd->indirect/PI;
377			return(10log10(0.478sum));
378			}
379
380
381			double
382	greg	1.2	guth_dgr(gd) /* compute Guth discomfort glare rating */
383	greg	1.1	struct glare_dir *gd;
384			{
385			#define q(w) (20.4w+1.52pow(w,.2)-.075)
386			register struct glare_src *gs;
387	greg	1.7	FVECT mydir;
388	greg	1.1	double p;
389			double sum;
390	greg	1.3	double wtot, brsum;
391	greg	1.1	int n;
392	greg	2.4
393	greg	1.7	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
394	greg	1.3	sum = wtot = brsum = 0.0; n = 0;
395	greg	1.1	for (gs = all_srcs; gs != NULL; gs = gs->next) {
396	greg	1.7	p = posindex(gs->dir, mydir, midview.vup);
397	greg	1.1	if (p <= FTINY)
398			continue;
399			sum += gs->lum * q(gs->dom) / p;
400	greg	1.3	brsum += gs->lum * gs->dom;
401			wtot += gs->dom;
402	greg	1.1	n++;
403			}
404			if (n == 0)
405			return(0.0);
406	greg	1.8	return( pow(.5sum/pow((brsum+(5.-wtot)gd->indirect/PI)/5.,.44),
407	greg	1.1	pow((double)n, -.0914) ) );
408			#undef q
409			}
410	greg	2.4
411
412	greg	1.1	#ifndef M_SQRT2
413			#define M_SQRT2 1.41421356237309504880
414			#endif
415	greg	2.4
416	greg	1.1	#define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
417	greg	2.4
418
419	greg	1.1	double
420	greg	1.2	guth_vcp(gd) /* compute Guth visual comfort probability */
421	greg	1.1	struct glare_dir *gd;
422			{
423	greg	2.5	extern double erfc();
424	greg	1.7	double dgr;
425	greg	2.4
426	greg	1.7	dgr = guth_dgr(gd);
427			if (dgr <= FTINY)
428			return(100.0);
429			return(100.norm_integral(6.374-1.3227log(dgr)));
430	greg	1.8	}
431	greg	2.4
432
433	greg	1.8	double
434			cie_cgi(gd) /* compute CIE Glare Index */
435			struct glare_dir *gd;
436			{
437			register struct glare_src *gs;
438			FVECT mydir;
439			double dillum;
440			double p;
441			double sum;
442	greg	2.4
443	greg	1.8	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
444			sum = 0.0;
445			for (gs = all_srcs; gs != NULL; gs = gs->next) {
446			p = posindex(gs->dir, mydir, midview.vup);
447			if (p <= FTINY)
448			continue;
449			sum += gs->lumgs->lum gs->dom / (p*p);
450			}
451			if (sum <= FTINY)
452			return(0.0);
453			dillum = direct(gd);
454			return(8.log10(2.sum*(1.+dillum/500.)/(dillum+gd->indirect)));
455	greg	2.4	}
456
457
458			double
459			ugr(gd) /* compute Unified Glare Rating */
460			struct glare_dir *gd;
461			{
462			register struct glare_src *gs;
463			FVECT mydir;
464			double p;
465			double sum;
466
467			spinvector(mydir, midview.vdir, midview.vup, gd->ang);
468			sum = 0.0;
469			for (gs = all_srcs; gs != NULL; gs = gs->next) {
470			p = posindex(gs->dir, mydir, midview.vup);
471			if (p <= FTINY)
472			continue;
473			sum += gs->lumgs->lum gs->dom / (p*p);
474			}
475			if (sum <= FTINY)
476			return(0.0);
477			return(8.log10(0.25sum*PI/gd->indirect));
478	greg	1.1	}