ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.11
Committed: Fri Apr 11 20:27:23 2014 UTC (9 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R2P1
Changes since 2.10: +14 -13 lines
Log Message:
Made reading of vectors more consistent with FVFORMAT macro

File Contents

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