Back to index

radiance  4R0+20100331
glarendx.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: glarendx.c,v 2.10 2004/01/02 12:48:36 schorsch Exp $";
00003 #endif
00004 /*
00005  * Compute Glare Index given by program name or -t option:
00006  *
00007  *     dgi -         Daylight Glare Index
00008  *     brs_gi -      Building Research Station Glare Index (Petherbridge
00009  *                                                      & Hopkinson)
00010  *     ugr -         Unified Glare Rating System (Fischer)
00011  *     guth_dgr -    Guth discomfort glare rating
00012  *     guth_vcp -    Guth visual comfort probability
00013  *     cie_cgi -     CIE Glare Index (1983, due to Einhorn)
00014  *     vert_dir -    Direct vertical illuminance
00015  *     vert_ind -    Indirect vertical illuminance (from input)
00016  *     vert_ill -    Total vertical illuminance
00017  *
00018  *            12 April 1991 Greg Ward     EPFL
00019  *            19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
00020  */
00021 
00022 #include <string.h>
00023 
00024 #include "standard.h"
00025 #include "view.h"
00026 
00027 
00028 struct glare_src {
00029        FVECT  dir;          /* source direction */
00030        double dom;          /* solid angle */
00031        double lum;          /* average luminance */
00032        struct glare_src     *next;
00033 } *all_srcs = NULL;
00034 
00035 struct glare_dir {
00036        double ang;          /* angle (in radians) */
00037        double indirect;     /* indirect illuminance */
00038        struct glare_dir     *next;
00039 } *all_dirs = NULL;
00040 
00041 typedef double gdfun(struct glare_dir *gd);
00042 
00043 static gdfun dgi;
00044 static gdfun brs_gi;
00045 static gdfun ugr;
00046 static gdfun guth_dgr;
00047 static gdfun guth_vcp;
00048 static gdfun cie_cgi;
00049 static gdfun direct;
00050 static gdfun indirect;
00051 static gdfun total;
00052 
00053 static gethfunc headline;
00054 static void init(void);
00055 static void read_input(void);
00056 static void print_values(gdfun *func);
00057 static double posindex(FVECT sd, FVECT vd, FVECT vu);
00058 
00059 
00060 struct named_func {
00061        char   *name;
00062        gdfun  *func;
00063        char   *descrip;
00064 } all_funcs[] = {
00065        {"dgi", dgi, "Daylight Glare Index"},
00066        {"brs_gi", brs_gi, "BRS Glare Index"},
00067        {"ugr", ugr, "Unified Glare Rating"},
00068        {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
00069        {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
00070        {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
00071        {"vert_dir", direct, "Direct Vertical Illuminance"},
00072        {"vert_ill", total, "Total Vertical Illuminance"},
00073        {"vert_ind", indirect, "Indirect Vertical Illuminance"},
00074        {NULL}
00075 };
00076 
00077 #define newp(type)   (type *)malloc(sizeof(type))
00078 
00079 char   *progname;
00080 int    print_header = 1;
00081 
00082 VIEW   midview = STDVIEW;
00083 
00084 int    wrongformat = 0;
00085 
00086 
00087 int
00088 main(
00089        int    argc,
00090        char   *argv[]
00091 )
00092 {
00093        struct named_func    *funp;
00094        char   *progtail;
00095        int    i;
00096                                    /* get program name */
00097        progname = argv[0];
00098        progtail = strrchr(progname, '/'); /* final component */
00099        if (progtail == NULL)
00100               progtail = progname;
00101        else
00102               progtail++;
00103                                    /* get options */
00104        for (i = 1; i < argc && argv[i][0] == '-'; i++)
00105               switch (argv[i][1]) {
00106               case 't':
00107                      progtail = argv[++i];
00108                      break;
00109               case 'h':
00110                      print_header = 0;
00111                      break;
00112               default:
00113                      goto userr;
00114               }
00115        if (i < argc-1)
00116               goto userr;
00117        if (i == argc-1)            /* open file */
00118               if (freopen(argv[i], "r", stdin) == NULL) {
00119                      perror(argv[i]);
00120                      exit(1);
00121               }
00122                                    /* find and run calculation */
00123        for (funp = all_funcs; funp->name != NULL; funp++)
00124               if (!strcmp(funp->name, progtail)) {
00125                      init();
00126                      read_input();
00127                      if (print_header) {
00128                             printargs(i, argv, stdout);
00129                             putchar('\n');
00130                      }
00131                      print_values(funp->func);
00132                      exit(0);             /* we're done */
00133               }
00134                                    /* invalid function */
00135 userr:
00136        fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
00137        fprintf(stderr, "\twhere 'type' is one of the following:\n");
00138        for (funp = all_funcs; funp->name != NULL; funp++)
00139               fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
00140        exit(1);
00141 }
00142 
00143 
00144 static int
00145 headline(                   /* get line from header */
00146        char   *s,
00147        void   *p
00148 )
00149 {
00150        char   fmt[32];
00151 
00152        if (print_header)           /* copy to output */
00153               fputs(s, stdout);
00154        if (isview(s))
00155               sscanview(&midview, s);
00156        else if (isformat(s)) {
00157               formatval(fmt, s);
00158               wrongformat = strcmp(fmt, "ascii");
00159        }
00160        return(0);
00161 }
00162 
00163 
00164 static void
00165 init(void)                         /* initialize calculation */
00166 {
00167                                    /* read header */
00168        getheader(stdin, headline, NULL);
00169        if (wrongformat) {
00170               fprintf(stderr, "%s: bad input format\n", progname);
00171               exit(1);
00172        }
00173                                    /* set view */
00174        if (setview(&midview) != NULL) {
00175               fprintf(stderr, "%s: bad view information in input\n",
00176                             progname);
00177               exit(1);
00178        }
00179 }
00180 
00181 
00182 static void
00183 read_input(void)                   /* read glare sources from stdin */
00184 {
00185 #define       S_SEARCH      0
00186 #define S_SOURCE     1
00187 #define S_DIREC             2
00188        int    state = S_SEARCH;
00189        char   buf[128];
00190        register struct glare_src   *gs;
00191        register struct glare_dir   *gd;
00192 
00193        while (fgets(buf, sizeof(buf), stdin) != NULL)
00194               switch (state) {
00195               case S_SEARCH:
00196                      if (!strcmp(buf, "BEGIN glare source\n"))
00197                             state = S_SOURCE;
00198                      else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
00199                             state = S_DIREC;
00200                      break;
00201               case S_SOURCE:
00202                      if (!strncmp(buf, "END", 3)) {
00203                             state = S_SEARCH;
00204                             break;
00205                      }
00206                      if ((gs = newp(struct glare_src)) == NULL)
00207                             goto memerr;
00208                      if (sscanf(buf, "%lf %lf %lf %lf %lf",
00209                                    &gs->dir[0], &gs->dir[1], &gs->dir[2],
00210                                    &gs->dom, &gs->lum) != 5)
00211                             goto readerr;
00212                      normalize(gs->dir);
00213                      gs->next = all_srcs;
00214                      all_srcs = gs;
00215                      break;
00216               case S_DIREC:
00217                      if (!strncmp(buf, "END", 3)) {
00218                             state = S_SEARCH;
00219                             break;
00220                      }
00221                      if ((gd = newp(struct glare_dir)) == NULL)
00222                             goto memerr;
00223                      if (sscanf(buf, "%lf %lf",
00224                                    &gd->ang, &gd->indirect) != 2)
00225                             goto readerr;
00226                      gd->ang *= PI/180.0; /* convert to radians */
00227                      gd->next = all_dirs;
00228                      all_dirs = gd;
00229                      break;
00230               }
00231        return;
00232 memerr:
00233        perror(progname);
00234        exit(1);
00235 readerr:
00236        fprintf(stderr, "%s: read error on input\n", progname);
00237        exit(1);
00238 #undef S_SEARCH
00239 #undef S_SOURCE
00240 #undef S_DIREC
00241 }
00242 
00243 
00244 static void
00245 print_values(        /* print out calculations */
00246        gdfun *func
00247 )
00248 {
00249        register struct glare_dir   *gd;
00250 
00251        for (gd = all_dirs; gd != NULL; gd = gd->next)
00252               printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
00253 }
00254 
00255 
00256 static double
00257 direct(                     /* compute direct vertical illuminance */
00258        struct glare_dir     *gd
00259 )
00260 {
00261        FVECT  mydir;
00262        double d, dval;
00263        register struct glare_src   *gs;
00264 
00265        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
00266        dval = 0.0;
00267        for (gs = all_srcs; gs != NULL; gs = gs->next) {
00268               d = DOT(mydir,gs->dir);
00269               if (d > FTINY)
00270                      dval += d * gs->dom * gs->lum;
00271        }
00272        return(dval);
00273 }
00274 
00275 
00276 static double
00277 indirect(                   /* return indirect vertical illuminance */
00278        struct glare_dir     *gd
00279 )
00280 {
00281        return(gd->indirect);
00282 }
00283 
00284 
00285 static double
00286 total(               /* return total vertical illuminance */
00287        struct glare_dir     *gd
00288 )
00289 {
00290        return(direct(gd)+gd->indirect);
00291 }
00292 
00293 
00294 /*
00295  * posindex - compute glare position index from:
00296  *
00297  *     Source Direction
00298  *     View Direction
00299  *     View Up Direction
00300  *
00301  * All vectors are assumed to be normalized.
00302  * This function is an implementation of the method proposed by
00303  * Robert Levin in his 1975 JIES article.
00304  * This calculation presumes the view direction and up vectors perpendicular.
00305  * We return a value less than zero for improper positions.
00306  */
00307 
00308 static double
00309 posindex(                   /* compute position index */
00310        FVECT  sd,
00311        FVECT  vd,
00312        FVECT  vu
00313 )
00314 {
00315        double sigma, tau;
00316        double d;
00317 
00318        d = DOT(sd,vd);
00319        if (d <= 0.0)
00320               return(-1.0);
00321        if (d >= 1.0)
00322               return(1.0);
00323        sigma = acos(d) * (180./PI);
00324        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
00325        if (d >= 1.0)
00326               tau = 0.0;
00327        else
00328               tau = acos(d) * (180./PI);
00329        return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
00330                      + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
00331               ) );
00332 }
00333 
00334 
00335 static double
00336 dgi(          /* compute Daylight Glare Index */
00337        struct glare_dir     *gd
00338 )
00339 {
00340        register struct glare_src   *gs;
00341        FVECT  mydir,testdir[7],vhor;
00342        double r,posn,omega,p[7],sum;
00343        int    i,n;
00344 
00345        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
00346        sum = 0.0; n = 0;
00347        for (gs = all_srcs; gs != NULL; gs = gs->next) {
00348 
00349               /* compute 1/p^2 weighted solid angle of the source */
00350               r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
00351               fcross(vhor,gs->dir,midview.vup);
00352               normalize(vhor);
00353               VCOPY(testdir[0],gs->dir);
00354               fvsum(testdir[1],gs->dir,vhor,r);
00355               fvsum(testdir[2],gs->dir,vhor,0.5*r);
00356               fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
00357               fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
00358               fvsum(testdir[3],gs->dir,vhor,-r);
00359               fvsum(testdir[4],gs->dir,vhor,-0.5*r);
00360               fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
00361               fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
00362               for (i = 0; i < 7; i++) {
00363                      normalize(testdir[i]);
00364                      posn = posindex(testdir[i],mydir,midview.vup);
00365                      if (posn <= FTINY)
00366                             p[i] = 0.0;
00367                      else
00368                             p[i] = 1./(posn*posn);
00369               }
00370               r = 1-gs->dom/2./PI;
00371               omega = gs->dom*p[0];
00372               omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
00373               omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
00374                        +0.3334*(p[2]+p[4]+p[5]+p[6]));
00375 
00376               sum += pow(gs->lum,1.6) * pow(omega,0.8) /
00377                      (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
00378               n++;
00379        }
00380        if (n == 0)
00381               return(0.0);
00382        return( 10*log10(0.478*sum) );
00383 }
00384 
00385 
00386 static double
00387 brs_gi(              /* compute BRS Glare Index */
00388        struct glare_dir     *gd
00389 )
00390 {
00391        register struct glare_src   *gs;
00392        FVECT  mydir;
00393        double p;
00394        double sum;
00395 
00396        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
00397        sum = 0.0;
00398        for (gs = all_srcs; gs != NULL; gs = gs->next) {
00399               p = posindex(gs->dir, mydir, midview.vup);
00400               if (p <= FTINY)
00401                      continue;
00402               sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
00403        }
00404        if (sum <= FTINY)
00405               return(0.0);
00406        sum /= gd->indirect/PI;
00407        return(10*log10(0.478*sum));
00408 }
00409 
00410 
00411 static double
00412 guth_dgr(            /* compute Guth discomfort glare rating */
00413        struct glare_dir     *gd
00414 )
00415 {
00416 #define q(w)  (20.4*w+1.52*pow(w,.2)-.075)
00417        register struct glare_src   *gs;
00418        FVECT  mydir;
00419        double p;
00420        double sum;
00421        double wtot, brsum;
00422        int    n;
00423 
00424        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
00425        sum = wtot = brsum = 0.0; n = 0;
00426        for (gs = all_srcs; gs != NULL; gs = gs->next) {
00427               p = posindex(gs->dir, mydir, midview.vup);
00428               if (p <= FTINY)
00429                      continue;
00430               sum += gs->lum * q(gs->dom) / p;
00431               brsum += gs->lum * gs->dom;
00432               wtot += gs->dom;
00433               n++;
00434        }
00435        if (n == 0)
00436               return(0.0);
00437        return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
00438                      pow((double)n, -.0914) ) );
00439 #undef q
00440 }
00441 
00442 
00443 #ifndef M_SQRT2
00444 #define       M_SQRT2       1.41421356237309504880
00445 #endif
00446 
00447 #define norm_integral(z)    (1.-.5*erfc((z)/M_SQRT2))
00448 
00449 
00450 static double
00451 guth_vcp(            /* compute Guth visual comfort probability */
00452        struct glare_dir     *gd
00453 )
00454 {
00455        extern double erfc();
00456        double dgr;
00457 
00458        dgr = guth_dgr(gd);
00459        if (dgr <= FTINY)
00460               return(100.0);
00461        return(100.*norm_integral(6.374-1.3227*log(dgr)));
00462 }
00463 
00464 
00465 static double
00466 cie_cgi(             /* compute CIE Glare Index */
00467        struct glare_dir     *gd
00468 )
00469 {
00470        register struct glare_src   *gs;
00471        FVECT  mydir;
00472        double dillum;
00473        double p;
00474        double sum;
00475 
00476        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
00477        sum = 0.0;
00478        for (gs = all_srcs; gs != NULL; gs = gs->next) {
00479               p = posindex(gs->dir, mydir, midview.vup);
00480               if (p <= FTINY)
00481                      continue;
00482               sum += gs->lum*gs->lum * gs->dom / (p*p);
00483        }
00484        if (sum <= FTINY)
00485               return(0.0);
00486        dillum = direct(gd);
00487        return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
00488 }
00489 
00490 
00491 static double
00492 ugr(          /* compute Unified Glare Rating */
00493        struct glare_dir     *gd
00494 )
00495 {
00496        register struct glare_src   *gs;
00497        FVECT  mydir;
00498        double p;
00499        double sum;
00500 
00501        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
00502        sum = 0.0;
00503        for (gs = all_srcs; gs != NULL; gs = gs->next) {
00504               p = posindex(gs->dir, mydir, midview.vup);
00505               if (p <= FTINY)
00506                      continue;
00507               sum += gs->lum*gs->lum * gs->dom / (p*p);
00508        }
00509        if (sum <= FTINY)
00510               return(0.0);
00511        return(8.*log10(0.25*sum*PI/gd->indirect));
00512 }