--- ray/src/util/glarendx.c 1991/04/16 16:40:22 1.1 +++ ray/src/util/glarendx.c 2003/02/22 02:07:30 2.7 @@ -1,59 +1,73 @@ -/* Copyright (c) 1991 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: glarendx.c,v 2.7 2003/02/22 02:07:30 greg Exp $"; #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 + * dgi - Daylight Glare Index + * brs_gi - Building Research Station Glare Index (Petherbridge + * & Hopkinson) + * ugr - Unified Glare Rating System (Fischer) + * 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 April 1993 R. Compagnon EPFL (added dgi, brs_gi, ugr) */ - + #include "standard.h" #include "view.h" - -extern double erfc(); - + + double posindex(); -int headline(); - -double direct(), gnuth_dgr(), gnuth_vcp(); - + +double direct(), guth_dgr(), guth_vcp(), cie_cgi(), + indirect(), total(), dgi(), brs_gi(), ugr(); + struct named_func { char *name; double (*func)(); + char *descrip; } all_funcs[] = { - {"direct", direct}, - {"gnuth_dgr", gnuth_dgr}, - {"gnuth_vcp", gnuth_vcp}, + {"dgi", dgi, "Daylight Glare Index"}, + {"brs_gi", brs_gi, "BRS Glare Index"}, + {"ugr", ugr, "Unified Glare Rating"}, + {"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} }; - + struct glare_src { FVECT dir; /* source direction */ double dom; /* solid angle */ double lum; /* average luminance */ struct glare_src *next; } *all_srcs = NULL; - + struct glare_dir { double ang; /* angle (in radians) */ double indirect; /* indirect illuminance */ struct glare_dir *next; } *all_dirs = NULL; - + #define newp(type) (type *)malloc(sizeof(type)) - + char *progname; int print_header = 1; - + VIEW midview = STDVIEW; - - + +int wrongformat = 0; + + main(argc, argv) int argc; char *argv[]; @@ -88,44 +102,63 @@ 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); } - - + + +int 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"); + } + return(0); } - - + + +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 @@ -135,7 +168,7 @@ read_input() /* read glare sources from stdin */ char buf[128]; register struct glare_src *gs; register struct glare_dir *gd; - + while (fgets(buf, sizeof(buf), stdin) != NULL) switch (state) { case S_SEARCH: @@ -155,6 +188,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; @@ -184,26 +218,26 @@ readerr: #undef S_SOURCE #undef S_DIREC } - - + + print_values(funp) /* print out calculations */ double (*funp)(); { register struct glare_dir *gd; - + for (gd = all_dirs; gd != NULL; gd = gd->next) printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd)); } - - + + double -direct(gd) /* compute direct illuminance */ +direct(gd) /* compute direct vertical illuminance */ struct glare_dir *gd; { FVECT mydir; double d, dval; register struct glare_src *gs; - + spinvector(mydir, midview.vdir, midview.vup, gd->ang); dval = 0.0; for (gs = all_srcs; gs != NULL; gs = gs->next) { @@ -213,8 +247,24 @@ struct glare_dir *gd; } return(dval); } - - + + +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,67 +275,200 @@ 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. */ - + double posindex(sd, vd, vu) /* compute position index */ FVECT sd, vd, vu; { double sigma, tau; double d; - + 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 = fabs(DOT(sd,vu)/sqrt(1.0-d*d)); + if (d >= 1.0) + tau = 0.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 ) ) ); } + + +double +dgi(gd) /* compute Daylight Glare Index */ +struct glare_dir *gd; +{ + register struct glare_src *gs; + FVECT mydir,testdir[7],vhor; + double r,omega,p[7],sum; + int i,n; + + spinvector(mydir, midview.vdir, midview.vup, gd->ang); + sum = 0.0; n = 0; + for (gs = all_srcs; gs != NULL; gs = gs->next) { + /* compute 1/p^2 weighted solid angle of the source */ + r = sqrt(1 - pow(1.-gs->dom/2./PI,2.)); + fcross(vhor,gs->dir,midview.vup); + normalize(vhor); + VCOPY(testdir[0],gs->dir); + fvsum(testdir[1],gs->dir,vhor,r); + fvsum(testdir[2],gs->dir,vhor,0.5*r); + fvsum(testdir[5],testdir[2],midview.vup,-0.866*r); + fvsum(testdir[2],testdir[2],midview.vup,0.866*r); + fvsum(testdir[3],gs->dir,vhor,-r); + fvsum(testdir[4],gs->dir,vhor,-0.5*r); + fvsum(testdir[6],testdir[4],midview.vup,0.866*r); + fvsum(testdir[4],testdir[4],midview.vup,-0.866*r); + for (i = 0; i < 7; i++) { + normalize(testdir[i]); + p[i] = pow(posindex(testdir[i],mydir,midview.vup),-2.0); + if (p[i] <= FTINY) p[i] = 0.0; + } + r = 1-gs->dom/2./PI; + omega = gs->dom*p[0]; + omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5); + omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3]) + +0.3334*(p[2]+p[4]+p[5]+p[6])); + sum += pow(gs->lum,1.6) * pow(omega,0.8) / + (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum); + n++; + } + if (n == 0) + return(0.0); + return( 10*log10(0.478*sum) ); +} + + double -gnuth_dgr(gd) /* compute Gnuth discomfort glare rating */ +brs_gi(gd) /* compute BRS Glare Index */ struct glare_dir *gd; { + register struct glare_src *gs; + FVECT mydir; + 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 += pow(gs->lum/p,1.6) * pow(gs->dom,0.8); + } + if (sum <= FTINY) + return(0.0); + sum /= gd->indirect/PI; + return(10*log10(0.478*sum)); +} + + +double +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 } - - -extern double erf(), erfc(); - + + #ifndef M_SQRT2 #define M_SQRT2 1.41421356237309504880 #endif - + #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2)) + + +double +guth_vcp(gd) /* compute Guth visual comfort probability */ +struct glare_dir *gd; +{ + extern double erfc(); + 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))); +} double -gnuth_vcp(gd) /* compute Gnuth visual comfort probability */ +ugr(gd) /* compute Unified Glare Rating */ struct glare_dir *gd; { - return(100.*norm_integral(-6.374+1.3227*log(gnuth_dgr(gd)))); + register struct glare_src *gs; + FVECT mydir; + 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); + return(8.*log10(0.25*sum*PI/gd->indirect)); }