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.15 by greg, Mon Feb 17 19:35:24 2020 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:
5 > * Compute Glare Index given by program name or -t option:
6   *
7 < *      gnuth_dgr -     Gnuth discomfort glare rating
8 < *      gnuth_vcp -     Gnuth visual comfort probability
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)
14 > *      vert_dir -      Direct vertical illuminance
15 > *      vert_ind -      Indirect vertical illuminance (from input)
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  
22 + #include <string.h>
23 +
24   #include "standard.h"
25   #include "view.h"
26  
19 – extern double   erfc();
27  
21 – double  posindex();
22 – int     headline();
23 –
24 – double  direct(), gnuth_dgr(), gnuth_vcp();
25 –
26 – struct named_func {
27 –        char    *name;
28 –        double  (*func)();
29 – } all_funcs[] = {
30 –        {"direct", direct},
31 –        {"gnuth_dgr", gnuth_dgr},
32 –        {"gnuth_vcp", gnuth_vcp},
33 –        {NULL}
34 – };
35 –
28   struct glare_src {
29          FVECT   dir;            /* source direction */
30          double  dom;            /* solid angle */
# Line 46 | Line 38 | struct glare_dir {
38          struct glare_dir        *next;
39   } *all_dirs = NULL;
40  
41 + typedef double gdfun(struct glare_dir *gd);
42 +
43 + static gdfun dgi;
44 + static gdfun brs_gi;
45 + static gdfun ugr;
46 + static gdfun guth_dgr;
47 + static gdfun guth_vcp;
48 + static gdfun cie_cgi;
49 + static gdfun direct;
50 + static gdfun indirect;
51 + static gdfun total;
52 +
53 + static gethfunc headline;
54 + static void init(void);
55 + static void read_input(void);
56 + static void print_values(gdfun *func);
57 + static double posindex(FVECT sd, FVECT vd, FVECT vu);
58 +
59 +
60 + struct named_func {
61 +        char    *name;
62 +        gdfun   *func;
63 +        char    *descrip;
64 + } all_funcs[] = {
65 +        {"dgi", dgi, "Daylight Glare Index"},
66 +        {"brs_gi", brs_gi, "BRS Glare Index"},
67 +        {"ugr", ugr, "Unified Glare Rating"},
68 +        {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
69 +        {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
70 +        {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
71 +        {"vert_dir", direct, "Direct Vertical Illuminance"},
72 +        {"vert_ill", total, "Total Vertical Illuminance"},
73 +        {"vert_ind", indirect, "Indirect Vertical Illuminance"},
74 +        {NULL}
75 + };
76 +
77   #define newp(type)      (type *)malloc(sizeof(type))
78  
79   char    *progname;
# Line 53 | Line 81 | int    print_header = 1;
81  
82   VIEW    midview = STDVIEW;
83  
84 + int     wrongformat = 0;
85  
86 < main(argc, argv)
87 < int     argc;
88 < char    *argv[];
86 >
87 > int
88 > main(
89 >        int     argc,
90 >        char    *argv[]
91 > )
92   {
61 –        extern char     *rindex();
93          struct named_func       *funp;
94          char    *progtail;
95          int     i;
96                                          /* get program name */
97          progname = argv[0];
98 <        progtail = rindex(progname, '/');       /* final component */
98 >        progtail = strrchr(progname, '/');      /* final component */
99          if (progtail == NULL)
100                  progtail = progname;
101          else
# Line 88 | Line 119 | char   *argv[];
119                          perror(argv[i]);
120                          exit(1);
121                  }
122 <                                        /* 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 */
122 >                                        /* find and run calculation */
123          for (funp = all_funcs; funp->name != NULL; funp++)
124                  if (!strcmp(funp->name, progtail)) {
125 +                        init();
126 +                        read_input();
127 +                        if (print_header) {
128 +                                printargs(i, argv, stdout);
129 +                                putchar('\n');
130 +                        }
131                          print_values(funp->func);
132                          exit(0);                /* we're done */
133                  }
134                                          /* invalid function */
111 –        fprintf(stderr, "%s: unknown function!\n", progtail);
112 –        exit(1);
135   userr:
136 <        fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
136 >        fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
137 >        fprintf(stderr, "\twhere 'type' is one of the following:\n");
138 >        for (funp = all_funcs; funp->name != NULL; funp++)
139 >                fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
140          exit(1);
141   }
142  
143  
144 < headline(s)                     /* get line from header */
145 < char    *s;
144 > static int
145 > headline(                       /* get line from header */
146 >        char    *s,
147 >        void    *p
148 > )
149   {
150 +        char    fmt[MAXFMTLEN];
151 +
152          if (print_header)               /* copy to output */
153                  fputs(s, stdout);
154 <        if (!strncmp(s, VIEWSTR, VIEWSTRL))
155 <                sscanview(&midview, s+VIEWSTRL);
154 >        if (isview(s))
155 >                sscanview(&midview, s);
156 >        else if (isformat(s)) {
157 >                formatval(fmt, s);
158 >                wrongformat = strcmp(fmt, "ascii");
159 >        }
160 >        return(0);
161   }
162  
163  
164 < read_input()                    /* read glare sources from stdin */
164 > static void
165 > init(void)                              /* initialize calculation */
166   {
167 +        char    *err;
168 +                                        /* read header */
169 +        getheader(stdin, headline, NULL);
170 +        if (wrongformat) {
171 +                fprintf(stderr, "%s: bad input format\n", progname);
172 +                exit(1);
173 +        }
174 +                                        /* set view */
175 +        if ((err = setview(&midview)) != NULL) {
176 +                fprintf(stderr, "%s: %s\n",
177 +                                progname, err);
178 +                exit(1);
179 +        }
180 +        if (fabs(DOT(midview.vdir, midview.vup)) > 1e-4) {
181 +                fprintf(stderr,
182 +                "%s: Input view direction not perpendicular to up vector\n",
183 +                                progname);
184 +                exit(1);
185 +        }
186 + }
187 +
188 +
189 + static void
190 + read_input(void)                        /* read glare sources from stdin */
191 + {
192   #define S_SEARCH        0
193   #define S_SOURCE        1
194   #define S_DIREC         2
195          int     state = S_SEARCH;
196          char    buf[128];
197 <        register struct glare_src       *gs;
198 <        register struct glare_dir       *gd;
197 >        struct glare_src        *gs;
198 >        struct glare_dir        *gd;
199  
200          while (fgets(buf, sizeof(buf), stdin) != NULL)
201                  switch (state) {
# Line 151 | Line 212 | read_input()                   /* read glare sources from stdin */
212                          }
213                          if ((gs = newp(struct glare_src)) == NULL)
214                                  goto memerr;
215 <                        if (sscanf(buf, "%lf %lf %lf %lf %lf",
216 <                                        &gs->dir[0], &gs->dir[1], &gs->dir[2],
217 <                                        &gs->dom, &gs->lum) != 5)
215 >                        if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
216 >                                                &gs->dir[2]) != 3 ||
217 >                                        sscanf(sskip2(buf, 3), "%lf %lf",
218 >                                                &gs->dom, &gs->lum) != 2)
219                                  goto readerr;
220 +                        normalize(gs->dir);
221                          gs->next = all_srcs;
222                          all_srcs = gs;
223                          break;
# Line 186 | Line 249 | readerr:
249   }
250  
251  
252 < print_values(funp)              /* print out calculations */
253 < double  (*funp)();
252 > static void
253 > print_values(           /* print out calculations */
254 >        gdfun *func
255 > )
256   {
257 <        register struct glare_dir       *gd;
257 >        struct glare_dir        *gd;
258  
259          for (gd = all_dirs; gd != NULL; gd = gd->next)
260 <                printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
260 >                printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
261   }
262  
263  
264 < double
265 < direct(gd)                      /* compute direct illuminance */
266 < struct glare_dir        *gd;
264 > static double
265 > direct(                 /* compute direct vertical illuminance */
266 >        struct glare_dir        *gd
267 > )
268   {
269          FVECT   mydir;
270          double  d, dval;
271 <        register struct glare_src       *gs;
271 >        struct glare_src        *gs;
272  
273          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
274          dval = 0.0;
# Line 215 | Line 281 | struct glare_dir       *gd;
281   }
282  
283  
284 + static double
285 + indirect(                       /* return indirect vertical illuminance */
286 +        struct glare_dir        *gd
287 + )
288 + {
289 +        return(gd->indirect);
290 + }
291 +
292 +
293 + static double
294 + total(                  /* return total vertical illuminance */
295 +        struct glare_dir        *gd
296 + )
297 + {
298 +        return(direct(gd)+gd->indirect);
299 + }
300 +
301 +
302   /*
303   * posindex -   compute glare position index from:
304   *
# Line 225 | Line 309 | struct glare_dir       *gd;
309   * All vectors are assumed to be normalized.
310   * This function is an implementation of the method proposed by
311   * Robert Levin in his 1975 JIES article.
312 + * This calculation presumes the view direction and up vectors perpendicular.
313   * We return a value less than zero for improper positions.
314   */
315  
316 < double
317 < posindex(sd, vd, vu)                    /* compute position index */
318 < FVECT   sd, vd, vu;
316 > static double
317 > posindex(                       /* compute position index */
318 >        FVECT   sd,
319 >        FVECT   vd,
320 >        FVECT   vu
321 > )
322   {
323          double  sigma, tau;
324          double  d;
# Line 238 | Line 326 | FVECT  sd, vd, vu;
326          d = DOT(sd,vd);
327          if (d <= 0.0)
328                  return(-1.0);
329 +        if (d >= 1.0-FTINY)
330 +                return(1.0);
331          sigma = acos(d) * (180./PI);
332 <        tau = acos(DOT(sd,vu)/sqrt(1.0-d*d)) * (180./PI);
332 >        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
333 >        tau = Acos(d) * (180./PI);
334          return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
335                          + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
336                  ) );
337   }
338  
339  
340 < double
341 < gnuth_dgr(gd)           /* compute Gnuth discomfort glare rating */
342 < struct glare_dir        *gd;
340 > static double
341 > dgi(            /* compute Daylight Glare Index */
342 >        struct glare_dir        *gd
343 > )
344   {
345 +        struct glare_src        *gs;
346 +        FVECT   mydir,testdir[7],vhor;
347 +        double  r,posn,omega,p[7],sum;
348 +        int     i;
349 +
350 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
351 +        sum = 0.0;
352 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
353 +
354 +                /* compute 1/p^2 weighted solid angle of the source */
355 +                r = 1. - gs->dom*(.5/PI);
356 +                r = sqrt(1. - r*r);
357 +                fcross(vhor,gs->dir,midview.vup);
358 +                if (normalize(vhor) == 0.)
359 +                        continue;
360 +                VCOPY(testdir[0],gs->dir);
361 +                fvsum(testdir[1],gs->dir,vhor,r);
362 +                fvsum(testdir[2],gs->dir,vhor,0.5*r);
363 +                fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
364 +                fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
365 +                fvsum(testdir[3],gs->dir,vhor,-r);
366 +                fvsum(testdir[4],gs->dir,vhor,-0.5*r);
367 +                fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
368 +                fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
369 +                for (i = 0; i < 7; i++) {
370 +                        p[i] = 0.0;
371 +                        if (normalize(testdir[i]) == 0.)
372 +                                continue;
373 +                        posn = posindex(testdir[i],mydir,midview.vup);
374 +                        if (posn > FTINY)
375 +                                p[i] = 1./(posn*posn);
376 +                }
377 +                r = 1. - gs->dom*(.5/PI);
378 +                omega = gs->dom*p[0];
379 +                omega += (r*PI*(1.+1./(r*r))-2.*PI)*(-p[0]+(p[1]+p[2])*0.5);
380 +                omega += (2.*PI-r*PI*(1.+1./(r*r)))*(-p[0]-0.1667*(p[1]+p[3])
381 +                          +0.3334*(p[2]+p[4]+p[5]+p[6]));
382 +                if (omega <= 0.)
383 +                        continue;
384 +                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
385 +                       (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
386 +        }
387 +        if (sum <= FTINY)
388 +                return(0.0);
389 +        return( 10.*log10(0.478*sum) );
390 + }
391 +
392 +
393 + static double
394 + brs_gi(         /* compute BRS Glare Index */
395 +        struct glare_dir        *gd
396 + )
397 + {
398 +        struct glare_src        *gs;
399 +        FVECT   mydir;
400 +        double  p;
401 +        double  sum;
402 +
403 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
404 +        sum = 0.0;
405 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
406 +                p = posindex(gs->dir, mydir, midview.vup);
407 +                if (p <= FTINY)
408 +                        continue;
409 +                sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
410 +        }
411 +        if (sum <= FTINY)
412 +                return(0.0);
413 +        sum *= PI/gd->indirect;
414 +        return(10*log10(0.478*sum));
415 + }
416 +
417 +
418 + static double
419 + guth_dgr(               /* compute Guth discomfort glare rating */
420 +        struct glare_dir        *gd
421 + )
422 + {
423   #define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
424 <        register struct glare_src       *gs;
424 >        struct glare_src        *gs;
425 >        FVECT   mydir;
426          double  p;
427          double  sum;
428 +        double  wtot, brsum;
429          int     n;
430  
431 <        sum = 0.0; n = 0;
431 >        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
432 >        sum = wtot = brsum = 0.0; n = 0;
433          for (gs = all_srcs; gs != NULL; gs = gs->next) {
434 <                p = posindex(gs->dir, midview.vdir, midview.vup);
434 >                p = posindex(gs->dir, mydir, midview.vup);
435                  if (p <= FTINY)
436                          continue;
437                  sum += gs->lum * q(gs->dom) / p;
438 +                brsum += gs->lum * gs->dom;
439 +                wtot += gs->dom;
440                  n++;
441          }
442          if (n == 0)
443                  return(0.0);
444 <        else
270 <                return( pow(
271 <                        .5*sum/pow(direct(gd)+gd->indirect,.44),
444 >        return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
445                          pow((double)n, -.0914) ) );
446   #undef q
447   }
448  
449  
277 – extern double   erf(), erfc();
278 –
450   #ifndef M_SQRT2
451   #define M_SQRT2 1.41421356237309504880
452   #endif
# Line 283 | Line 454 | extern double  erf(), erfc();
454   #define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))
455  
456  
457 < double
458 < gnuth_vcp(gd)           /* compute Gnuth visual comfort probability */
459 < struct glare_dir        *gd;
457 > static double
458 > guth_vcp(               /* compute Guth visual comfort probability */
459 >        struct glare_dir        *gd
460 > )
461   {
462 <        return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd))));
462 >        extern double   erfc();
463 >        double  dgr;
464 >
465 >        dgr = guth_dgr(gd);
466 >        if (dgr <= FTINY)
467 >                return(100.0);
468 >        return(100.*norm_integral(6.374-1.3227*log(dgr)));
469 > }
470 >
471 >
472 > static double
473 > cie_cgi(                /* compute CIE Glare Index */
474 >        struct glare_dir        *gd
475 > )
476 > {
477 >        struct glare_src        *gs;
478 >        FVECT   mydir;
479 >        double  dillum;
480 >        double  p;
481 >        double  sum;
482 >
483 >        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
484 >        sum = 0.0;
485 >        for (gs = all_srcs; gs != NULL; gs = gs->next) {
486 >                p = posindex(gs->dir, mydir, midview.vup);
487 >                if (p <= FTINY)
488 >                        continue;
489 >                sum += gs->lum*gs->lum * gs->dom / (p*p);
490 >        }
491 >        if (sum <= FTINY)
492 >                return(0.0);
493 >        dillum = direct(gd);
494 >        return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
495 > }
496 >
497 >
498 > static double
499 > ugr(            /* compute Unified Glare Rating */
500 >        struct glare_dir        *gd
501 > )
502 > {
503 >        struct glare_src        *gs;
504 >        FVECT   mydir;
505 >        double  p;
506 >        double  sum;
507 >
508 >        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
509 >        sum = 0.0;
510 >        for (gs = all_srcs; gs != NULL; gs = gs->next) {
511 >                p = posindex(gs->dir, mydir, midview.vup);
512 >                if (p <= FTINY)
513 >                        continue;
514 >                sum += gs->lum*gs->lum * gs->dom / (p*p);
515 >        }
516 >        if (sum <= FTINY)
517 >                return(0.0);
518 >        return(8.*log10(0.25*sum*PI/gd->indirect));
519   }

Diff Legend

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