ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/glarendx.c
(Generate patch)

Comparing ray/src/util/glarendx.c (file contents):
Revision 1.8 by greg, Thu May 2 13:17:41 1991 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:
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)
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 "paths.h"
26   #include "view.h"
27  
20 – extern double   erfc();
28  
22 – double  posindex();
23 – int     headline();
24 –
25 – double  direct(), guth_dgr(), guth_vcp(), cie_cgi();
26 –
27 – struct named_func {
28 –        char    *name;
29 –        double  (*func)();
30 – } all_funcs[] = {
31 –        {"direct", direct},
32 –        {"guth_dgr", guth_dgr},
33 –        {"guth_vcp", guth_vcp},
34 –        {"cie_cgi", cie_cgi},
35 –        {NULL}
36 – };
37 –
29   struct glare_src {
30          FVECT   dir;            /* source direction */
31          double  dom;            /* solid angle */
# Line 48 | Line 39 | struct glare_dir {
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 +        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"},
72 +        {"vert_dir", direct, "Direct Vertical Illuminance"},
73 +        {"vert_ill", total, "Total Vertical Illuminance"},
74 +        {"vert_ind", indirect, "Indirect Vertical Illuminance"},
75 +        {NULL}
76 + };
77 +
78   #define newp(type)      (type *)malloc(sizeof(type))
79  
80   char    *progname;
# Line 58 | Line 85 | VIEW   midview = STDVIEW;
85   int     wrongformat = 0;
86  
87  
88 < main(argc, argv)
89 < int     argc;
90 < char    *argv[];
88 > int
89 > main(
90 >        int     argc,
91 >        char    *argv[]
92 > )
93   {
65 –        extern char     *rindex();
94          struct named_func       *funp;
67 –        char    *progtail;
95          int     i;
96                                          /* get program name */
97 <        progname = argv[0];
71 <        progtail = rindex(progname, '/');       /* final component */
72 <        if (progtail == NULL)
73 <                progtail = progname;
74 <        else
75 <                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 92 | Line 114 | char   *argv[];
114                          perror(argv[i]);
115                          exit(1);
116                  }
117 <                                        /* read header */
96 <        getheader(stdin, headline, NULL);
97 <        if (wrongformat) {
98 <                fprintf(stderr, "%s: bad input format\n", progname);
99 <                exit(1);
100 <        }
101 <        if (print_header) {             /* add to header */
102 <                printargs(i, argv, stdout);
103 <                putchar('\n');
104 <        }
105 <                                        /* set view */
106 <        if (setview(&midview) != NULL) {
107 <                fprintf(stderr, "%s: bad view information in input\n",
108 <                                progname);
109 <                exit(1);
110 <        }
111 <                                        /* get findglare data */
112 <        read_input();
113 <                                        /* find calculation */
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) {
123 >                                printargs(i, argv, stdout);
124 >                                putchar('\n');
125 >                        }
126                          print_values(funp->func);
127                          exit(0);                /* we're done */
128                  }
129                                          /* invalid function */
120 –        fprintf(stderr, "%s: unknown function!\n", progtail);
121 –        exit(1);
130   userr:
131 <        fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
131 >        fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
132 >        fprintf(stderr, "\twhere 'type' is one of the following:\n");
133 >        for (funp = all_funcs; funp->name != NULL; funp++)
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;
139 > static int
140 > headline(                       /* get line from header */
141 >        char    *s,
142 >        void    *p
143 > )
144   {
145 <        char    fmt[32];
145 >        char    fmt[MAXFMTLEN];
146  
147          if (print_header)               /* copy to output */
148                  fputs(s, stdout);
149 <        if (!strncmp(s, VIEWSTR, VIEWSTRL))
150 <                sscanview(&midview, s+VIEWSTRL);
149 >        if (isview(s))
150 >                sscanview(&midview, s);
151          else if (isformat(s)) {
152                  formatval(fmt, s);
153                  wrongformat = strcmp(fmt, "ascii");
154          }
155 +        return(0);
156   }
157  
158  
159 < read_input()                    /* read glare sources from stdin */
159 > static void
160 > init(void)                              /* initialize calculation */
161   {
162 +        char    *err;
163 +                                        /* read header */
164 +        getheader(stdin, headline, NULL);
165 +        if (wrongformat) {
166 +                fprintf(stderr, "%s: bad input format\n", progname);
167 +                exit(1);
168 +        }
169 +                                        /* set view */
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 + 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;
192 >        struct glare_src        *gs;
193 >        struct glare_dir        *gd;
194  
195          while (fgets(buf, sizeof(buf), stdin) != NULL)
196                  switch (state) {
# Line 166 | 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 202 | Line 244 | readerr:
244   }
245  
246  
247 < print_values(funp)              /* print out calculations */
248 < double  (*funp)();
247 > static void
248 > print_values(           /* print out calculations */
249 >        gdfun *func
250 > )
251   {
252 <        register struct glare_dir       *gd;
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 illuminance */
261 < struct glare_dir        *gd;
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;
266 >        struct glare_src        *gs;
267  
268          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
269          dval = 0.0;
# Line 231 | Line 276 | struct glare_dir       *gd;
276   }
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 + static double
289 + total(                  /* return total vertical illuminance */
290 +        struct glare_dir        *gd
291 + )
292 + {
293 +        return(direct(gd)+gd->indirect);
294 + }
295 +
296 +
297   /*
298   * posindex -   compute glare position index from:
299   *
# Line 245 | Line 308 | struct glare_dir       *gd;
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;
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;
# Line 255 | Line 321 | FVECT  sd, vd, vu;
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 = DOT(sd,vu)/sqrt(1.0-d*d);
328 <        if (d >= 1.0)
263 <                tau = 0.0;
264 <        else
265 <                tau = acos(d) * (180./PI);
327 >        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
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 < guth_dgr(gd)            /* compute Guth discomfort glare rating */
337 < struct glare_dir        *gd;
335 > static double
336 > dgi(            /* compute Daylight Glare Index */
337 >        struct glare_dir        *gd
338 > )
339   {
340 +        struct glare_src        *gs;
341 +        FVECT   mydir,testdir[7],vhor;
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;
347 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
348 +
349 +                /* compute 1/p^2 weighted solid angle of the source */
350 +                r = 1. - gs->dom*(.5/PI);
351 +                r = sqrt(1. - r*r);
352 +                fcross(vhor,gs->dir,midview.vup);
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);
358 +                fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
359 +                fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
360 +                fvsum(testdir[3],gs->dir,vhor,-r);
361 +                fvsum(testdir[4],gs->dir,vhor,-0.5*r);
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 +                        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*(.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])
376 +                          +0.3334*(p[2]+p[4]+p[5]+p[6]));
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 (sum <= FTINY)
383 +                return(0.0);
384 +        return( 10.*log10(0.478*sum) );
385 + }
386 +
387 +
388 + static double
389 + brs_gi(         /* compute BRS Glare Index */
390 +        struct glare_dir        *gd
391 + )
392 + {
393 +        struct glare_src        *gs;
394 +        FVECT   mydir;
395 +        double  p;
396 +        double  sum;
397 +
398 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
399 +        sum = 0.0;
400 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
401 +                p = posindex(gs->dir, mydir, midview.vup);
402 +                if (p <= FTINY)
403 +                        continue;
404 +                sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
405 +        }
406 +        if (sum <= FTINY)
407 +                return(0.0);
408 +        sum *= PI/gd->indirect;
409 +        return(10*log10(0.478*sum));
410 + }
411 +
412 +
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;
# Line 294 | Line 436 | struct glare_dir       *gd;
436          }
437          if (n == 0)
438                  return(0.0);
297 –
439          return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
440                          pow((double)n, -.0914) ) );
441   #undef q
442   }
443  
444  
304 – extern double   erf(), erfc();
305 –
445   #ifndef M_SQRT2
446   #define M_SQRT2 1.41421356237309504880
447   #endif
# Line 310 | Line 449 | extern double  erf(), erfc();
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;
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  
460          dgr = guth_dgr(gd);
# Line 323 | Line 464 | struct glare_dir       *gd;
464   }
465  
466  
467 < double
468 < cie_cgi(gd)             /* compute CIE Glare Index */
469 < struct glare_dir        *gd;
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;
# Line 345 | Line 487 | struct glare_dir       *gd;
487                  return(0.0);
488          dillum = direct(gd);
489          return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
490 + }
491 +
492 +
493 + static double
494 + ugr(            /* compute Unified Glare Rating */
495 +        struct glare_dir        *gd
496 + )
497 + {
498 +        struct glare_src        *gs;
499 +        FVECT   mydir;
500 +        double  p;
501 +        double  sum;
502 +
503 +        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
504 +        sum = 0.0;
505 +        for (gs = all_srcs; gs != NULL; gs = gs->next) {
506 +                p = posindex(gs->dir, mydir, midview.vup);
507 +                if (p <= FTINY)
508 +                        continue;
509 +                sum += gs->lum*gs->lum * gs->dom / (p*p);
510 +        }
511 +        if (sum <= FTINY)
512 +                return(0.0);
513 +        return(8.*log10(0.25*sum*PI/gd->indirect));
514   }

Diff Legend

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