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:
 *
 *      gnuth_dgr -     Gnuth discomfort glare rating
 *      gnuth_vcp -     Gnuth visual comfort probability
 *
 *              12 April 1991   Greg Ward       EPFL
 */

#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 */
        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;


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);
                }
                                        /* 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 */
        for (funp = all_funcs; funp->name != NULL; funp++)
                if (!strcmp(funp->name, progtail)) {
                        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);
        exit(1);
}


headline(s)                     /* get line from header */
char    *s;
{
        if (print_header)               /* copy to output */
                fputs(s, stdout);
        if (!strncmp(s, VIEWSTR, VIEWSTRL))
                sscanview(&midview, s+VIEWSTRL);
}


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


/*
 * 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.
 * 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);
        sigma = acos(d) * (180./PI);
        tau = acos(DOT(sd,vu)/sqrt(1.0-d*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;
{
#define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
        register struct glare_src       *gs;
        double  p;
        double  sum;
        int     n;

        sum = 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;
                n++;
        }
        if (n == 0)
                return(0.0);
        else
                return( pow(
                        .5*sum/pow(direct(gd)+gd->indirect,.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
gnuth_vcp(gd)           /* compute Gnuth visual comfort probability */
struct glare_dir        *gd;
{
        return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd))));
}
Revision:	1.1
Committed:	Tue Apr 16 16:40:22 1991 UTC (33 years ago) by greg
Content type:	text/plain
Branch:	MAIN
Log Message:	Initial revision
#	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			* Compute Glare Index given by program name:
9			*
10			* gnuth_dgr - Gnuth discomfort glare rating
11			* gnuth_vcp - Gnuth visual comfort probability
12			*
13			* 12 April 1991 Greg Ward EPFL
14			*/
15
16			#include "standard.h"
17			#include "view.h"
18
19			extern double erfc();
20
21			double posindex();
22			int headline();
23
24			double direct(), gnuth_dgr(), gnuth_vcp();
25
26			struct named_func {
27			char *name;
28			double (*func)();
29			} all_funcs[] = {
30			{"direct", direct},
31			{"gnuth_dgr", gnuth_dgr},
32			{"gnuth_vcp", gnuth_vcp},
33			{NULL}
34			};
35
36			struct glare_src {
37			FVECT dir; /* source direction */
38			double dom; /* solid angle */
39			double lum; /* average luminance */
40			struct glare_src *next;
41			} *all_srcs = NULL;
42
43			struct glare_dir {
44			double ang; /* angle (in radians) */
45			double indirect; /* indirect illuminance */
46			struct glare_dir *next;
47			} *all_dirs = NULL;
48
49			#define newp(type) (type *)malloc(sizeof(type))
50
51			char *progname;
52			int print_header = 1;
53
54			VIEW midview = STDVIEW;
55
56
57			main(argc, argv)
58			int argc;
59			char *argv[];
60			{
61			extern char *rindex();
62			struct named_func *funp;
63			char *progtail;
64			int i;
65			/* get program name */
66			progname = argv[0];
67			progtail = rindex(progname, '/'); /* final component */
68			if (progtail == NULL)
69			progtail = progname;
70			else
71			progtail++;
72			/* get options */
73			for (i = 1; i < argc && argv[i][0] == '-'; i++)
74			switch (argv[i][1]) {
75			case 't':
76			progtail = argv[++i];
77			break;
78			case 'h':
79			print_header = 0;
80			break;
81			default:
82			goto userr;
83			}
84			if (i < argc-1)
85			goto userr;
86			if (i == argc-1) /* open file */
87			if (freopen(argv[i], "r", stdin) == NULL) {
88			perror(argv[i]);
89			exit(1);
90			}
91			/* read header */
92			getheader(stdin, headline);
93			if (print_header) { /* add to header */
94			printargs(i, argv, stdout);
95			putchar('\n');
96			}
97			/* set view */
98			if (setview(&midview) != NULL) {
99			fprintf(stderr, "%s: bad view information in input\n");
100			exit(1);
101			}
102			/* get findglare data */
103			read_input();
104			/* find calculation */
105			for (funp = all_funcs; funp->name != NULL; funp++)
106			if (!strcmp(funp->name, progtail)) {
107			print_values(funp->func);
108			exit(0); /* we're done */
109			}
110			/* invalid function */
111			fprintf(stderr, "%s: unknown function!\n", progtail);
112			exit(1);
113			userr:
114			fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
115			exit(1);
116			}
117
118
119			headline(s) /* get line from header */
120			char *s;
121			{
122			if (print_header) /* copy to output */
123			fputs(s, stdout);
124			if (!strncmp(s, VIEWSTR, VIEWSTRL))
125			sscanview(&midview, s+VIEWSTRL);
126			}
127
128
129			read_input() /* read glare sources from stdin */
130			{
131			#define S_SEARCH 0
132			#define S_SOURCE 1
133			#define S_DIREC 2
134			int state = S_SEARCH;
135			char buf[128];
136			register struct glare_src *gs;
137			register struct glare_dir *gd;
138
139			while (fgets(buf, sizeof(buf), stdin) != NULL)
140			switch (state) {
141			case S_SEARCH:
142			if (!strcmp(buf, "BEGIN glare source\n"))
143			state = S_SOURCE;
144			else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
145			state = S_DIREC;
146			break;
147			case S_SOURCE:
148			if (!strncmp(buf, "END", 3)) {
149			state = S_SEARCH;
150			break;
151			}
152			if ((gs = newp(struct glare_src)) == NULL)
153			goto memerr;
154			if (sscanf(buf, "%lf %lf %lf %lf %lf",
155			&gs->dir[0], &gs->dir[1], &gs->dir[2],
156			&gs->dom, &gs->lum) != 5)
157			goto readerr;
158			gs->next = all_srcs;
159			all_srcs = gs;
160			break;
161			case S_DIREC:
162			if (!strncmp(buf, "END", 3)) {
163			state = S_SEARCH;
164			break;
165			}
166			if ((gd = newp(struct glare_dir)) == NULL)
167			goto memerr;
168			if (sscanf(buf, "%lf %lf",
169			&gd->ang, &gd->indirect) != 2)
170			goto readerr;
171			gd->ang = PI/180.0; / convert to radians */
172			gd->next = all_dirs;
173			all_dirs = gd;
174			break;
175			}
176			return;
177			memerr:
178			perror(progname);
179			exit(1);
180			readerr:
181			fprintf(stderr, "%s: read error on input\n", progname);
182			exit(1);
183			#undef S_SEARCH
184			#undef S_SOURCE
185			#undef S_DIREC
186			}
187
188
189			print_values(funp) /* print out calculations */
190			double (*funp)();
191			{
192			register struct glare_dir *gd;
193
194			for (gd = all_dirs; gd != NULL; gd = gd->next)
195			printf("%f\t%f\n", gd->ang(180.0/PI), (funp)(gd));
196			}
197
198
199			double
200			direct(gd) /* compute direct illuminance */
201			struct glare_dir *gd;
202			{
203			FVECT mydir;
204			double d, dval;
205			register struct glare_src *gs;
206
207			spinvector(mydir, midview.vdir, midview.vup, gd->ang);
208			dval = 0.0;
209			for (gs = all_srcs; gs != NULL; gs = gs->next) {
210			d = DOT(mydir,gs->dir);
211			if (d > FTINY)
212			dval += d * gs->dom * gs->lum;
213			}
214			return(dval);
215			}
216
217
218			/*
219			* posindex - compute glare position index from:
220			*
221			* Source Direction
222			* View Direction
223			* View Up Direction
224			*
225			* All vectors are assumed to be normalized.
226			* This function is an implementation of the method proposed by
227			* Robert Levin in his 1975 JIES article.
228			* We return a value less than zero for improper positions.
229			*/
230
231			double
232			posindex(sd, vd, vu) /* compute position index */
233			FVECT sd, vd, vu;
234			{
235			double sigma, tau;
236			double d;
237
238			d = DOT(sd,vd);
239			if (d <= 0.0)
240			return(-1.0);
241			sigma = acos(d) * (180./PI);
242			tau = acos(DOT(sd,vu)/sqrt(1.0-dd)) (180./PI);
243			return( exp( sigma( (35.2 - tau.31889 - 1.22exp(-.22222tau))*1e-3
244			+ sigma(21. + tau(.26667 + tau-.002963))1e-5 )
245			) );
246			}
247
248
249			double
250			gnuth_dgr(gd) /* compute Gnuth discomfort glare rating */
251			struct glare_dir *gd;
252			{
253			#define q(w) (20.4w+1.52pow(w,.2)-.075)
254			register struct glare_src *gs;
255			double p;
256			double sum;
257			int n;
258
259			sum = 0.0; n = 0;
260			for (gs = all_srcs; gs != NULL; gs = gs->next) {
261			p = posindex(gs->dir, midview.vdir, midview.vup);
262			if (p <= FTINY)
263			continue;
264			sum += gs->lum * q(gs->dom) / p;
265			n++;
266			}
267			if (n == 0)
268			return(0.0);
269			else
270			return( pow(
271			.5*sum/pow(direct(gd)+gd->indirect,.44),
272			pow((double)n, -.0914) ) );
273			#undef q
274			}
275
276
277			extern double erf(), erfc();
278
279			#ifndef M_SQRT2
280			#define M_SQRT2 1.41421356237309504880
281			#endif
282
283			#define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
284
285
286			double
287			gnuth_vcp(gd) /* compute Gnuth visual comfort probability */
288			struct glare_dir *gd;
289			{
290			return(100.norm_integral(-6.374+1.3227log(gnuth_dgr(gd))));
291			}