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.7 by greg, Sat Feb 22 02:07:30 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1991 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * Compute Glare Index given by program name or -t option:
6   *
7 + *      dgi -           Daylight Glare Index
8 + *      brs_gi -        Building Research Station Glare Index (Petherbridge
9 + *                                                             & Hopkinson)
10 + *      ugr -           Unified Glare Rating System (Fischer)
11   *      guth_dgr -      Guth discomfort glare rating
12   *      guth_vcp -      Guth visual comfort probability
13   *      cie_cgi -       CIE Glare Index (1983, due to Einhorn)
# Line 15 | Line 16 | static char SCCSid[] = "$SunId$ LBL";
16   *      vert_ill -      Total vertical illuminance
17   *
18   *              12 April 1991   Greg Ward       EPFL
19 + *              19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
20   */
21 <
21 >
22   #include "standard.h"
23   #include "view.h"
24 <
25 < extern double   erfc();
24 <
24 >
25 >
26   double  posindex();
27 <
27 >
28   double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
29 <        indirect(), total();
30 <
29 >        indirect(), total(), dgi(), brs_gi(), ugr();
30 >
31   struct named_func {
32          char    *name;
33          double  (*func)();
34          char    *descrip;
35   } all_funcs[] = {
36 +        {"dgi", dgi, "Daylight Glare Index"},
37 +        {"brs_gi", brs_gi, "BRS Glare Index"},
38 +        {"ugr", ugr, "Unified Glare Rating"},
39          {"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"},
# Line 40 | Line 44 | struct named_func {
44          {"vert_ind", indirect, "Indirect Vertical Illuminance"},
45          {NULL}
46   };
47 <
47 >
48   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 <
54 >
55   struct glare_dir {
56          double  ang;            /* angle (in radians) */
57          double  indirect;       /* indirect illuminance */
58          struct glare_dir        *next;
59   } *all_dirs = NULL;
60 <
60 >
61   #define newp(type)      (type *)malloc(sizeof(type))
62 <
62 >
63   char    *progname;
64   int     print_header = 1;
65 <
65 >
66   VIEW    midview = STDVIEW;
67 <
67 >
68   int     wrongformat = 0;
69 <
70 <
69 >
70 >
71   main(argc, argv)
72   int     argc;
73   char    *argv[];
# Line 118 | Line 122 | userr:
122                  fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
123          exit(1);
124   }
125 <
126 <
125 >
126 >
127 > int
128   headline(s)                     /* get line from header */
129   char    *s;
130   {
131          char    fmt[32];
132 <
132 >
133          if (print_header)               /* copy to output */
134                  fputs(s, stdout);
135          if (isview(s))
# Line 133 | Line 138 | char   *s;
138                  formatval(fmt, s);
139                  wrongformat = strcmp(fmt, "ascii");
140          }
141 +        return(0);
142   }
143 <
144 <
143 >
144 >
145   init()                          /* initialize calculation */
146   {
147                                          /* read header */
# Line 151 | Line 157 | init()                         /* initialize calculation */
157                  exit(1);
158          }
159   }
160 <
161 <
160 >
161 >
162   read_input()                    /* read glare sources from stdin */
163   {
164   #define S_SEARCH        0
# Line 162 | Line 168 | read_input()                   /* read glare sources from stdin */
168          char    buf[128];
169          register struct glare_src       *gs;
170          register struct glare_dir       *gd;
171 <
171 >
172          while (fgets(buf, sizeof(buf), stdin) != NULL)
173                  switch (state) {
174                  case S_SEARCH:
# Line 212 | Line 218 | readerr:
218   #undef S_SOURCE
219   #undef S_DIREC
220   }
221 <
222 <
221 >
222 >
223   print_values(funp)              /* print out calculations */
224   double  (*funp)();
225   {
226          register struct glare_dir       *gd;
227 <
227 >
228          for (gd = all_dirs; gd != NULL; gd = gd->next)
229                  printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
230   }
231 <
232 <
231 >
232 >
233   double
234   direct(gd)                      /* compute direct vertical illuminance */
235   struct glare_dir        *gd;
# Line 231 | Line 237 | struct glare_dir       *gd;
237          FVECT   mydir;
238          double  d, dval;
239          register struct glare_src       *gs;
240 <
240 >
241          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
242          dval = 0.0;
243          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 241 | Line 247 | struct glare_dir       *gd;
247          }
248          return(dval);
249   }
250 <
251 <
250 >
251 >
252   double
253   indirect(gd)                    /* return indirect vertical illuminance */
254   struct glare_dir        *gd;
255   {
256          return(gd->indirect);
257   }
258 <
259 <
258 >
259 >
260   double
261   total(gd)                       /* return total vertical illuminance */
262   struct glare_dir        *gd;
263   {
264          return(direct(gd)+gd->indirect);
265   }
266 <
267 <
266 >
267 >
268   /*
269   * posindex -   compute glare position index from:
270   *
# Line 272 | Line 278 | struct glare_dir       *gd;
278   * This calculation presumes the view direction and up vectors perpendicular.
279   * We return a value less than zero for improper positions.
280   */
281 <
281 >
282   double
283   posindex(sd, vd, vu)                    /* compute position index */
284   FVECT   sd, vd, vu;
285   {
286          double  sigma, tau;
287          double  d;
288 <
288 >
289          d = DOT(sd,vd);
290          if (d <= 0.0)
291                  return(-1.0);
292          if (d >= 1.0)
293                  return(1.0);
294          sigma = acos(d) * (180./PI);
295 <        d = DOT(sd,vu)/sqrt(1.0-d*d);
295 >        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
296          if (d >= 1.0)
297                  tau = 0.0;
292 –        else if (d <= -1.0)
293 –                tau = 180.0;
298          else
299                  tau = acos(d) * (180./PI);
300          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 +
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  
319 +                /* 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  
343 +                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   double
354 + 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   guth_dgr(gd)            /* compute Guth discomfort glare rating */
379   struct glare_dir        *gd;
380   {
# Line 310 | Line 385 | struct glare_dir       *gd;
385          double  sum;
386          double  wtot, brsum;
387          int     n;
388 <
388 >
389          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
390          sum = wtot = brsum = 0.0; n = 0;
391          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 324 | Line 399 | struct glare_dir       *gd;
399          }
400          if (n == 0)
401                  return(0.0);
327 –
402          return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
403                          pow((double)n, -.0914) ) );
404   #undef q
405   }
406 <
407 <
334 < extern double   erf(), erfc();
335 <
406 >
407 >
408   #ifndef M_SQRT2
409   #define M_SQRT2 1.41421356237309504880
410   #endif
411 <
411 >
412   #define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))
413 <
414 <
413 >
414 >
415   double
416   guth_vcp(gd)            /* compute Guth visual comfort probability */
417   struct glare_dir        *gd;
418   {
419 +        extern double   erfc();
420          double  dgr;
421 <
421 >
422          dgr = guth_dgr(gd);
423          if (dgr <= FTINY)
424                  return(100.0);
425          return(100.*norm_integral(6.374-1.3227*log(dgr)));
426   }
427 <
428 <
427 >
428 >
429   double
430   cie_cgi(gd)             /* compute CIE Glare Index */
431   struct glare_dir        *gd;
# Line 362 | Line 435 | struct glare_dir       *gd;
435          double  dillum;
436          double  p;
437          double  sum;
438 <
438 >
439          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
440          sum = 0.0;
441          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 375 | Line 448 | struct glare_dir       *gd;
448                  return(0.0);
449          dillum = direct(gd);
450          return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
451 + }
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   }

Diff Legend

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