ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.8
Committed: Mon Jun 30 14:59:13 2003 UTC (20 years, 9 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.7: +4 -3 lines
Log Message:
Replaced most outdated BSD function calls with their posix equivalents, and cleaned up a few other platform dependencies.

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 schorsch 2.8 static const char RCSid[] = "$Id: glarendx.c,v 2.7 2003/02/22 02:07:30 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 greg 2.4
22 schorsch 2.8 #include <string.h>
23    
24 greg 1.1 #include "standard.h"
25     #include "view.h"
26 greg 2.4
27    
28 greg 1.1 double posindex();
29 greg 2.4
30 greg 1.9 double direct(), guth_dgr(), guth_vcp(), cie_cgi(),
31 greg 2.4 indirect(), total(), dgi(), brs_gi(), ugr();
32    
33 greg 1.1 struct named_func {
34     char *name;
35     double (*func)();
36 greg 1.9 char *descrip;
37 greg 1.1 } all_funcs[] = {
38 greg 2.4 {"dgi", dgi, "Daylight Glare Index"},
39     {"brs_gi", brs_gi, "BRS Glare Index"},
40     {"ugr", ugr, "Unified Glare Rating"},
41 greg 1.9 {"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 greg 1.1 {NULL}
48     };
49 greg 2.4
50 greg 1.1 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 greg 2.4
57 greg 1.1 struct glare_dir {
58     double ang; /* angle (in radians) */
59     double indirect; /* indirect illuminance */
60     struct glare_dir *next;
61     } *all_dirs = NULL;
62 greg 2.4
63 greg 1.1 #define newp(type) (type *)malloc(sizeof(type))
64 greg 2.4
65 greg 1.1 char *progname;
66     int print_header = 1;
67 greg 2.4
68 greg 1.1 VIEW midview = STDVIEW;
69 greg 2.4
70 greg 1.4 int wrongformat = 0;
71 greg 2.4
72    
73 greg 1.1 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 schorsch 2.8 progtail = strrchr(progname, '/'); /* final component */
83 greg 1.1 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 greg 1.9 /* find and run calculation */
107 greg 1.1 for (funp = all_funcs; funp->name != NULL; funp++)
108     if (!strcmp(funp->name, progtail)) {
109 greg 1.9 init();
110     read_input();
111     if (print_header) {
112     printargs(i, argv, stdout);
113     putchar('\n');
114     }
115 greg 1.1 print_values(funp->func);
116     exit(0); /* we're done */
117     }
118     /* invalid function */
119     userr:
120 greg 1.9 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 greg 1.1 exit(1);
125     }
126 greg 2.4
127    
128 gwlarson 2.6 int
129 greg 1.1 headline(s) /* get line from header */
130     char *s;
131     {
132 greg 1.4 char fmt[32];
133 greg 2.4
134 greg 1.1 if (print_header) /* copy to output */
135     fputs(s, stdout);
136 greg 2.2 if (isview(s))
137     sscanview(&midview, s);
138 greg 1.4 else if (isformat(s)) {
139     formatval(fmt, s);
140 greg 1.6 wrongformat = strcmp(fmt, "ascii");
141 greg 1.4 }
142 gwlarson 2.6 return(0);
143 greg 1.1 }
144 greg 2.4
145    
146 greg 1.9 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 greg 2.4
162    
163 greg 1.1 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 greg 2.4
173 greg 1.1 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 greg 1.2 normalize(gs->dir);
193 greg 1.1 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 greg 2.4
223    
224 greg 1.1 print_values(funp) /* print out calculations */
225     double (*funp)();
226     {
227     register struct glare_dir *gd;
228 greg 2.4
229 greg 1.1 for (gd = all_dirs; gd != NULL; gd = gd->next)
230     printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
231     }
232 greg 2.4
233    
234 greg 1.1 double
235 greg 1.9 direct(gd) /* compute direct vertical illuminance */
236 greg 1.1 struct glare_dir *gd;
237     {
238     FVECT mydir;
239     double d, dval;
240     register struct glare_src *gs;
241 greg 2.4
242 greg 1.1 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 greg 1.9 }
251 greg 2.4
252    
253 greg 1.9 double
254     indirect(gd) /* return indirect vertical illuminance */
255     struct glare_dir *gd;
256     {
257     return(gd->indirect);
258     }
259 greg 2.4
260    
261 greg 1.9 double
262     total(gd) /* return total vertical illuminance */
263     struct glare_dir *gd;
264     {
265     return(direct(gd)+gd->indirect);
266 greg 1.1 }
267 greg 2.4
268    
269 greg 1.1 /*
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 greg 1.5 * This calculation presumes the view direction and up vectors perpendicular.
280 greg 1.1 * We return a value less than zero for improper positions.
281     */
282 greg 2.4
283 greg 1.1 double
284     posindex(sd, vd, vu) /* compute position index */
285     FVECT sd, vd, vu;
286     {
287     double sigma, tau;
288     double d;
289 greg 2.4
290 greg 1.1 d = DOT(sd,vd);
291     if (d <= 0.0)
292     return(-1.0);
293 greg 1.5 if (d >= 1.0)
294     return(1.0);
295 greg 1.1 sigma = acos(d) * (180./PI);
296 greg 2.4 d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
297 greg 1.5 if (d >= 1.0)
298     tau = 0.0;
299     else
300     tau = acos(d) * (180./PI);
301 greg 1.1 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 greg 2.4
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,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 greg 1.1
320 greg 2.4 /* 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     p[i] = pow(posindex(testdir[i],mydir,midview.vup),-2.0);
336     if (p[i] <= FTINY) p[i] = 0.0;
337     }
338     r = 1-gs->dom/2./PI;
339     omega = gs->dom*p[0];
340     omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
341     omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
342     +0.3334*(p[2]+p[4]+p[5]+p[6]));
343 greg 1.1
344 greg 2.4 sum += pow(gs->lum,1.6) * pow(omega,0.8) /
345     (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
346     n++;
347     }
348     if (n == 0)
349     return(0.0);
350     return( 10*log10(0.478*sum) );
351     }
352    
353    
354 greg 1.1 double
355 greg 2.4 brs_gi(gd) /* compute BRS Glare Index */
356     struct glare_dir *gd;
357     {
358     register struct glare_src *gs;
359     FVECT mydir;
360     double p;
361     double sum;
362    
363     spinvector(mydir, midview.vdir, midview.vup, gd->ang);
364     sum = 0.0;
365     for (gs = all_srcs; gs != NULL; gs = gs->next) {
366     p = posindex(gs->dir, mydir, midview.vup);
367     if (p <= FTINY)
368     continue;
369     sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
370     }
371     if (sum <= FTINY)
372     return(0.0);
373     sum /= gd->indirect/PI;
374     return(10*log10(0.478*sum));
375     }
376    
377    
378     double
379 greg 1.2 guth_dgr(gd) /* compute Guth discomfort glare rating */
380 greg 1.1 struct glare_dir *gd;
381     {
382     #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
383     register struct glare_src *gs;
384 greg 1.7 FVECT mydir;
385 greg 1.1 double p;
386     double sum;
387 greg 1.3 double wtot, brsum;
388 greg 1.1 int n;
389 greg 2.4
390 greg 1.7 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
391 greg 1.3 sum = wtot = brsum = 0.0; n = 0;
392 greg 1.1 for (gs = all_srcs; gs != NULL; gs = gs->next) {
393 greg 1.7 p = posindex(gs->dir, mydir, midview.vup);
394 greg 1.1 if (p <= FTINY)
395     continue;
396     sum += gs->lum * q(gs->dom) / p;
397 greg 1.3 brsum += gs->lum * gs->dom;
398     wtot += gs->dom;
399 greg 1.1 n++;
400     }
401     if (n == 0)
402     return(0.0);
403 greg 1.8 return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
404 greg 1.1 pow((double)n, -.0914) ) );
405     #undef q
406     }
407 greg 2.4
408    
409 greg 1.1 #ifndef M_SQRT2
410     #define M_SQRT2 1.41421356237309504880
411     #endif
412 greg 2.4
413 greg 1.1 #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
414 greg 2.4
415    
416 greg 1.1 double
417 greg 1.2 guth_vcp(gd) /* compute Guth visual comfort probability */
418 greg 1.1 struct glare_dir *gd;
419     {
420 greg 2.5 extern double erfc();
421 greg 1.7 double dgr;
422 greg 2.4
423 greg 1.7 dgr = guth_dgr(gd);
424     if (dgr <= FTINY)
425     return(100.0);
426     return(100.*norm_integral(6.374-1.3227*log(dgr)));
427 greg 1.8 }
428 greg 2.4
429    
430 greg 1.8 double
431     cie_cgi(gd) /* compute CIE Glare Index */
432     struct glare_dir *gd;
433     {
434     register struct glare_src *gs;
435     FVECT mydir;
436     double dillum;
437     double p;
438     double sum;
439 greg 2.4
440 greg 1.8 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
441     sum = 0.0;
442     for (gs = all_srcs; gs != NULL; gs = gs->next) {
443     p = posindex(gs->dir, mydir, midview.vup);
444     if (p <= FTINY)
445     continue;
446     sum += gs->lum*gs->lum * gs->dom / (p*p);
447     }
448     if (sum <= FTINY)
449     return(0.0);
450     dillum = direct(gd);
451     return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
452 greg 2.4 }
453    
454    
455     double
456     ugr(gd) /* compute Unified Glare Rating */
457     struct glare_dir *gd;
458     {
459     register struct glare_src *gs;
460     FVECT mydir;
461     double p;
462     double sum;
463    
464     spinvector(mydir, midview.vdir, midview.vup, gd->ang);
465     sum = 0.0;
466     for (gs = all_srcs; gs != NULL; gs = gs->next) {
467     p = posindex(gs->dir, mydir, midview.vup);
468     if (p <= FTINY)
469     continue;
470     sum += gs->lum*gs->lum * gs->dom / (p*p);
471     }
472     if (sum <= FTINY)
473     return(0.0);
474     return(8.*log10(0.25*sum*PI/gd->indirect));
475 greg 1.1 }