ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.10
Committed: Fri Jan 2 12:48:36 2004 UTC (20 years, 3 months ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 2.9: +143 -109 lines
Log Message:
Ansification and fixed typing/prototype of getheader() and its callback.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: glarendx.c,v 2.9 2003/11/19 16:21:28 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 "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 register struct glare_src *gs;
191 register 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, "%lf %lf %lf %lf %lf",
209 &gs->dir[0], &gs->dir[1], &gs->dir[2],
210 &gs->dom, &gs->lum) != 5)
211 goto readerr;
212 normalize(gs->dir);
213 gs->next = all_srcs;
214 all_srcs = gs;
215 break;
216 case S_DIREC:
217 if (!strncmp(buf, "END", 3)) {
218 state = S_SEARCH;
219 break;
220 }
221 if ((gd = newp(struct glare_dir)) == NULL)
222 goto memerr;
223 if (sscanf(buf, "%lf %lf",
224 &gd->ang, &gd->indirect) != 2)
225 goto readerr;
226 gd->ang *= PI/180.0; /* convert to radians */
227 gd->next = all_dirs;
228 all_dirs = gd;
229 break;
230 }
231 return;
232 memerr:
233 perror(progname);
234 exit(1);
235 readerr:
236 fprintf(stderr, "%s: read error on input\n", progname);
237 exit(1);
238 #undef S_SEARCH
239 #undef S_SOURCE
240 #undef S_DIREC
241 }
242
243
244 static void
245 print_values( /* print out calculations */
246 gdfun *func
247 )
248 {
249 register struct glare_dir *gd;
250
251 for (gd = all_dirs; gd != NULL; gd = gd->next)
252 printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
253 }
254
255
256 static double
257 direct( /* compute direct vertical illuminance */
258 struct glare_dir *gd
259 )
260 {
261 FVECT mydir;
262 double d, dval;
263 register struct glare_src *gs;
264
265 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
266 dval = 0.0;
267 for (gs = all_srcs; gs != NULL; gs = gs->next) {
268 d = DOT(mydir,gs->dir);
269 if (d > FTINY)
270 dval += d * gs->dom * gs->lum;
271 }
272 return(dval);
273 }
274
275
276 static double
277 indirect( /* return indirect vertical illuminance */
278 struct glare_dir *gd
279 )
280 {
281 return(gd->indirect);
282 }
283
284
285 static double
286 total( /* return total vertical illuminance */
287 struct glare_dir *gd
288 )
289 {
290 return(direct(gd)+gd->indirect);
291 }
292
293
294 /*
295 * posindex - compute glare position index from:
296 *
297 * Source Direction
298 * View Direction
299 * View Up Direction
300 *
301 * All vectors are assumed to be normalized.
302 * This function is an implementation of the method proposed by
303 * Robert Levin in his 1975 JIES article.
304 * This calculation presumes the view direction and up vectors perpendicular.
305 * We return a value less than zero for improper positions.
306 */
307
308 static double
309 posindex( /* compute position index */
310 FVECT sd,
311 FVECT vd,
312 FVECT vu
313 )
314 {
315 double sigma, tau;
316 double d;
317
318 d = DOT(sd,vd);
319 if (d <= 0.0)
320 return(-1.0);
321 if (d >= 1.0)
322 return(1.0);
323 sigma = acos(d) * (180./PI);
324 d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
325 if (d >= 1.0)
326 tau = 0.0;
327 else
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 register struct glare_src *gs;
341 FVECT mydir,testdir[7],vhor;
342 double r,posn,omega,p[7],sum;
343 int i,n;
344
345 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
346 sum = 0.0; n = 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 = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
351 fcross(vhor,gs->dir,midview.vup);
352 normalize(vhor);
353 VCOPY(testdir[0],gs->dir);
354 fvsum(testdir[1],gs->dir,vhor,r);
355 fvsum(testdir[2],gs->dir,vhor,0.5*r);
356 fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
357 fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
358 fvsum(testdir[3],gs->dir,vhor,-r);
359 fvsum(testdir[4],gs->dir,vhor,-0.5*r);
360 fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
361 fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
362 for (i = 0; i < 7; i++) {
363 normalize(testdir[i]);
364 posn = posindex(testdir[i],mydir,midview.vup);
365 if (posn <= FTINY)
366 p[i] = 0.0;
367 else
368 p[i] = 1./(posn*posn);
369 }
370 r = 1-gs->dom/2./PI;
371 omega = gs->dom*p[0];
372 omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
373 omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
374 +0.3334*(p[2]+p[4]+p[5]+p[6]));
375
376 sum += pow(gs->lum,1.6) * pow(omega,0.8) /
377 (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
378 n++;
379 }
380 if (n == 0)
381 return(0.0);
382 return( 10*log10(0.478*sum) );
383 }
384
385
386 static double
387 brs_gi( /* compute BRS Glare Index */
388 struct glare_dir *gd
389 )
390 {
391 register struct glare_src *gs;
392 FVECT mydir;
393 double p;
394 double sum;
395
396 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
397 sum = 0.0;
398 for (gs = all_srcs; gs != NULL; gs = gs->next) {
399 p = posindex(gs->dir, mydir, midview.vup);
400 if (p <= FTINY)
401 continue;
402 sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
403 }
404 if (sum <= FTINY)
405 return(0.0);
406 sum /= gd->indirect/PI;
407 return(10*log10(0.478*sum));
408 }
409
410
411 static double
412 guth_dgr( /* compute Guth discomfort glare rating */
413 struct glare_dir *gd
414 )
415 {
416 #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
417 register struct glare_src *gs;
418 FVECT mydir;
419 double p;
420 double sum;
421 double wtot, brsum;
422 int n;
423
424 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
425 sum = wtot = brsum = 0.0; n = 0;
426 for (gs = all_srcs; gs != NULL; gs = gs->next) {
427 p = posindex(gs->dir, mydir, midview.vup);
428 if (p <= FTINY)
429 continue;
430 sum += gs->lum * q(gs->dom) / p;
431 brsum += gs->lum * gs->dom;
432 wtot += gs->dom;
433 n++;
434 }
435 if (n == 0)
436 return(0.0);
437 return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
438 pow((double)n, -.0914) ) );
439 #undef q
440 }
441
442
443 #ifndef M_SQRT2
444 #define M_SQRT2 1.41421356237309504880
445 #endif
446
447 #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
448
449
450 static double
451 guth_vcp( /* compute Guth visual comfort probability */
452 struct glare_dir *gd
453 )
454 {
455 extern double erfc();
456 double dgr;
457
458 dgr = guth_dgr(gd);
459 if (dgr <= FTINY)
460 return(100.0);
461 return(100.*norm_integral(6.374-1.3227*log(dgr)));
462 }
463
464
465 static double
466 cie_cgi( /* compute CIE Glare Index */
467 struct glare_dir *gd
468 )
469 {
470 register struct glare_src *gs;
471 FVECT mydir;
472 double dillum;
473 double p;
474 double sum;
475
476 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
477 sum = 0.0;
478 for (gs = all_srcs; gs != NULL; gs = gs->next) {
479 p = posindex(gs->dir, mydir, midview.vup);
480 if (p <= FTINY)
481 continue;
482 sum += gs->lum*gs->lum * gs->dom / (p*p);
483 }
484 if (sum <= FTINY)
485 return(0.0);
486 dillum = direct(gd);
487 return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
488 }
489
490
491 static double
492 ugr( /* compute Unified Glare Rating */
493 struct glare_dir *gd
494 )
495 {
496 register struct glare_src *gs;
497 FVECT mydir;
498 double p;
499 double sum;
500
501 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
502 sum = 0.0;
503 for (gs = all_srcs; gs != NULL; gs = gs->next) {
504 p = posindex(gs->dir, mydir, midview.vup);
505 if (p <= FTINY)
506 continue;
507 sum += gs->lum*gs->lum * gs->dom / (p*p);
508 }
509 if (sum <= FTINY)
510 return(0.0);
511 return(8.*log10(0.25*sum*PI/gd->indirect));
512 }