src/util/glarendx.c

/* Copyright (c) 1991 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Compute Glare Index given by program name or -t option:
 *
 *      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
 */

#include "standard.h"
#include "view.h"

extern double   erfc();

double  posindex();

double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
        indirect(), total();

struct named_func {
        char    *name;
        double  (*func)();
        char    *descrip;
} all_funcs[] = {
        {"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[];
{
        extern char     *rindex();
        struct named_func       *funp;
        char    *progtail;
        int     i;
                                        /* get program name */
        progname = argv[0];
        progtail = rindex(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);
}


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");
        }
}


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 = 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
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
}


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;
{
        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)));
}
Revision:	2.3
Committed:	Thu Nov 19 09:59:08 1992 UTC (31 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.2:	+2 -0 lines
Log Message:	added additional protection against acos() domain error
#	User	Rev	Content
1	greg	1.1	/* Copyright (c) 1991 Regents of the University of California */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8	greg	1.9	* Compute Glare Index given by program name or -t option:
9	greg	1.1	*
10	greg	1.2	* guth_dgr - Guth discomfort glare rating
11			* guth_vcp - Guth visual comfort probability
12	greg	1.9	* cie_cgi - CIE Glare Index (1983, due to Einhorn)
13			* vert_dir - Direct vertical illuminance
14			* vert_ind - Indirect vertical illuminance (from input)
15			* vert_ill - Total vertical illuminance
16	greg	1.1	*
17			* 12 April 1991 Greg Ward EPFL
18			*/
19
20			#include "standard.h"
21			#include "view.h"
22
23			extern double erfc();
24
25			double posindex();
26
27	greg	1.9	double direct(), guth_dgr(), guth_vcp(), cie_cgi(),
28			indirect(), total();
29	greg	1.1
30			struct named_func {
31			char *name;
32			double (*func)();
33	greg	1.9	char *descrip;
34	greg	1.1	} all_funcs[] = {
35	greg	1.9	{"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
36			{"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
37			{"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
38			{"vert_dir", direct, "Direct Vertical Illuminance"},
39			{"vert_ill", total, "Total Vertical Illuminance"},
40			{"vert_ind", indirect, "Indirect Vertical Illuminance"},
41	greg	1.1	{NULL}
42			};
43
44			struct glare_src {
45			FVECT dir; /* source direction */
46			double dom; /* solid angle */
47			double lum; /* average luminance */
48			struct glare_src *next;
49			} *all_srcs = NULL;
50
51			struct glare_dir {
52			double ang; /* angle (in radians) */
53			double indirect; /* indirect illuminance */
54			struct glare_dir *next;
55			} *all_dirs = NULL;
56
57			#define newp(type) (type *)malloc(sizeof(type))
58
59			char *progname;
60			int print_header = 1;
61
62			VIEW midview = STDVIEW;
63
64	greg	1.4	int wrongformat = 0;
65	greg	1.1
66	greg	1.4
67	greg	1.1	main(argc, argv)
68			int argc;
69			char *argv[];
70			{
71			extern char *rindex();
72			struct named_func *funp;
73			char *progtail;
74			int i;
75			/* get program name */
76			progname = argv[0];
77			progtail = rindex(progname, '/'); /* final component */
78			if (progtail == NULL)
79			progtail = progname;
80			else
81			progtail++;
82			/* get options */
83			for (i = 1; i < argc && argv[i][0] == '-'; i++)
84			switch (argv[i][1]) {
85			case 't':
86			progtail = argv[++i];
87			break;
88			case 'h':
89			print_header = 0;
90			break;
91			default:
92			goto userr;
93			}
94			if (i < argc-1)
95			goto userr;
96			if (i == argc-1) /* open file */
97			if (freopen(argv[i], "r", stdin) == NULL) {
98			perror(argv[i]);
99			exit(1);
100			}
101	greg	1.9	/* find and run calculation */
102	greg	1.1	for (funp = all_funcs; funp->name != NULL; funp++)
103			if (!strcmp(funp->name, progtail)) {
104	greg	1.9	init();
105			read_input();
106			if (print_header) {
107			printargs(i, argv, stdout);
108			putchar('\n');
109			}
110	greg	1.1	print_values(funp->func);
111			exit(0); /* we're done */
112			}
113			/* invalid function */
114			userr:
115	greg	1.9	fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
116			fprintf(stderr, "\twhere 'type' is one of the following:\n");
117			for (funp = all_funcs; funp->name != NULL; funp++)
118			fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
119	greg	1.1	exit(1);
120			}
121
122
123			headline(s) /* get line from header */
124			char *s;
125			{
126	greg	1.4	char fmt[32];
127
128	greg	1.1	if (print_header) /* copy to output */
129			fputs(s, stdout);
130	greg	2.2	if (isview(s))
131			sscanview(&midview, s);
132	greg	1.4	else if (isformat(s)) {
133			formatval(fmt, s);
134	greg	1.6	wrongformat = strcmp(fmt, "ascii");
135	greg	1.4	}
136	greg	1.1	}
137
138
139	greg	1.9	init() /* initialize calculation */
140			{
141			/* read header */
142			getheader(stdin, headline, NULL);
143			if (wrongformat) {
144			fprintf(stderr, "%s: bad input format\n", progname);
145			exit(1);
146			}
147			/* set view */
148			if (setview(&midview) != NULL) {
149			fprintf(stderr, "%s: bad view information in input\n",
150			progname);
151			exit(1);
152			}
153			}
154
155
156	greg	1.1	read_input() /* read glare sources from stdin */
157			{
158			#define S_SEARCH 0
159			#define S_SOURCE 1
160			#define S_DIREC 2
161			int state = S_SEARCH;
162			char buf[128];
163			register struct glare_src *gs;
164			register struct glare_dir *gd;
165
166			while (fgets(buf, sizeof(buf), stdin) != NULL)
167			switch (state) {
168			case S_SEARCH:
169			if (!strcmp(buf, "BEGIN glare source\n"))
170			state = S_SOURCE;
171			else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
172			state = S_DIREC;
173			break;
174			case S_SOURCE:
175			if (!strncmp(buf, "END", 3)) {
176			state = S_SEARCH;
177			break;
178			}
179			if ((gs = newp(struct glare_src)) == NULL)
180			goto memerr;
181			if (sscanf(buf, "%lf %lf %lf %lf %lf",
182			&gs->dir[0], &gs->dir[1], &gs->dir[2],
183			&gs->dom, &gs->lum) != 5)
184			goto readerr;
185	greg	1.2	normalize(gs->dir);
186	greg	1.1	gs->next = all_srcs;
187			all_srcs = gs;
188			break;
189			case S_DIREC:
190			if (!strncmp(buf, "END", 3)) {
191			state = S_SEARCH;
192			break;
193			}
194			if ((gd = newp(struct glare_dir)) == NULL)
195			goto memerr;
196			if (sscanf(buf, "%lf %lf",
197			&gd->ang, &gd->indirect) != 2)
198			goto readerr;
199			gd->ang = PI/180.0; / convert to radians */
200			gd->next = all_dirs;
201			all_dirs = gd;
202			break;
203			}
204			return;
205			memerr:
206			perror(progname);
207			exit(1);
208			readerr:
209			fprintf(stderr, "%s: read error on input\n", progname);
210			exit(1);
211			#undef S_SEARCH
212			#undef S_SOURCE
213			#undef S_DIREC
214			}
215
216
217			print_values(funp) /* print out calculations */
218			double (*funp)();
219			{
220			register struct glare_dir *gd;
221
222			for (gd = all_dirs; gd != NULL; gd = gd->next)
223			printf("%f\t%f\n", gd->ang(180.0/PI), (funp)(gd));
224			}
225
226
227			double
228	greg	1.9	direct(gd) /* compute direct vertical illuminance */
229	greg	1.1	struct glare_dir *gd;
230			{
231			FVECT mydir;
232			double d, dval;
233			register struct glare_src *gs;
234
235			spinvector(mydir, midview.vdir, midview.vup, gd->ang);
236			dval = 0.0;
237			for (gs = all_srcs; gs != NULL; gs = gs->next) {
238			d = DOT(mydir,gs->dir);
239			if (d > FTINY)
240			dval += d * gs->dom * gs->lum;
241			}
242			return(dval);
243	greg	1.9	}
244
245
246			double
247			indirect(gd) /* return indirect vertical illuminance */
248			struct glare_dir *gd;
249			{
250			return(gd->indirect);
251			}
252
253
254			double
255			total(gd) /* return total vertical illuminance */
256			struct glare_dir *gd;
257			{
258			return(direct(gd)+gd->indirect);
259	greg	1.1	}
260
261
262			/*
263			* posindex - compute glare position index from:
264			*
265			* Source Direction
266			* View Direction
267			* View Up Direction
268			*
269			* All vectors are assumed to be normalized.
270			* This function is an implementation of the method proposed by
271			* Robert Levin in his 1975 JIES article.
272	greg	1.5	* This calculation presumes the view direction and up vectors perpendicular.
273	greg	1.1	* We return a value less than zero for improper positions.
274			*/
275
276			double
277			posindex(sd, vd, vu) /* compute position index */
278			FVECT sd, vd, vu;
279			{
280			double sigma, tau;
281			double d;
282
283			d = DOT(sd,vd);
284			if (d <= 0.0)
285			return(-1.0);
286	greg	1.5	if (d >= 1.0)
287			return(1.0);
288	greg	1.1	sigma = acos(d) * (180./PI);
289	greg	1.5	d = DOT(sd,vu)/sqrt(1.0-d*d);
290			if (d >= 1.0)
291			tau = 0.0;
292	greg	2.3	else if (d <= -1.0)
293			tau = 180.0;
294	greg	1.5	else
295			tau = acos(d) * (180./PI);
296	greg	1.1	return( exp( sigma( (35.2 - tau.31889 - 1.22exp(-.22222tau))*1e-3
297			+ sigma(21. + tau(.26667 + tau-.002963))1e-5 )
298			) );
299			}
300
301
302			double
303	greg	1.2	guth_dgr(gd) /* compute Guth discomfort glare rating */
304	greg	1.1	struct glare_dir *gd;
305			{
306			#define q(w) (20.4w+1.52pow(w,.2)-.075)
307			register struct glare_src *gs;
308	greg	1.7	FVECT mydir;
309	greg	1.1	double p;
310			double sum;
311	greg	1.3	double wtot, brsum;
312	greg	1.1	int n;
313
314	greg	1.7	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
315	greg	1.3	sum = wtot = brsum = 0.0; n = 0;
316	greg	1.1	for (gs = all_srcs; gs != NULL; gs = gs->next) {
317	greg	1.7	p = posindex(gs->dir, mydir, midview.vup);
318	greg	1.1	if (p <= FTINY)
319			continue;
320			sum += gs->lum * q(gs->dom) / p;
321	greg	1.3	brsum += gs->lum * gs->dom;
322			wtot += gs->dom;
323	greg	1.1	n++;
324			}
325			if (n == 0)
326			return(0.0);
327	greg	1.8
328			return( pow(.5sum/pow((brsum+(5.-wtot)gd->indirect/PI)/5.,.44),
329	greg	1.1	pow((double)n, -.0914) ) );
330			#undef q
331			}
332
333
334			extern double erf(), erfc();
335
336			#ifndef M_SQRT2
337			#define M_SQRT2 1.41421356237309504880
338			#endif
339
340			#define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
341
342
343			double
344	greg	1.2	guth_vcp(gd) /* compute Guth visual comfort probability */
345	greg	1.1	struct glare_dir *gd;
346			{
347	greg	1.7	double dgr;
348
349			dgr = guth_dgr(gd);
350			if (dgr <= FTINY)
351			return(100.0);
352			return(100.norm_integral(6.374-1.3227log(dgr)));
353	greg	1.8	}
354
355
356			double
357			cie_cgi(gd) /* compute CIE Glare Index */
358			struct glare_dir *gd;
359			{
360			register struct glare_src *gs;
361			FVECT mydir;
362			double dillum;
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 += gs->lumgs->lum gs->dom / (p*p);
373			}
374			if (sum <= FTINY)
375			return(0.0);
376			dillum = direct(gd);
377			return(8.log10(2.sum*(1.+dillum/500.)/(dillum+gd->indirect)));
378	greg	1.1	}