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 1.9 by greg, Thu May 2 14:00:58 1991 UTC

# Line 5 | Line 5 | static char SCCSid[] = "$SunId$ LBL";
5   #endif
6  
7   /*
8 < * Compute Glare Index given by program name:
8 > * Compute Glare Index given by program name or -t option:
9   *
10 < *      gnuth_dgr -     Gnuth discomfort glare rating
11 < *      gnuth_vcp -     Gnuth visual comfort probability
10 > *      guth_dgr -      Guth discomfort glare rating
11 > *      guth_vcp -      Guth visual comfort probability
12 > *      cie_cgi -       CIE Glare Index (1983, due to Einhorn)
13 > *      vert_dir -      Direct vertical illuminance
14 > *      vert_ind -      Indirect vertical illuminance (from input)
15 > *      vert_ill -      Total vertical illuminance
16   *
17   *              12 April 1991   Greg Ward       EPFL
18   */
# Line 19 | Line 23 | static char SCCSid[] = "$SunId$ LBL";
23   extern double   erfc();
24  
25   double  posindex();
22 int     headline();
26  
27 < double  direct(), gnuth_dgr(), gnuth_vcp();
27 > double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
28 >        indirect(), total();
29  
30   struct named_func {
31          char    *name;
32          double  (*func)();
33 +        char    *descrip;
34   } all_funcs[] = {
35 <        {"direct", direct},
36 <        {"gnuth_dgr", gnuth_dgr},
37 <        {"gnuth_vcp", gnuth_vcp},
35 >        {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
36 >        {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
37 >        {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
38 >        {"vert_dir", direct, "Direct Vertical Illuminance"},
39 >        {"vert_ill", total, "Total Vertical Illuminance"},
40 >        {"vert_ind", indirect, "Indirect Vertical Illuminance"},
41          {NULL}
42   };
43  
# Line 53 | Line 61 | int    print_header = 1;
61  
62   VIEW    midview = STDVIEW;
63  
64 + int     wrongformat = 0;
65  
66 +
67   main(argc, argv)
68   int     argc;
69   char    *argv[];
# Line 88 | Line 98 | char   *argv[];
98                          perror(argv[i]);
99                          exit(1);
100                  }
101 <                                        /* 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 */
101 >                                        /* find and run calculation */
102          for (funp = all_funcs; funp->name != NULL; funp++)
103                  if (!strcmp(funp->name, progtail)) {
104 +                        init();
105 +                        read_input();
106 +                        if (print_header) {
107 +                                printargs(i, argv, stdout);
108 +                                putchar('\n');
109 +                        }
110                          print_values(funp->func);
111                          exit(0);                /* we're done */
112                  }
113                                          /* invalid function */
111        fprintf(stderr, "%s: unknown function!\n", progtail);
112        exit(1);
114   userr:
115 <        fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
115 >        fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
116 >        fprintf(stderr, "\twhere 'type' is one of the following:\n");
117 >        for (funp = all_funcs; funp->name != NULL; funp++)
118 >                fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
119          exit(1);
120   }
121  
# Line 119 | Line 123 | userr:
123   headline(s)                     /* get line from header */
124   char    *s;
125   {
126 +        char    fmt[32];
127 +
128          if (print_header)               /* copy to output */
129                  fputs(s, stdout);
130          if (!strncmp(s, VIEWSTR, VIEWSTRL))
131                  sscanview(&midview, s+VIEWSTRL);
132 +        else if (isformat(s)) {
133 +                formatval(fmt, s);
134 +                wrongformat = strcmp(fmt, "ascii");
135 +        }
136   }
137  
138  
139 + init()                          /* initialize calculation */
140 + {
141 +                                        /* read header */
142 +        getheader(stdin, headline, NULL);
143 +        if (wrongformat) {
144 +                fprintf(stderr, "%s: bad input format\n", progname);
145 +                exit(1);
146 +        }
147 +                                        /* set view */
148 +        if (setview(&midview) != NULL) {
149 +                fprintf(stderr, "%s: bad view information in input\n",
150 +                                progname);
151 +                exit(1);
152 +        }
153 + }
154 +
155 +
156   read_input()                    /* read glare sources from stdin */
157   {
158   #define S_SEARCH        0
# Line 155 | Line 182 | read_input()                   /* read glare sources from stdin */
182                                          &gs->dir[0], &gs->dir[1], &gs->dir[2],
183                                          &gs->dom, &gs->lum) != 5)
184                                  goto readerr;
185 +                        normalize(gs->dir);
186                          gs->next = all_srcs;
187                          all_srcs = gs;
188                          break;
# Line 197 | Line 225 | double (*funp)();
225  
226  
227   double
228 < direct(gd)                      /* compute direct illuminance */
228 > direct(gd)                      /* compute direct vertical illuminance */
229   struct glare_dir        *gd;
230   {
231          FVECT   mydir;
# Line 215 | Line 243 | struct glare_dir       *gd;
243   }
244  
245  
246 + double
247 + indirect(gd)                    /* return indirect vertical illuminance */
248 + struct glare_dir        *gd;
249 + {
250 +        return(gd->indirect);
251 + }
252 +
253 +
254 + double
255 + total(gd)                       /* return total vertical illuminance */
256 + struct glare_dir        *gd;
257 + {
258 +        return(direct(gd)+gd->indirect);
259 + }
260 +
261 +
262   /*
263   * posindex -   compute glare position index from:
264   *
# Line 225 | Line 269 | struct glare_dir       *gd;
269   * All vectors are assumed to be normalized.
270   * This function is an implementation of the method proposed by
271   * Robert Levin in his 1975 JIES article.
272 + * This calculation presumes the view direction and up vectors perpendicular.
273   * We return a value less than zero for improper positions.
274   */
275  
# Line 238 | Line 283 | FVECT  sd, vd, vu;
283          d = DOT(sd,vd);
284          if (d <= 0.0)
285                  return(-1.0);
286 +        if (d >= 1.0)
287 +                return(1.0);
288          sigma = acos(d) * (180./PI);
289 <        tau = acos(DOT(sd,vu)/sqrt(1.0-d*d)) * (180./PI);
289 >        d = DOT(sd,vu)/sqrt(1.0-d*d);
290 >        if (d >= 1.0)
291 >                tau = 0.0;
292 >        else
293 >                tau = acos(d) * (180./PI);
294          return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
295                          + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
296                  ) );
# Line 247 | Line 298 | FVECT  sd, vd, vu;
298  
299  
300   double
301 < gnuth_dgr(gd)           /* compute Gnuth discomfort glare rating */
301 > guth_dgr(gd)            /* compute Guth discomfort glare rating */
302   struct glare_dir        *gd;
303   {
304   #define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
305          register struct glare_src       *gs;
306 +        FVECT   mydir;
307          double  p;
308          double  sum;
309 +        double  wtot, brsum;
310          int     n;
311  
312 <        sum = 0.0; n = 0;
312 >        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
313 >        sum = wtot = brsum = 0.0; n = 0;
314          for (gs = all_srcs; gs != NULL; gs = gs->next) {
315 <                p = posindex(gs->dir, midview.vdir, midview.vup);
315 >                p = posindex(gs->dir, mydir, midview.vup);
316                  if (p <= FTINY)
317                          continue;
318                  sum += gs->lum * q(gs->dom) / p;
319 +                brsum += gs->lum * gs->dom;
320 +                wtot += gs->dom;
321                  n++;
322          }
323          if (n == 0)
324                  return(0.0);
325 <        else
326 <                return( pow(
271 <                        .5*sum/pow(direct(gd)+gd->indirect,.44),
325 >
326 >        return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
327                          pow((double)n, -.0914) ) );
328   #undef q
329   }
# Line 284 | Line 339 | extern double  erf(), erfc();
339  
340  
341   double
342 < gnuth_vcp(gd)           /* compute Gnuth visual comfort probability */
342 > guth_vcp(gd)            /* compute Guth visual comfort probability */
343   struct glare_dir        *gd;
344   {
345 <        return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd))));
345 >        double  dgr;
346 >
347 >        dgr = guth_dgr(gd);
348 >        if (dgr <= FTINY)
349 >                return(100.0);
350 >        return(100.*norm_integral(6.374-1.3227*log(dgr)));
351 > }
352 >
353 >
354 > double
355 > cie_cgi(gd)             /* compute CIE Glare Index */
356 > struct glare_dir        *gd;
357 > {
358 >        register struct glare_src       *gs;
359 >        FVECT   mydir;
360 >        double  dillum;
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 += gs->lum*gs->lum * gs->dom / (p*p);
371 >        }
372 >        if (sum <= FTINY)
373 >                return(0.0);
374 >        dillum = direct(gd);
375 >        return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
376   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines