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

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.16 static const char RCSid[] = "$Id: glarendx.c,v 2.15 2020/02/17 19:35:24 greg Exp $";
3 greg 1.1 #endif
4     /*
5 greg 1.9 * Compute Glare Index given by program name or -t option:
6 greg 1.1 *
7 greg 2.4 * dgi - Daylight Glare Index
8     * brs_gi - Building Research Station Glare Index (Petherbridge
9     * & Hopkinson)
10     * ugr - Unified Glare Rating System (Fischer)
11 greg 1.2 * guth_dgr - Guth discomfort glare rating
12     * guth_vcp - Guth visual comfort probability
13 greg 1.9 * 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 greg 1.1 *
18     * 12 April 1991 Greg Ward EPFL
19 greg 2.4 * 19 April 1993 R. Compagnon EPFL (added dgi, brs_gi, ugr)
20 greg 1.1 */
21 schorsch 2.10
22 schorsch 2.8 #include <string.h>
23    
24 greg 1.1 #include "standard.h"
25 greg 2.16 #include "paths.h"
26 greg 1.1 #include "view.h"
27 schorsch 2.10
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 greg 1.1 struct named_func {
62     char *name;
63 schorsch 2.10 gdfun *func;
64 greg 1.9 char *descrip;
65 greg 1.1 } all_funcs[] = {
66 greg 2.4 {"dgi", dgi, "Daylight Glare Index"},
67     {"brs_gi", brs_gi, "BRS Glare Index"},
68     {"ugr", ugr, "Unified Glare Rating"},
69 greg 1.9 {"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 greg 1.1 {NULL}
76     };
77 schorsch 2.10
78 greg 1.1 #define newp(type) (type *)malloc(sizeof(type))
79 schorsch 2.10
80 greg 1.1 char *progname;
81     int print_header = 1;
82 schorsch 2.10
83 greg 1.1 VIEW midview = STDVIEW;
84 schorsch 2.10
85 greg 1.4 int wrongformat = 0;
86 schorsch 2.10
87    
88     int
89     main(
90     int argc,
91     char *argv[]
92     )
93 greg 1.1 {
94     struct named_func *funp;
95     int i;
96     /* get program name */
97 greg 2.16 progname = fixargv0(argv[0]);
98 greg 1.1 /* get options */
99     for (i = 1; i < argc && argv[i][0] == '-'; i++)
100     switch (argv[i][1]) {
101     case 't':
102 greg 2.16 progname = argv[++i];
103 greg 1.1 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 greg 1.9 /* find and run calculation */
118 greg 1.1 for (funp = all_funcs; funp->name != NULL; funp++)
119 greg 2.16 if (!strcmp(funp->name, progname)) {
120 greg 1.9 init();
121     read_input();
122     if (print_header) {
123     printargs(i, argv, stdout);
124     putchar('\n');
125     }
126 greg 1.1 print_values(funp->func);
127     exit(0); /* we're done */
128     }
129     /* invalid function */
130     userr:
131 greg 1.9 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 greg 1.1 exit(1);
136     }
137 schorsch 2.10
138    
139     static int
140     headline( /* get line from header */
141     char *s,
142     void *p
143     )
144 greg 1.1 {
145 greg 2.13 char fmt[MAXFMTLEN];
146 schorsch 2.10
147 greg 1.1 if (print_header) /* copy to output */
148     fputs(s, stdout);
149 greg 2.2 if (isview(s))
150     sscanview(&midview, s);
151 greg 1.4 else if (isformat(s)) {
152     formatval(fmt, s);
153 greg 1.6 wrongformat = strcmp(fmt, "ascii");
154 greg 1.4 }
155 gwlarson 2.6 return(0);
156 greg 1.1 }
157 schorsch 2.10
158    
159     static void
160     init(void) /* initialize calculation */
161 greg 1.9 {
162 greg 2.14 char *err;
163 greg 1.9 /* 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 greg 2.14 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 greg 1.9 progname);
179     exit(1);
180     }
181     }
182 schorsch 2.10
183    
184     static void
185     read_input(void) /* read glare sources from stdin */
186 greg 1.1 {
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 greg 2.11 struct glare_src *gs;
193     struct glare_dir *gd;
194 schorsch 2.10
195 greg 1.1 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 greg 2.11 if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
211     &gs->dir[2]) != 3 ||
212 greg 2.12 sscanf(sskip2(buf, 3), "%lf %lf",
213 greg 2.11 &gs->dom, &gs->lum) != 2)
214 greg 1.1 goto readerr;
215 greg 1.2 normalize(gs->dir);
216 greg 1.1 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 schorsch 2.10
246    
247     static void
248     print_values( /* print out calculations */
249     gdfun *func
250     )
251 greg 1.1 {
252 greg 2.11 struct glare_dir *gd;
253 schorsch 2.10
254 greg 1.1 for (gd = all_dirs; gd != NULL; gd = gd->next)
255 schorsch 2.10 printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
256 greg 1.1 }
257 schorsch 2.10
258    
259     static double
260     direct( /* compute direct vertical illuminance */
261     struct glare_dir *gd
262     )
263 greg 1.1 {
264     FVECT mydir;
265     double d, dval;
266 greg 2.11 struct glare_src *gs;
267 schorsch 2.10
268 greg 1.1 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 greg 1.9 }
277 schorsch 2.10
278    
279     static double
280     indirect( /* return indirect vertical illuminance */
281     struct glare_dir *gd
282     )
283 greg 1.9 {
284     return(gd->indirect);
285     }
286 schorsch 2.10
287    
288     static double
289     total( /* return total vertical illuminance */
290     struct glare_dir *gd
291     )
292 greg 1.9 {
293     return(direct(gd)+gd->indirect);
294 greg 1.1 }
295 schorsch 2.10
296    
297 greg 1.1 /*
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 greg 1.5 * This calculation presumes the view direction and up vectors perpendicular.
308 greg 1.1 * We return a value less than zero for improper positions.
309     */
310 schorsch 2.10
311     static double
312     posindex( /* compute position index */
313     FVECT sd,
314     FVECT vd,
315     FVECT vu
316     )
317 greg 1.1 {
318     double sigma, tau;
319     double d;
320 schorsch 2.10
321 greg 1.1 d = DOT(sd,vd);
322     if (d <= 0.0)
323     return(-1.0);
324 greg 2.14 if (d >= 1.0-FTINY)
325 greg 1.5 return(1.0);
326 greg 1.1 sigma = acos(d) * (180./PI);
327 greg 2.4 d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
328 greg 2.14 tau = Acos(d) * (180./PI);
329 greg 1.1 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 schorsch 2.10
334    
335     static double
336     dgi( /* compute Daylight Glare Index */
337     struct glare_dir *gd
338     )
339 greg 2.4 {
340 greg 2.11 struct glare_src *gs;
341 greg 2.4 FVECT mydir,testdir[7],vhor;
342 greg 2.9 double r,posn,omega,p[7],sum;
343 greg 2.14 int i;
344 schorsch 2.10
345 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
346 greg 2.14 sum = 0.0;
347 greg 2.4 for (gs = all_srcs; gs != NULL; gs = gs->next) {
348 greg 1.1
349 greg 2.4 /* compute 1/p^2 weighted solid angle of the source */
350 greg 2.14 r = 1. - gs->dom*(.5/PI);
351     r = sqrt(1. - r*r);
352 greg 2.4 fcross(vhor,gs->dir,midview.vup);
353 greg 2.14 if (normalize(vhor) == 0.)
354     continue;
355 greg 2.4 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 greg 2.14 p[i] = 0.0;
366     if (normalize(testdir[i]) == 0.)
367     continue;
368 greg 2.9 posn = posindex(testdir[i],mydir,midview.vup);
369 greg 2.14 if (posn > FTINY)
370 greg 2.9 p[i] = 1./(posn*posn);
371 greg 2.4 }
372 greg 2.14 r = 1. - gs->dom*(.5/PI);
373 greg 2.4 omega = gs->dom*p[0];
374 greg 2.14 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 greg 2.4 +0.3334*(p[2]+p[4]+p[5]+p[6]));
377 greg 2.14 if (omega <= 0.)
378 greg 2.15 continue;
379 schorsch 2.10 sum += pow(gs->lum,1.6) * pow(omega,0.8) /
380 greg 2.14 (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
381 greg 2.4 }
382 greg 2.14 if (sum <= FTINY)
383 greg 2.4 return(0.0);
384 greg 2.14 return( 10.*log10(0.478*sum) );
385 greg 2.4 }
386    
387 schorsch 2.10
388     static double
389     brs_gi( /* compute BRS Glare Index */
390     struct glare_dir *gd
391     )
392 greg 2.4 {
393 greg 2.11 struct glare_src *gs;
394 greg 2.4 FVECT mydir;
395     double p;
396     double sum;
397 schorsch 2.10
398 greg 2.4 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 greg 2.14 sum *= PI/gd->indirect;
409 greg 2.4 return(10*log10(0.478*sum));
410     }
411    
412    
413 schorsch 2.10 static double
414     guth_dgr( /* compute Guth discomfort glare rating */
415     struct glare_dir *gd
416     )
417 greg 1.1 {
418     #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
419 greg 2.11 struct glare_src *gs;
420 greg 1.7 FVECT mydir;
421 greg 1.1 double p;
422     double sum;
423 greg 1.3 double wtot, brsum;
424 greg 1.1 int n;
425 schorsch 2.10
426 greg 1.7 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
427 greg 1.3 sum = wtot = brsum = 0.0; n = 0;
428 greg 1.1 for (gs = all_srcs; gs != NULL; gs = gs->next) {
429 greg 1.7 p = posindex(gs->dir, mydir, midview.vup);
430 greg 1.1 if (p <= FTINY)
431     continue;
432     sum += gs->lum * q(gs->dom) / p;
433 greg 1.3 brsum += gs->lum * gs->dom;
434     wtot += gs->dom;
435 greg 1.1 n++;
436     }
437     if (n == 0)
438     return(0.0);
439 greg 1.8 return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
440 greg 1.1 pow((double)n, -.0914) ) );
441     #undef q
442     }
443 schorsch 2.10
444    
445 greg 1.1 #ifndef M_SQRT2
446     #define M_SQRT2 1.41421356237309504880
447     #endif
448 schorsch 2.10
449 greg 1.1 #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
450 schorsch 2.10
451    
452     static double
453     guth_vcp( /* compute Guth visual comfort probability */
454     struct glare_dir *gd
455     )
456 greg 1.1 {
457 greg 2.5 extern double erfc();
458 greg 1.7 double dgr;
459 schorsch 2.10
460 greg 1.7 dgr = guth_dgr(gd);
461     if (dgr <= FTINY)
462     return(100.0);
463     return(100.*norm_integral(6.374-1.3227*log(dgr)));
464 greg 1.8 }
465 schorsch 2.10
466    
467     static double
468     cie_cgi( /* compute CIE Glare Index */
469     struct glare_dir *gd
470     )
471 greg 1.8 {
472 greg 2.11 struct glare_src *gs;
473 greg 1.8 FVECT mydir;
474     double dillum;
475     double p;
476     double sum;
477 schorsch 2.10
478 greg 1.8 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 greg 2.4 }
491    
492    
493 schorsch 2.10 static double
494     ugr( /* compute Unified Glare Rating */
495     struct glare_dir *gd
496     )
497 greg 2.4 {
498 greg 2.11 struct glare_src *gs;
499 greg 2.4 FVECT mydir;
500     double p;
501     double sum;
502 schorsch 2.10
503 greg 2.4 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 greg 1.1 }