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 2.3 by greg, Thu Nov 19 09:59:08 1992 UTC vs.
Revision 2.4 by greg, Thu Mar 24 11:06:59 1994 UTC

# Line 7 | Line 7 | static char SCCSid[] = "$SunId$ LBL";
7   /*
8   * Compute Glare Index given by program name or -t option:
9   *
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)
# Line 15 | Line 19 | static char SCCSid[] = "$SunId$ LBL";
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 <
27 >
28   extern double   erfc();
29 <
29 >
30   double  posindex();
31 <
31 >
32   double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
33 <        indirect(), total();
34 <
33 >        indirect(), total(), dgi(), brs_gi(), ugr();
34 >
35   struct named_func {
36          char    *name;
37          double  (*func)();
38          char    *descrip;
39   } all_funcs[] = {
40 +        {"dgi", dgi, "Daylight Glare Index"},
41 +        {"brs_gi", brs_gi, "BRS Glare Index"},
42 +        {"ugr", ugr, "Unified Glare Rating"},
43          {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
44          {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
45          {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
# Line 40 | Line 48 | struct named_func {
48          {"vert_ind", indirect, "Indirect Vertical Illuminance"},
49          {NULL}
50   };
51 <
51 >
52   struct glare_src {
53          FVECT   dir;            /* source direction */
54          double  dom;            /* solid angle */
55          double  lum;            /* average luminance */
56          struct glare_src        *next;
57   } *all_srcs = NULL;
58 <
58 >
59   struct glare_dir {
60          double  ang;            /* angle (in radians) */
61          double  indirect;       /* indirect illuminance */
62          struct glare_dir        *next;
63   } *all_dirs = NULL;
64 <
64 >
65   #define newp(type)      (type *)malloc(sizeof(type))
66 <
66 >
67   char    *progname;
68   int     print_header = 1;
69 <
69 >
70   VIEW    midview = STDVIEW;
71 <
71 >
72   int     wrongformat = 0;
73 <
74 <
73 >
74 >
75   main(argc, argv)
76   int     argc;
77   char    *argv[];
# Line 118 | Line 126 | userr:
126                  fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
127          exit(1);
128   }
129 <
130 <
129 >
130 >
131   headline(s)                     /* get line from header */
132   char    *s;
133   {
134          char    fmt[32];
135 <
135 >
136          if (print_header)               /* copy to output */
137                  fputs(s, stdout);
138          if (isview(s))
# Line 134 | Line 142 | char   *s;
142                  wrongformat = strcmp(fmt, "ascii");
143          }
144   }
145 <
146 <
145 >
146 >
147   init()                          /* initialize calculation */
148   {
149                                          /* read header */
# Line 151 | Line 159 | init()                         /* initialize calculation */
159                  exit(1);
160          }
161   }
162 <
163 <
162 >
163 >
164   read_input()                    /* read glare sources from stdin */
165   {
166   #define S_SEARCH        0
# Line 162 | Line 170 | read_input()                   /* read glare sources from stdin */
170          char    buf[128];
171          register struct glare_src       *gs;
172          register struct glare_dir       *gd;
173 <
173 >
174          while (fgets(buf, sizeof(buf), stdin) != NULL)
175                  switch (state) {
176                  case S_SEARCH:
# Line 212 | Line 220 | readerr:
220   #undef S_SOURCE
221   #undef S_DIREC
222   }
223 <
224 <
223 >
224 >
225   print_values(funp)              /* print out calculations */
226   double  (*funp)();
227   {
228          register struct glare_dir       *gd;
229 <
229 >
230          for (gd = all_dirs; gd != NULL; gd = gd->next)
231                  printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
232   }
233 <
234 <
233 >
234 >
235   double
236   direct(gd)                      /* compute direct vertical illuminance */
237   struct glare_dir        *gd;
# Line 231 | Line 239 | struct glare_dir       *gd;
239          FVECT   mydir;
240          double  d, dval;
241          register struct glare_src       *gs;
242 <
242 >
243          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
244          dval = 0.0;
245          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 241 | Line 249 | struct glare_dir       *gd;
249          }
250          return(dval);
251   }
252 <
253 <
252 >
253 >
254   double
255   indirect(gd)                    /* return indirect vertical illuminance */
256   struct glare_dir        *gd;
257   {
258          return(gd->indirect);
259   }
260 <
261 <
260 >
261 >
262   double
263   total(gd)                       /* return total vertical illuminance */
264   struct glare_dir        *gd;
265   {
266          return(direct(gd)+gd->indirect);
267   }
268 <
269 <
268 >
269 >
270   /*
271   * posindex -   compute glare position index from:
272   *
# Line 272 | Line 280 | struct glare_dir       *gd;
280   * This calculation presumes the view direction and up vectors perpendicular.
281   * We return a value less than zero for improper positions.
282   */
283 <
283 >
284   double
285   posindex(sd, vd, vu)                    /* compute position index */
286   FVECT   sd, vd, vu;
287   {
288          double  sigma, tau;
289          double  d;
290 <
290 >
291          d = DOT(sd,vd);
292          if (d <= 0.0)
293                  return(-1.0);
294          if (d >= 1.0)
295                  return(1.0);
296          sigma = acos(d) * (180./PI);
297 <        d = DOT(sd,vu)/sqrt(1.0-d*d);
297 >        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
298          if (d >= 1.0)
299                  tau = 0.0;
292 –        else if (d <= -1.0)
293 –                tau = 180.0;
300          else
301                  tau = acos(d) * (180./PI);
302          return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
303                          + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
304                  ) );
305   }
306 +
307 +
308 + double
309 + dgi(gd)         /* compute Daylight Glare Index */
310 + struct glare_dir        *gd;
311 + {
312 +        register struct glare_src       *gs;
313 +        FVECT   mydir,testdir[7],vhor;
314 +        double  r,omega,p[7],sum;
315 +        int     i,n;
316 +
317 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
318 +        sum = 0.0; n = 0;
319 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
320  
321 +                /* compute 1/p^2 weighted solid angle of the source */
322 +                r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
323 +                fcross(vhor,gs->dir,midview.vup);
324 +                normalize(vhor);
325 +                VCOPY(testdir[0],gs->dir);
326 +                fvsum(testdir[1],gs->dir,vhor,r);
327 +                fvsum(testdir[2],gs->dir,vhor,0.5*r);
328 +                fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
329 +                fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
330 +                fvsum(testdir[3],gs->dir,vhor,-r);
331 +                fvsum(testdir[4],gs->dir,vhor,-0.5*r);
332 +                fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
333 +                fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
334 +                for (i = 0; i < 7; i++) {
335 +                        normalize(testdir[i]);
336 +                        p[i] = pow(posindex(testdir[i],mydir,midview.vup),-2.0);
337 +                        if (p[i] <= FTINY) p[i] = 0.0;
338 +                }
339 +                r = 1-gs->dom/2./PI;
340 +                omega = gs->dom*p[0];
341 +                omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
342 +                omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
343 +                          +0.3334*(p[2]+p[4]+p[5]+p[6]));
344  
345 +                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
346 +                       (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
347 +                n++;
348 +        }
349 +        if (n == 0)
350 +                return(0.0);
351 +        return( 10*log10(0.478*sum) );
352 + }
353 +
354 +
355   double
356 + brs_gi(gd)              /* compute BRS Glare Index */
357 + struct glare_dir        *gd;
358 + {
359 +        register struct glare_src       *gs;
360 +        FVECT   mydir;
361 +        double  p;
362 +        double  sum;
363 +
364 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
365 +        sum = 0.0;
366 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
367 +                p = posindex(gs->dir, mydir, midview.vup);
368 +                if (p <= FTINY)
369 +                        continue;
370 +                sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
371 +        }
372 +        if (sum <= FTINY)
373 +                return(0.0);
374 +        sum /= gd->indirect/PI;
375 +        return(10*log10(0.478*sum));
376 + }
377 +
378 +
379 + double
380   guth_dgr(gd)            /* compute Guth discomfort glare rating */
381   struct glare_dir        *gd;
382   {
# Line 310 | Line 387 | struct glare_dir       *gd;
387          double  sum;
388          double  wtot, brsum;
389          int     n;
390 <
390 >
391          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
392          sum = wtot = brsum = 0.0; n = 0;
393          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 324 | Line 401 | struct glare_dir       *gd;
401          }
402          if (n == 0)
403                  return(0.0);
327 –
404          return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
405                          pow((double)n, -.0914) ) );
406   #undef q
407   }
408 <
409 <
408 >
409 >
410   extern double   erf(), erfc();
411 <
411 >
412   #ifndef M_SQRT2
413   #define M_SQRT2 1.41421356237309504880
414   #endif
415 <
415 >
416   #define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))
417 <
418 <
417 >
418 >
419   double
420   guth_vcp(gd)            /* compute Guth visual comfort probability */
421   struct glare_dir        *gd;
422   {
423          double  dgr;
424 <
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 <
430 >
431 >
432   double
433   cie_cgi(gd)             /* compute CIE Glare Index */
434   struct glare_dir        *gd;
# Line 362 | Line 438 | struct glare_dir       *gd;
438          double  dillum;
439          double  p;
440          double  sum;
441 <
441 >
442          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
443          sum = 0.0;
444          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 375 | Line 451 | struct glare_dir       *gd;
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 + ugr(gd)         /* compute Unified Glare Rating */
459 + struct glare_dir        *gd;
460 + {
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