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.17 by greg, Tue Jun 3 21:31:51 2025 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  
22 + #include <string.h>
23 +
24   #include "standard.h"
25 + #include "paths.h"
26   #include "view.h"
27  
23 – extern double   erfc();
28  
29 < double  posindex();
29 > struct glare_src {
30 >        FVECT   dir;            /* source direction */
31 >        double  dom;            /* solid angle */
32 >        double  lum;            /* average luminance */
33 >        struct glare_src        *next;
34 > } *all_srcs = NULL;
35  
36 < double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
37 <        indirect(), total();
36 > struct glare_dir {
37 >        double  ang;            /* angle (in radians) */
38 >        double  indirect;       /* indirect illuminance */
39 >        struct glare_dir        *next;
40 > } *all_dirs = NULL;
41  
42 + typedef double gdfun(struct glare_dir *gd);
43 +
44 + static gdfun dgi;
45 + static gdfun brs_gi;
46 + static gdfun ugr;
47 + static gdfun guth_dgr;
48 + static gdfun guth_vcp;
49 + static gdfun cie_cgi;
50 + static gdfun direct;
51 + static gdfun indirect;
52 + static gdfun total;
53 +
54 + static gethfunc headline;
55 + static void init(void);
56 + static void read_input(void);
57 + static void print_values(gdfun *func);
58 + static double posindex(FVECT sd, FVECT vd, FVECT vu);
59 +
60 +
61   struct named_func {
62          char    *name;
63 <        double  (*func)();
63 >        gdfun   *func;
64          char    *descrip;
65   } all_funcs[] = {
66 +        {"dgi", dgi, "Daylight Glare Index"},
67 +        {"brs_gi", brs_gi, "BRS Glare Index"},
68 +        {"ugr", ugr, "Unified Glare Rating"},
69          {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
70          {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
71          {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
# Line 41 | Line 75 | struct named_func {
75          {NULL}
76   };
77  
44 – struct glare_src {
45 –        FVECT   dir;            /* source direction */
46 –        double  dom;            /* solid angle */
47 –        double  lum;            /* average luminance */
48 –        struct glare_src        *next;
49 – } *all_srcs = NULL;
50 –
51 – struct glare_dir {
52 –        double  ang;            /* angle (in radians) */
53 –        double  indirect;       /* indirect illuminance */
54 –        struct glare_dir        *next;
55 – } *all_dirs = NULL;
56 –
78   #define newp(type)      (type *)malloc(sizeof(type))
79  
59 – char    *progname;
80   int     print_header = 1;
81  
82   VIEW    midview = STDVIEW;
# Line 64 | Line 84 | VIEW   midview = STDVIEW;
84   int     wrongformat = 0;
85  
86  
87 < main(argc, argv)
88 < int     argc;
89 < char    *argv[];
87 > int
88 > main(
89 >        int     argc,
90 >        char    *argv[]
91 > )
92   {
71 –        extern char     *rindex();
93          struct named_func       *funp;
73 –        char    *progtail;
94          int     i;
95                                          /* get program name */
96 <        progname = argv[0];
77 <        progtail = rindex(progname, '/');       /* final component */
78 <        if (progtail == NULL)
79 <                progtail = progname;
80 <        else
81 <                progtail++;
96 >        progname = fixargv0(argv[0]);
97                                          /* get options */
98          for (i = 1; i < argc && argv[i][0] == '-'; i++)
99                  switch (argv[i][1]) {
100                  case 't':
101 <                        progtail = argv[++i];
101 >                        progname = argv[++i];
102                          break;
103                  case 'h':
104                          print_header = 0;
# Line 100 | Line 115 | char   *argv[];
115                  }
116                                          /* find and run calculation */
117          for (funp = all_funcs; funp->name != NULL; funp++)
118 <                if (!strcmp(funp->name, progtail)) {
118 >                if (!strcmp(funp->name, progname)) {
119                          init();
120                          read_input();
121                          if (print_header) {
# Line 120 | Line 135 | userr:
135   }
136  
137  
138 < headline(s)                     /* get line from header */
139 < char    *s;
138 > static int
139 > headline(                       /* get line from header */
140 >        char    *s,
141 >        void    *p
142 > )
143   {
144 <        char    fmt[32];
144 >        char    fmt[MAXFMTLEN];
145  
146          if (print_header)               /* copy to output */
147                  fputs(s, stdout);
# Line 133 | Line 151 | char   *s;
151                  formatval(fmt, s);
152                  wrongformat = strcmp(fmt, "ascii");
153          }
154 +        return(0);
155   }
156  
157  
158 < init()                          /* initialize calculation */
158 > static void
159 > init(void)                              /* initialize calculation */
160   {
161 +        char    *err;
162                                          /* read header */
163          getheader(stdin, headline, NULL);
164          if (wrongformat) {
# Line 145 | Line 166 | init()                         /* initialize calculation */
166                  exit(1);
167          }
168                                          /* set view */
169 <        if (setview(&midview) != NULL) {
170 <                fprintf(stderr, "%s: bad view information in input\n",
169 >        if ((err = setview(&midview)) != NULL) {
170 >                fprintf(stderr, "%s: %s\n",
171 >                                progname, err);
172 >                exit(1);
173 >        }
174 >        if (fabs(DOT(midview.vdir, midview.vup)) > 1e-4) {
175 >                fprintf(stderr,
176 >                "%s: Input view direction not perpendicular to up vector\n",
177                                  progname);
178                  exit(1);
179          }
180   }
181  
182  
183 < read_input()                    /* read glare sources from stdin */
183 > static void
184 > read_input(void)                        /* read glare sources from stdin */
185   {
186   #define S_SEARCH        0
187   #define S_SOURCE        1
188   #define S_DIREC         2
189          int     state = S_SEARCH;
190          char    buf[128];
191 <        register struct glare_src       *gs;
192 <        register struct glare_dir       *gd;
191 >        struct glare_src        *gs;
192 >        struct glare_dir        *gd;
193  
194          while (fgets(buf, sizeof(buf), stdin) != NULL)
195                  switch (state) {
# Line 178 | Line 206 | read_input()                   /* read glare sources from stdin */
206                          }
207                          if ((gs = newp(struct glare_src)) == NULL)
208                                  goto memerr;
209 <                        if (sscanf(buf, "%lf %lf %lf %lf %lf",
210 <                                        &gs->dir[0], &gs->dir[1], &gs->dir[2],
211 <                                        &gs->dom, &gs->lum) != 5)
209 >                        if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
210 >                                                &gs->dir[2]) != 3 ||
211 >                                        sscanf(sskip2(buf, 3), "%lf %lf",
212 >                                                &gs->dom, &gs->lum) != 2)
213                                  goto readerr;
214                          normalize(gs->dir);
215                          gs->next = all_srcs;
# Line 214 | Line 243 | readerr:
243   }
244  
245  
246 < print_values(funp)              /* print out calculations */
247 < double  (*funp)();
246 > static void
247 > print_values(           /* print out calculations */
248 >        gdfun *func
249 > )
250   {
251 <        register struct glare_dir       *gd;
251 >        struct glare_dir        *gd;
252  
253          for (gd = all_dirs; gd != NULL; gd = gd->next)
254 <                printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
254 >                printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
255   }
256  
257  
258 < double
259 < direct(gd)                      /* compute direct vertical illuminance */
260 < struct glare_dir        *gd;
258 > static double
259 > direct(                 /* compute direct vertical illuminance */
260 >        struct glare_dir        *gd
261 > )
262   {
263          FVECT   mydir;
264          double  d, dval;
265 <        register struct glare_src       *gs;
265 >        struct glare_src        *gs;
266  
267          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
268          dval = 0.0;
# Line 243 | Line 275 | struct glare_dir       *gd;
275   }
276  
277  
278 < double
279 < indirect(gd)                    /* return indirect vertical illuminance */
280 < struct glare_dir        *gd;
278 > static double
279 > indirect(                       /* return indirect vertical illuminance */
280 >        struct glare_dir        *gd
281 > )
282   {
283          return(gd->indirect);
284   }
285  
286  
287 < double
288 < total(gd)                       /* return total vertical illuminance */
289 < struct glare_dir        *gd;
287 > static double
288 > total(                  /* return total vertical illuminance */
289 >        struct glare_dir        *gd
290 > )
291   {
292          return(direct(gd)+gd->indirect);
293   }
# Line 273 | Line 307 | struct glare_dir       *gd;
307   * We return a value less than zero for improper positions.
308   */
309  
310 < double
311 < posindex(sd, vd, vu)                    /* compute position index */
312 < FVECT   sd, vd, vu;
310 > static double
311 > posindex(                       /* compute position index */
312 >        FVECT   sd,
313 >        FVECT   vd,
314 >        FVECT   vu
315 > )
316   {
317          double  sigma, tau;
318          double  d;
# Line 283 | Line 320 | FVECT  sd, vd, vu;
320          d = DOT(sd,vd);
321          if (d <= 0.0)
322                  return(-1.0);
323 <        if (d >= 1.0)
323 >        if (d >= 1.0-FTINY)
324                  return(1.0);
325          sigma = acos(d) * (180./PI);
326 <        d = DOT(sd,vu)/sqrt(1.0-d*d);
327 <        if (d >= 1.0)
291 <                tau = 0.0;
292 <        else if (d <= -1.0)
293 <                tau = 180.0;
294 <        else
295 <                tau = acos(d) * (180./PI);
326 >        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
327 >        tau = Acos(d) * (180./PI);
328          return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
329                          + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
330                  ) );
331   }
332  
333  
334 < double
335 < guth_dgr(gd)            /* compute Guth discomfort glare rating */
336 < struct glare_dir        *gd;
334 > static double
335 > dgi(            /* compute Daylight Glare Index */
336 >        struct glare_dir        *gd
337 > )
338   {
339 +        struct glare_src        *gs;
340 +        FVECT   mydir,testdir[7],vhor;
341 +        double  r,posn,omega,p[7],sum;
342 +        int     i;
343 +
344 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
345 +        sum = 0.0;
346 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
347 +
348 +                /* compute 1/p^2 weighted solid angle of the source */
349 +                r = 1. - gs->dom*(.5/PI);
350 +                r = sqrt(1. - r*r);
351 +                fcross(vhor,gs->dir,midview.vup);
352 +                if (normalize(vhor) == 0.)
353 +                        continue;
354 +                VCOPY(testdir[0],gs->dir);
355 +                fvsum(testdir[1],gs->dir,vhor,r);
356 +                fvsum(testdir[2],gs->dir,vhor,0.5*r);
357 +                fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
358 +                fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
359 +                fvsum(testdir[3],gs->dir,vhor,-r);
360 +                fvsum(testdir[4],gs->dir,vhor,-0.5*r);
361 +                fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
362 +                fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
363 +                for (i = 0; i < 7; i++) {
364 +                        p[i] = 0.0;
365 +                        if (normalize(testdir[i]) == 0.)
366 +                                continue;
367 +                        posn = posindex(testdir[i],mydir,midview.vup);
368 +                        if (posn > FTINY)
369 +                                p[i] = 1./(posn*posn);
370 +                }
371 +                r = 1. - gs->dom*(.5/PI);
372 +                omega = gs->dom*p[0];
373 +                omega += (r*PI*(1.+1./(r*r))-2.*PI)*(-p[0]+(p[1]+p[2])*0.5);
374 +                omega += (2.*PI-r*PI*(1.+1./(r*r)))*(-p[0]-0.1667*(p[1]+p[3])
375 +                          +0.3334*(p[2]+p[4]+p[5]+p[6]));
376 +                if (omega <= 0.)
377 +                        continue;
378 +                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
379 +                       (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
380 +        }
381 +        if (sum <= FTINY)
382 +                return(0.0);
383 +        return( 10.*log10(0.478*sum) );
384 + }
385 +
386 +
387 + static double
388 + brs_gi(         /* compute BRS Glare Index */
389 +        struct glare_dir        *gd
390 + )
391 + {
392 +        struct glare_src        *gs;
393 +        FVECT   mydir;
394 +        double  p;
395 +        double  sum;
396 +
397 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
398 +        sum = 0.0;
399 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
400 +                p = posindex(gs->dir, mydir, midview.vup);
401 +                if (p <= FTINY)
402 +                        continue;
403 +                sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
404 +        }
405 +        if (sum <= FTINY)
406 +                return(0.0);
407 +        sum *= PI/gd->indirect;
408 +        return(10*log10(0.478*sum));
409 + }
410 +
411 +
412 + static double
413 + guth_dgr(               /* compute Guth discomfort glare rating */
414 +        struct glare_dir        *gd
415 + )
416 + {
417   #define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
418 <        register struct glare_src       *gs;
418 >        struct glare_src        *gs;
419          FVECT   mydir;
420          double  p;
421          double  sum;
# Line 324 | Line 435 | struct glare_dir       *gd;
435          }
436          if (n == 0)
437                  return(0.0);
327 –
438          return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
439                          pow((double)n, -.0914) ) );
440   #undef q
441   }
442  
443  
334 – extern double   erf(), erfc();
335 –
444   #ifndef M_SQRT2
445   #define M_SQRT2 1.41421356237309504880
446   #endif
# Line 340 | Line 448 | extern double  erf(), erfc();
448   #define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))
449  
450  
451 < double
452 < guth_vcp(gd)            /* compute Guth visual comfort probability */
453 < struct glare_dir        *gd;
451 > static double
452 > guth_vcp(               /* compute Guth visual comfort probability */
453 >        struct glare_dir        *gd
454 > )
455   {
456 +        extern double   erfc();
457          double  dgr;
458  
459          dgr = guth_dgr(gd);
# Line 353 | Line 463 | struct glare_dir       *gd;
463   }
464  
465  
466 < double
467 < cie_cgi(gd)             /* compute CIE Glare Index */
468 < struct glare_dir        *gd;
466 > static double
467 > cie_cgi(                /* compute CIE Glare Index */
468 >        struct glare_dir        *gd
469 > )
470   {
471 <        register struct glare_src       *gs;
471 >        struct glare_src        *gs;
472          FVECT   mydir;
473          double  dillum;
474          double  p;
# Line 375 | Line 486 | struct glare_dir       *gd;
486                  return(0.0);
487          dillum = direct(gd);
488          return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
489 + }
490 +
491 +
492 + static double
493 + ugr(            /* compute Unified Glare Rating */
494 +        struct glare_dir        *gd
495 + )
496 + {
497 +        struct glare_src        *gs;
498 +        FVECT   mydir;
499 +        double  p;
500 +        double  sum;
501 +
502 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
503 +        sum = 0.0;
504 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
505 +                p = posindex(gs->dir, mydir, midview.vup);
506 +                if (p <= FTINY)
507 +                        continue;
508 +                sum += gs->lum*gs->lum * gs->dom / (p*p);
509 +        }
510 +        if (sum <= FTINY)
511 +                return(0.0);
512 +        return(8.*log10(0.25*sum*PI/gd->indirect));
513   }

Diff Legend

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