ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.7
Committed: Sat Feb 22 02:07:30 2003 UTC (21 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 2.6: +1 -4 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

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