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.5 by greg, Sun Jun 4 10:58:38 1995 UTC vs.
Revision 2.16 by greg, Fri Jul 24 16:58:16 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 or -t option:
6   *
# Line 21 | Line 18 | static char SCCSid[] = "$SunId$ LBL";
18   *              12 April 1991   Greg Ward       EPFL
19   *              19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
20   */
21 <
21 >
22 > #include <string.h>
23 >
24   #include "standard.h"
25 + #include "paths.h"
26   #include "view.h"
27 <
28 <
29 < double  posindex();
30 <
31 < double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
32 <        indirect(), total(), dgi(), brs_gi(), ugr();
33 <
27 >
28 >
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 > 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"},
# Line 47 | Line 74 | struct named_func {
74          {"vert_ind", indirect, "Indirect Vertical Illuminance"},
75          {NULL}
76   };
77 <
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 <
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 <
77 >
78   #define newp(type)      (type *)malloc(sizeof(type))
79 <
79 >
80   char    *progname;
81   int     print_header = 1;
82 <
82 >
83   VIEW    midview = STDVIEW;
84 <
84 >
85   int     wrongformat = 0;
86 <
87 <
88 < main(argc, argv)
89 < int     argc;
90 < char    *argv[];
86 >
87 >
88 > int
89 > main(
90 >        int     argc,
91 >        char    *argv[]
92 > )
93   {
78 –        extern char     *rindex();
94          struct named_func       *funp;
80 –        char    *progtail;
95          int     i;
96                                          /* get program name */
97 <        progname = argv[0];
84 <        progtail = rindex(progname, '/');       /* final component */
85 <        if (progtail == NULL)
86 <                progtail = progname;
87 <        else
88 <                progtail++;
97 >        progname = fixargv0(argv[0]);
98                                          /* get options */
99          for (i = 1; i < argc && argv[i][0] == '-'; i++)
100                  switch (argv[i][1]) {
101                  case 't':
102 <                        progtail = argv[++i];
102 >                        progname = argv[++i];
103                          break;
104                  case 'h':
105                          print_header = 0;
# Line 107 | Line 116 | char   *argv[];
116                  }
117                                          /* find and run calculation */
118          for (funp = all_funcs; funp->name != NULL; funp++)
119 <                if (!strcmp(funp->name, progtail)) {
119 >                if (!strcmp(funp->name, progname)) {
120                          init();
121                          read_input();
122                          if (print_header) {
# Line 125 | Line 134 | userr:
134                  fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
135          exit(1);
136   }
137 <
138 <
139 < headline(s)                     /* get line from header */
140 < char    *s;
137 >
138 >
139 > static int
140 > headline(                       /* get line from header */
141 >        char    *s,
142 >        void    *p
143 > )
144   {
145 <        char    fmt[32];
146 <
145 >        char    fmt[MAXFMTLEN];
146 >
147          if (print_header)               /* copy to output */
148                  fputs(s, stdout);
149          if (isview(s))
# Line 140 | Line 152 | char   *s;
152                  formatval(fmt, s);
153                  wrongformat = strcmp(fmt, "ascii");
154          }
155 +        return(0);
156   }
157 <
158 <
159 < init()                          /* initialize calculation */
157 >
158 >
159 > static void
160 > init(void)                              /* initialize calculation */
161   {
162 +        char    *err;
163                                          /* read header */
164          getheader(stdin, headline, NULL);
165          if (wrongformat) {
# Line 152 | Line 167 | init()                         /* initialize calculation */
167                  exit(1);
168          }
169                                          /* set view */
170 <        if (setview(&midview) != NULL) {
171 <                fprintf(stderr, "%s: bad view information in input\n",
170 >        if ((err = setview(&midview)) != NULL) {
171 >                fprintf(stderr, "%s: %s\n",
172 >                                progname, err);
173 >                exit(1);
174 >        }
175 >        if (fabs(DOT(midview.vdir, midview.vup)) > 1e-4) {
176 >                fprintf(stderr,
177 >                "%s: Input view direction not perpendicular to up vector\n",
178                                  progname);
179                  exit(1);
180          }
181   }
182 <
183 <
184 < read_input()                    /* read glare sources from stdin */
182 >
183 >
184 > static void
185 > read_input(void)                        /* read glare sources from stdin */
186   {
187   #define S_SEARCH        0
188   #define S_SOURCE        1
189   #define S_DIREC         2
190          int     state = S_SEARCH;
191          char    buf[128];
192 <        register struct glare_src       *gs;
193 <        register struct glare_dir       *gd;
194 <
192 >        struct glare_src        *gs;
193 >        struct glare_dir        *gd;
194 >
195          while (fgets(buf, sizeof(buf), stdin) != NULL)
196                  switch (state) {
197                  case S_SEARCH:
# Line 185 | Line 207 | read_input()                   /* read glare sources from stdin */
207                          }
208                          if ((gs = newp(struct glare_src)) == NULL)
209                                  goto memerr;
210 <                        if (sscanf(buf, "%lf %lf %lf %lf %lf",
211 <                                        &gs->dir[0], &gs->dir[1], &gs->dir[2],
212 <                                        &gs->dom, &gs->lum) != 5)
210 >                        if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
211 >                                                &gs->dir[2]) != 3 ||
212 >                                        sscanf(sskip2(buf, 3), "%lf %lf",
213 >                                                &gs->dom, &gs->lum) != 2)
214                                  goto readerr;
215                          normalize(gs->dir);
216                          gs->next = all_srcs;
# Line 219 | Line 242 | readerr:
242   #undef S_SOURCE
243   #undef S_DIREC
244   }
245 <
246 <
247 < print_values(funp)              /* print out calculations */
248 < double  (*funp)();
245 >
246 >
247 > static void
248 > print_values(           /* print out calculations */
249 >        gdfun *func
250 > )
251   {
252 <        register struct glare_dir       *gd;
253 <
252 >        struct glare_dir        *gd;
253 >
254          for (gd = all_dirs; gd != NULL; gd = gd->next)
255 <                printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
255 >                printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
256   }
257 <
258 <
259 < double
260 < direct(gd)                      /* compute direct vertical illuminance */
261 < struct glare_dir        *gd;
257 >
258 >
259 > static double
260 > direct(                 /* compute direct vertical illuminance */
261 >        struct glare_dir        *gd
262 > )
263   {
264          FVECT   mydir;
265          double  d, dval;
266 <        register struct glare_src       *gs;
267 <
266 >        struct glare_src        *gs;
267 >
268          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
269          dval = 0.0;
270          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 248 | Line 274 | struct glare_dir       *gd;
274          }
275          return(dval);
276   }
277 <
278 <
279 < double
280 < indirect(gd)                    /* return indirect vertical illuminance */
281 < struct glare_dir        *gd;
277 >
278 >
279 > static double
280 > indirect(                       /* return indirect vertical illuminance */
281 >        struct glare_dir        *gd
282 > )
283   {
284          return(gd->indirect);
285   }
286 <
287 <
288 < double
289 < total(gd)                       /* return total vertical illuminance */
290 < struct glare_dir        *gd;
286 >
287 >
288 > static double
289 > total(                  /* return total vertical illuminance */
290 >        struct glare_dir        *gd
291 > )
292   {
293          return(direct(gd)+gd->indirect);
294   }
295 <
296 <
295 >
296 >
297   /*
298   * posindex -   compute glare position index from:
299   *
# Line 279 | Line 307 | struct glare_dir       *gd;
307   * This calculation presumes the view direction and up vectors perpendicular.
308   * We return a value less than zero for improper positions.
309   */
310 <
311 < double
312 < posindex(sd, vd, vu)                    /* compute position index */
313 < FVECT   sd, vd, vu;
310 >
311 > static double
312 > posindex(                       /* compute position index */
313 >        FVECT   sd,
314 >        FVECT   vd,
315 >        FVECT   vu
316 > )
317   {
318          double  sigma, tau;
319          double  d;
320 <
320 >
321          d = DOT(sd,vd);
322          if (d <= 0.0)
323                  return(-1.0);
324 <        if (d >= 1.0)
324 >        if (d >= 1.0-FTINY)
325                  return(1.0);
326          sigma = acos(d) * (180./PI);
327          d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
328 <        if (d >= 1.0)
298 <                tau = 0.0;
299 <        else
300 <                tau = acos(d) * (180./PI);
328 >        tau = Acos(d) * (180./PI);
329          return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
330                          + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
331                  ) );
332   }
333 <
334 <
335 < double
336 < dgi(gd)         /* compute Daylight Glare Index */
337 < struct glare_dir        *gd;
333 >
334 >
335 > static double
336 > dgi(            /* compute Daylight Glare Index */
337 >        struct glare_dir        *gd
338 > )
339   {
340 <        register struct glare_src       *gs;
340 >        struct glare_src        *gs;
341          FVECT   mydir,testdir[7],vhor;
342 <        double  r,omega,p[7],sum;
343 <        int     i,n;
344 <
342 >        double  r,posn,omega,p[7],sum;
343 >        int     i;
344 >
345          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
346 <        sum = 0.0; n = 0;
346 >        sum = 0.0;
347          for (gs = all_srcs; gs != NULL; gs = gs->next) {
348  
349                  /* compute 1/p^2 weighted solid angle of the source */
350 <                r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
350 >                r = 1. - gs->dom*(.5/PI);
351 >                r = sqrt(1. - r*r);
352                  fcross(vhor,gs->dir,midview.vup);
353 <                normalize(vhor);
353 >                if (normalize(vhor) == 0.)
354 >                        continue;
355                  VCOPY(testdir[0],gs->dir);
356                  fvsum(testdir[1],gs->dir,vhor,r);
357                  fvsum(testdir[2],gs->dir,vhor,0.5*r);
# Line 331 | Line 362 | struct glare_dir       *gd;
362                  fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
363                  fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
364                  for (i = 0; i < 7; i++) {
365 <                        normalize(testdir[i]);
366 <                        p[i] = pow(posindex(testdir[i],mydir,midview.vup),-2.0);
367 <                        if (p[i] <= FTINY) p[i] = 0.0;
365 >                        p[i] = 0.0;
366 >                        if (normalize(testdir[i]) == 0.)
367 >                                continue;
368 >                        posn = posindex(testdir[i],mydir,midview.vup);
369 >                        if (posn > FTINY)
370 >                                p[i] = 1./(posn*posn);
371                  }
372 <                r = 1-gs->dom/2./PI;
372 >                r = 1. - gs->dom*(.5/PI);
373                  omega = gs->dom*p[0];
374 <                omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
375 <                omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
374 >                omega += (r*PI*(1.+1./(r*r))-2.*PI)*(-p[0]+(p[1]+p[2])*0.5);
375 >                omega += (2.*PI-r*PI*(1.+1./(r*r)))*(-p[0]-0.1667*(p[1]+p[3])
376                            +0.3334*(p[2]+p[4]+p[5]+p[6]));
377 <
378 <                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
379 <                       (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
380 <                n++;
377 >                if (omega <= 0.)
378 >                        continue;
379 >                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
380 >                       (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
381          }
382 <        if (n == 0)
382 >        if (sum <= FTINY)
383                  return(0.0);
384 <        return( 10*log10(0.478*sum) );
384 >        return( 10.*log10(0.478*sum) );
385   }
352 –
386  
387 < double
388 < brs_gi(gd)              /* compute BRS Glare Index */
389 < struct glare_dir        *gd;
387 >
388 > static double
389 > brs_gi(         /* compute BRS Glare Index */
390 >        struct glare_dir        *gd
391 > )
392   {
393 <        register struct glare_src       *gs;
393 >        struct glare_src        *gs;
394          FVECT   mydir;
395          double  p;
396          double  sum;
397 <
397 >
398          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
399          sum = 0.0;
400          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 370 | Line 405 | struct glare_dir       *gd;
405          }
406          if (sum <= FTINY)
407                  return(0.0);
408 <        sum /= gd->indirect/PI;
408 >        sum *= PI/gd->indirect;
409          return(10*log10(0.478*sum));
410   }
411  
412  
413 < double
414 < guth_dgr(gd)            /* compute Guth discomfort glare rating */
415 < struct glare_dir        *gd;
413 > static double
414 > guth_dgr(               /* compute Guth discomfort glare rating */
415 >        struct glare_dir        *gd
416 > )
417   {
418   #define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
419 <        register struct glare_src       *gs;
419 >        struct glare_src        *gs;
420          FVECT   mydir;
421          double  p;
422          double  sum;
423          double  wtot, brsum;
424          int     n;
425 <
425 >
426          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
427          sum = wtot = brsum = 0.0; n = 0;
428          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 404 | Line 440 | struct glare_dir       *gd;
440                          pow((double)n, -.0914) ) );
441   #undef q
442   }
443 <
444 <
443 >
444 >
445   #ifndef M_SQRT2
446   #define M_SQRT2 1.41421356237309504880
447   #endif
448 <
448 >
449   #define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))
450 <
451 <
452 < double
453 < guth_vcp(gd)            /* compute Guth visual comfort probability */
454 < struct glare_dir        *gd;
450 >
451 >
452 > static double
453 > guth_vcp(               /* compute Guth visual comfort probability */
454 >        struct glare_dir        *gd
455 > )
456   {
457          extern double   erfc();
458          double  dgr;
459 <
459 >
460          dgr = guth_dgr(gd);
461          if (dgr <= FTINY)
462                  return(100.0);
463          return(100.*norm_integral(6.374-1.3227*log(dgr)));
464   }
465 <
466 <
467 < double
468 < cie_cgi(gd)             /* compute CIE Glare Index */
469 < struct glare_dir        *gd;
465 >
466 >
467 > static double
468 > cie_cgi(                /* compute CIE Glare Index */
469 >        struct glare_dir        *gd
470 > )
471   {
472 <        register struct glare_src       *gs;
472 >        struct glare_src        *gs;
473          FVECT   mydir;
474          double  dillum;
475          double  p;
476          double  sum;
477 <
477 >
478          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
479          sum = 0.0;
480          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 452 | Line 490 | struct glare_dir       *gd;
490   }
491  
492  
493 < double
494 < ugr(gd)         /* compute Unified Glare Rating */
495 < struct glare_dir        *gd;
493 > static double
494 > ugr(            /* compute Unified Glare Rating */
495 >        struct glare_dir        *gd
496 > )
497   {
498 <        register struct glare_src       *gs;
498 >        struct glare_src        *gs;
499          FVECT   mydir;
500          double  p;
501          double  sum;
502 <
502 >
503          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
504          sum = 0.0;
505          for (gs = all_srcs; gs != NULL; gs = gs->next) {

Diff Legend

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