ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
(Generate patch)

Comparing ray/src/util/glarendx.c (file contents):
Revision 1.1 by greg, Tue Apr 16 16:40:22 1991 UTC vs.
Revision 2.6 by gwlarson, Tue Oct 27 08:47:17 1998 UTC

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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines