--- ray/src/util/glarendx.c 1991/04/16 16:40:22 1.1 +++ ray/src/util/glarendx.c 1992/11/19 09:59:08 2.3 @@ -5,10 +5,14 @@ static char SCCSid[] = "$SunId$ LBL"; #endif /* - * Compute Glare Index given by program name: + * Compute Glare Index given by program name or -t option: * - * gnuth_dgr - Gnuth discomfort glare rating - * gnuth_vcp - Gnuth visual comfort probability + * guth_dgr - Guth discomfort glare rating + * guth_vcp - Guth visual comfort probability + * cie_cgi - CIE Glare Index (1983, due to Einhorn) + * vert_dir - Direct vertical illuminance + * vert_ind - Indirect vertical illuminance (from input) + * vert_ill - Total vertical illuminance * * 12 April 1991 Greg Ward EPFL */ @@ -19,17 +23,21 @@ static char SCCSid[] = "$SunId$ LBL"; extern double erfc(); double posindex(); -int headline(); -double direct(), gnuth_dgr(), gnuth_vcp(); +double direct(), guth_dgr(), guth_vcp(), cie_cgi(), + indirect(), total(); struct named_func { char *name; double (*func)(); + char *descrip; } all_funcs[] = { - {"direct", direct}, - {"gnuth_dgr", gnuth_dgr}, - {"gnuth_vcp", gnuth_vcp}, + {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"}, + {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"}, + {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"}, + {"vert_dir", direct, "Direct Vertical Illuminance"}, + {"vert_ill", total, "Total Vertical Illuminance"}, + {"vert_ind", indirect, "Indirect Vertical Illuminance"}, {NULL} }; @@ -53,7 +61,9 @@ int print_header = 1; VIEW midview = STDVIEW; +int wrongformat = 0; + main(argc, argv) int argc; char *argv[]; @@ -88,30 +98,24 @@ char *argv[]; perror(argv[i]); exit(1); } - /* read header */ - getheader(stdin, headline); - if (print_header) { /* add to header */ - printargs(i, argv, stdout); - putchar('\n'); - } - /* set view */ - if (setview(&midview) != NULL) { - fprintf(stderr, "%s: bad view information in input\n"); - exit(1); - } - /* get findglare data */ - read_input(); - /* find calculation */ + /* find and run calculation */ for (funp = all_funcs; funp->name != NULL; funp++) if (!strcmp(funp->name, progtail)) { + init(); + read_input(); + if (print_header) { + printargs(i, argv, stdout); + putchar('\n'); + } print_values(funp->func); exit(0); /* we're done */ } /* invalid function */ - fprintf(stderr, "%s: unknown function!\n", progtail); - exit(1); userr: - fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname); + fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname); + fprintf(stderr, "\twhere 'type' is one of the following:\n"); + for (funp = all_funcs; funp->name != NULL; funp++) + fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip); exit(1); } @@ -119,13 +123,36 @@ userr: headline(s) /* get line from header */ char *s; { + char fmt[32]; + if (print_header) /* copy to output */ fputs(s, stdout); - if (!strncmp(s, VIEWSTR, VIEWSTRL)) - sscanview(&midview, s+VIEWSTRL); + if (isview(s)) + sscanview(&midview, s); + else if (isformat(s)) { + formatval(fmt, s); + wrongformat = strcmp(fmt, "ascii"); + } } +init() /* initialize calculation */ +{ + /* read header */ + getheader(stdin, headline, NULL); + if (wrongformat) { + fprintf(stderr, "%s: bad input format\n", progname); + exit(1); + } + /* set view */ + if (setview(&midview) != NULL) { + fprintf(stderr, "%s: bad view information in input\n", + progname); + exit(1); + } +} + + read_input() /* read glare sources from stdin */ { #define S_SEARCH 0 @@ -155,6 +182,7 @@ read_input() /* read glare sources from stdin */ &gs->dir[0], &gs->dir[1], &gs->dir[2], &gs->dom, &gs->lum) != 5) goto readerr; + normalize(gs->dir); gs->next = all_srcs; all_srcs = gs; break; @@ -197,7 +225,7 @@ double (*funp)(); double -direct(gd) /* compute direct illuminance */ +direct(gd) /* compute direct vertical illuminance */ struct glare_dir *gd; { FVECT mydir; @@ -215,6 +243,22 @@ struct glare_dir *gd; } +double +indirect(gd) /* return indirect vertical illuminance */ +struct glare_dir *gd; +{ + return(gd->indirect); +} + + +double +total(gd) /* return total vertical illuminance */ +struct glare_dir *gd; +{ + return(direct(gd)+gd->indirect); +} + + /* * posindex - compute glare position index from: * @@ -225,6 +269,7 @@ struct glare_dir *gd; * All vectors are assumed to be normalized. * This function is an implementation of the method proposed by * Robert Levin in his 1975 JIES article. + * This calculation presumes the view direction and up vectors perpendicular. * We return a value less than zero for improper positions. */ @@ -238,8 +283,16 @@ FVECT sd, vd, vu; d = DOT(sd,vd); if (d <= 0.0) return(-1.0); + if (d >= 1.0) + return(1.0); sigma = acos(d) * (180./PI); - tau = acos(DOT(sd,vu)/sqrt(1.0-d*d)) * (180./PI); + d = DOT(sd,vu)/sqrt(1.0-d*d); + if (d >= 1.0) + tau = 0.0; + else if (d <= -1.0) + tau = 180.0; + else + tau = acos(d) * (180./PI); return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3 + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 ) ) ); @@ -247,28 +300,32 @@ FVECT sd, vd, vu; double -gnuth_dgr(gd) /* compute Gnuth discomfort glare rating */ +guth_dgr(gd) /* compute Guth discomfort glare rating */ struct glare_dir *gd; { #define q(w) (20.4*w+1.52*pow(w,.2)-.075) register struct glare_src *gs; + FVECT mydir; double p; double sum; + double wtot, brsum; int n; - sum = 0.0; n = 0; + spinvector(mydir, midview.vdir, midview.vup, gd->ang); + sum = wtot = brsum = 0.0; n = 0; for (gs = all_srcs; gs != NULL; gs = gs->next) { - p = posindex(gs->dir, midview.vdir, midview.vup); + p = posindex(gs->dir, mydir, midview.vup); if (p <= FTINY) continue; sum += gs->lum * q(gs->dom) / p; + brsum += gs->lum * gs->dom; + wtot += gs->dom; n++; } if (n == 0) return(0.0); - else - return( pow( - .5*sum/pow(direct(gd)+gd->indirect,.44), + + return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44), pow((double)n, -.0914) ) ); #undef q } @@ -284,8 +341,38 @@ extern double erf(), erfc(); double -gnuth_vcp(gd) /* compute Gnuth visual comfort probability */ +guth_vcp(gd) /* compute Guth visual comfort probability */ struct glare_dir *gd; { - return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd)))); + double dgr; + + dgr = guth_dgr(gd); + if (dgr <= FTINY) + return(100.0); + return(100.*norm_integral(6.374-1.3227*log(dgr))); +} + + +double +cie_cgi(gd) /* compute CIE Glare Index */ +struct glare_dir *gd; +{ + register struct glare_src *gs; + FVECT mydir; + double dillum; + double p; + double sum; + + spinvector(mydir, midview.vdir, midview.vup, gd->ang); + sum = 0.0; + for (gs = all_srcs; gs != NULL; gs = gs->next) { + p = posindex(gs->dir, mydir, midview.vup); + if (p <= FTINY) + continue; + sum += gs->lum*gs->lum * gs->dom / (p*p); + } + if (sum <= FTINY) + return(0.0); + dillum = direct(gd); + return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect))); }