--- ray/src/util/glarendx.c 1991/04/29 08:31:23 1.7 +++ ray/src/util/glarendx.c 2003/11/19 16:21:28 2.9 @@ -1,72 +1,85 @@ -/* 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.9 2003/11/19 16:21:28 greg Exp $"; #endif - /* - * Compute Glare Index given by program name: + * Compute Glare Index given by program name or -t option: * + * 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 #include "standard.h" #include "view.h" - -extern double erfc(); - + + double posindex(); -int headline(); - -double direct(), guth_dgr(), guth_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}, - {"guth_dgr", guth_dgr}, - {"guth_vcp", guth_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[]; { - extern char *rindex(); struct named_func *funp; char *progtail; int i; /* get program name */ progname = argv[0]; - progtail = rindex(progname, '/'); /* final component */ + progtail = strrchr(progname, '/'); /* final component */ if (progtail == NULL) progtail = progname; else @@ -90,55 +103,63 @@ char *argv[]; perror(argv[i]); exit(1); } - /* read header */ - getheader(stdin, headline, NULL); - if (wrongformat) { - fprintf(stderr, "%s: bad input format\n", progname); - exit(1); - } - 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", - progname); - 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 @@ -148,7 +169,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: @@ -198,26 +219,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) { @@ -227,8 +248,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: * @@ -242,21 +279,21 @@ struct glare_dir *gd; * 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); - d = DOT(sd,vu)/sqrt(1.0-d*d); + d = fabs(DOT(sd,vu)/sqrt(1.0-d*d)); if (d >= 1.0) tau = 0.0; else @@ -265,9 +302,83 @@ FVECT sd, vd, vu; + 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,posn,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]); + posn = posindex(testdir[i],mydir,midview.vup); + if (posn <= FTINY) + p[i] = 0.0; + else + p[i] = 1./(posn*posn); + } + 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 +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; { @@ -278,7 +389,7 @@ struct glare_dir *gd; double sum; double wtot, brsum; int n; - + spinvector(mydir, midview.vdir, midview.vup, gd->ang); sum = wtot = brsum = 0.0; n = 0; for (gs = all_srcs; gs != NULL; gs = gs->next) { @@ -292,31 +403,76 @@ struct glare_dir *gd; } if (n == 0) return(0.0); - else - return( pow( - .5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.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 +ugr(gd) /* compute Unified Glare Rating */ +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 += gs->lum*gs->lum * gs->dom / (p*p); + } + if (sum <= FTINY) + return(0.0); + return(8.*log10(0.25*sum*PI/gd->indirect)); }