ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.9
Committed: Wed Nov 19 16:21:28 2003 UTC (20 years, 4 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

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: glarendx.c,v 2.8 2003/06/30 14:59:13 schorsch Exp $";
3 #endif
4 /*
5 * Compute Glare Index given by program name or -t option:
6 *
7 * dgi - Daylight Glare Index
8 * brs_gi - Building Research Station Glare Index (Petherbridge
9 * & Hopkinson)
10 * ugr - Unified Glare Rating System (Fischer)
11 * guth_dgr - Guth discomfort glare rating
12 * guth_vcp - Guth visual comfort probability
13 * 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 *
18 * 12 April 1991 Greg Ward EPFL
19 * 19 April 1993 R. Compagnon EPFL (added dgi, brs_gi, ugr)
20 */
21
22 #include <string.h>
23
24 #include "standard.h"
25 #include "view.h"
26
27
28 double posindex();
29
30 double direct(), guth_dgr(), guth_vcp(), cie_cgi(),
31 indirect(), total(), dgi(), brs_gi(), ugr();
32
33 struct named_func {
34 char *name;
35 double (*func)();
36 char *descrip;
37 } all_funcs[] = {
38 {"dgi", dgi, "Daylight Glare Index"},
39 {"brs_gi", brs_gi, "BRS Glare Index"},
40 {"ugr", ugr, "Unified Glare Rating"},
41 {"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 {NULL}
48 };
49
50 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
57 struct glare_dir {
58 double ang; /* angle (in radians) */
59 double indirect; /* indirect illuminance */
60 struct glare_dir *next;
61 } *all_dirs = NULL;
62
63 #define newp(type) (type *)malloc(sizeof(type))
64
65 char *progname;
66 int print_header = 1;
67
68 VIEW midview = STDVIEW;
69
70 int wrongformat = 0;
71
72
73 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 progtail = strrchr(progname, '/'); /* final component */
83 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 /* find and run calculation */
107 for (funp = all_funcs; funp->name != NULL; funp++)
108 if (!strcmp(funp->name, progtail)) {
109 init();
110 read_input();
111 if (print_header) {
112 printargs(i, argv, stdout);
113 putchar('\n');
114 }
115 print_values(funp->func);
116 exit(0); /* we're done */
117 }
118 /* invalid function */
119 userr:
120 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 exit(1);
125 }
126
127
128 int
129 headline(s) /* get line from header */
130 char *s;
131 {
132 char fmt[32];
133
134 if (print_header) /* copy to output */
135 fputs(s, stdout);
136 if (isview(s))
137 sscanview(&midview, s);
138 else if (isformat(s)) {
139 formatval(fmt, s);
140 wrongformat = strcmp(fmt, "ascii");
141 }
142 return(0);
143 }
144
145
146 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
162
163 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
173 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 normalize(gs->dir);
193 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
223
224 print_values(funp) /* print out calculations */
225 double (*funp)();
226 {
227 register struct glare_dir *gd;
228
229 for (gd = all_dirs; gd != NULL; gd = gd->next)
230 printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
231 }
232
233
234 double
235 direct(gd) /* compute direct vertical illuminance */
236 struct glare_dir *gd;
237 {
238 FVECT mydir;
239 double d, dval;
240 register struct glare_src *gs;
241
242 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 }
251
252
253 double
254 indirect(gd) /* return indirect vertical illuminance */
255 struct glare_dir *gd;
256 {
257 return(gd->indirect);
258 }
259
260
261 double
262 total(gd) /* return total vertical illuminance */
263 struct glare_dir *gd;
264 {
265 return(direct(gd)+gd->indirect);
266 }
267
268
269 /*
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 * This calculation presumes the view direction and up vectors perpendicular.
280 * We return a value less than zero for improper positions.
281 */
282
283 double
284 posindex(sd, vd, vu) /* compute position index */
285 FVECT sd, vd, vu;
286 {
287 double sigma, tau;
288 double d;
289
290 d = DOT(sd,vd);
291 if (d <= 0.0)
292 return(-1.0);
293 if (d >= 1.0)
294 return(1.0);
295 sigma = acos(d) * (180./PI);
296 d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
297 if (d >= 1.0)
298 tau = 0.0;
299 else
300 tau = acos(d) * (180./PI);
301 return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
302 + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
303 ) );
304 }
305
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 double r,posn,omega,p[7],sum;
314 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
320 /* 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 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 }
341 r = 1-gs->dom/2./PI;
342 omega = gs->dom*p[0];
343 omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
344 omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
345 +0.3334*(p[2]+p[4]+p[5]+p[6]));
346
347 sum += pow(gs->lum,1.6) * pow(omega,0.8) /
348 (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
349 n++;
350 }
351 if (n == 0)
352 return(0.0);
353 return( 10*log10(0.478*sum) );
354 }
355
356
357 double
358 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(10*log10(0.478*sum));
378 }
379
380
381 double
382 guth_dgr(gd) /* compute Guth discomfort glare rating */
383 struct glare_dir *gd;
384 {
385 #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
386 register struct glare_src *gs;
387 FVECT mydir;
388 double p;
389 double sum;
390 double wtot, brsum;
391 int n;
392
393 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
394 sum = wtot = brsum = 0.0; n = 0;
395 for (gs = all_srcs; gs != NULL; gs = gs->next) {
396 p = posindex(gs->dir, mydir, midview.vup);
397 if (p <= FTINY)
398 continue;
399 sum += gs->lum * q(gs->dom) / p;
400 brsum += gs->lum * gs->dom;
401 wtot += gs->dom;
402 n++;
403 }
404 if (n == 0)
405 return(0.0);
406 return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
407 pow((double)n, -.0914) ) );
408 #undef q
409 }
410
411
412 #ifndef M_SQRT2
413 #define M_SQRT2 1.41421356237309504880
414 #endif
415
416 #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
417
418
419 double
420 guth_vcp(gd) /* compute Guth visual comfort probability */
421 struct glare_dir *gd;
422 {
423 extern double erfc();
424 double dgr;
425
426 dgr = guth_dgr(gd);
427 if (dgr <= FTINY)
428 return(100.0);
429 return(100.*norm_integral(6.374-1.3227*log(dgr)));
430 }
431
432
433 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
443 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->lum*gs->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 }
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->lum*gs->lum * gs->dom / (p*p);
474 }
475 if (sum <= FTINY)
476 return(0.0);
477 return(8.*log10(0.25*sum*PI/gd->indirect));
478 }