ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.16
Committed: Fri Jul 24 16:58:16 2020 UTC (3 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R3, HEAD
Changes since 2.15: +5 -10 lines
Log Message:
Hopeful fix to Windows test issues

File Contents

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