ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 1.3
Committed: Wed Apr 17 11:55:03 1991 UTC (32 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.2: +6 -3 lines
Log Message:
corrected calculation of Guth F (average field luminance)

File Contents

# User Rev Content
1 greg 1.1 /* Copyright (c) 1991 Regents of the University of California */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * Compute Glare Index given by program name:
9     *
10 greg 1.2 * guth_dgr - Guth discomfort glare rating
11     * guth_vcp - Guth visual comfort probability
12 greg 1.1 *
13     * 12 April 1991 Greg Ward EPFL
14     */
15    
16     #include "standard.h"
17     #include "view.h"
18    
19     extern double erfc();
20    
21     double posindex();
22     int headline();
23    
24 greg 1.2 double direct(), guth_dgr(), guth_vcp();
25 greg 1.1
26     struct named_func {
27     char *name;
28     double (*func)();
29     } all_funcs[] = {
30     {"direct", direct},
31 greg 1.2 {"guth_dgr", guth_dgr},
32     {"guth_vcp", guth_vcp},
33 greg 1.1 {NULL}
34     };
35    
36     struct glare_src {
37     FVECT dir; /* source direction */
38     double dom; /* solid angle */
39     double lum; /* average luminance */
40     struct glare_src *next;
41     } *all_srcs = NULL;
42    
43     struct glare_dir {
44     double ang; /* angle (in radians) */
45     double indirect; /* indirect illuminance */
46     struct glare_dir *next;
47     } *all_dirs = NULL;
48    
49     #define newp(type) (type *)malloc(sizeof(type))
50    
51     char *progname;
52     int print_header = 1;
53    
54     VIEW midview = STDVIEW;
55    
56    
57     main(argc, argv)
58     int argc;
59     char *argv[];
60     {
61     extern char *rindex();
62     struct named_func *funp;
63     char *progtail;
64     int i;
65     /* get program name */
66     progname = argv[0];
67     progtail = rindex(progname, '/'); /* final component */
68     if (progtail == NULL)
69     progtail = progname;
70     else
71     progtail++;
72     /* get options */
73     for (i = 1; i < argc && argv[i][0] == '-'; i++)
74     switch (argv[i][1]) {
75     case 't':
76     progtail = argv[++i];
77     break;
78     case 'h':
79     print_header = 0;
80     break;
81     default:
82     goto userr;
83     }
84     if (i < argc-1)
85     goto userr;
86     if (i == argc-1) /* open file */
87     if (freopen(argv[i], "r", stdin) == NULL) {
88     perror(argv[i]);
89     exit(1);
90     }
91     /* read header */
92     getheader(stdin, headline);
93     if (print_header) { /* add to header */
94     printargs(i, argv, stdout);
95     putchar('\n');
96     }
97     /* set view */
98     if (setview(&midview) != NULL) {
99     fprintf(stderr, "%s: bad view information in input\n");
100     exit(1);
101     }
102     /* get findglare data */
103     read_input();
104     /* find calculation */
105     for (funp = all_funcs; funp->name != NULL; funp++)
106     if (!strcmp(funp->name, progtail)) {
107     print_values(funp->func);
108     exit(0); /* we're done */
109     }
110     /* invalid function */
111     fprintf(stderr, "%s: unknown function!\n", progtail);
112     exit(1);
113     userr:
114     fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
115     exit(1);
116     }
117    
118    
119     headline(s) /* get line from header */
120     char *s;
121     {
122     if (print_header) /* copy to output */
123     fputs(s, stdout);
124     if (!strncmp(s, VIEWSTR, VIEWSTRL))
125     sscanview(&midview, s+VIEWSTRL);
126     }
127    
128    
129     read_input() /* read glare sources from stdin */
130     {
131     #define S_SEARCH 0
132     #define S_SOURCE 1
133     #define S_DIREC 2
134     int state = S_SEARCH;
135     char buf[128];
136     register struct glare_src *gs;
137     register struct glare_dir *gd;
138    
139     while (fgets(buf, sizeof(buf), stdin) != NULL)
140     switch (state) {
141     case S_SEARCH:
142     if (!strcmp(buf, "BEGIN glare source\n"))
143     state = S_SOURCE;
144     else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
145     state = S_DIREC;
146     break;
147     case S_SOURCE:
148     if (!strncmp(buf, "END", 3)) {
149     state = S_SEARCH;
150     break;
151     }
152     if ((gs = newp(struct glare_src)) == NULL)
153     goto memerr;
154     if (sscanf(buf, "%lf %lf %lf %lf %lf",
155     &gs->dir[0], &gs->dir[1], &gs->dir[2],
156     &gs->dom, &gs->lum) != 5)
157     goto readerr;
158 greg 1.2 normalize(gs->dir);
159 greg 1.1 gs->next = all_srcs;
160     all_srcs = gs;
161     break;
162     case S_DIREC:
163     if (!strncmp(buf, "END", 3)) {
164     state = S_SEARCH;
165     break;
166     }
167     if ((gd = newp(struct glare_dir)) == NULL)
168     goto memerr;
169     if (sscanf(buf, "%lf %lf",
170     &gd->ang, &gd->indirect) != 2)
171     goto readerr;
172     gd->ang *= PI/180.0; /* convert to radians */
173     gd->next = all_dirs;
174     all_dirs = gd;
175     break;
176     }
177     return;
178     memerr:
179     perror(progname);
180     exit(1);
181     readerr:
182     fprintf(stderr, "%s: read error on input\n", progname);
183     exit(1);
184     #undef S_SEARCH
185     #undef S_SOURCE
186     #undef S_DIREC
187     }
188    
189    
190     print_values(funp) /* print out calculations */
191     double (*funp)();
192     {
193     register struct glare_dir *gd;
194    
195     for (gd = all_dirs; gd != NULL; gd = gd->next)
196     printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
197     }
198    
199    
200     double
201     direct(gd) /* compute direct illuminance */
202     struct glare_dir *gd;
203     {
204     FVECT mydir;
205     double d, dval;
206     register struct glare_src *gs;
207    
208     spinvector(mydir, midview.vdir, midview.vup, gd->ang);
209     dval = 0.0;
210     for (gs = all_srcs; gs != NULL; gs = gs->next) {
211     d = DOT(mydir,gs->dir);
212     if (d > FTINY)
213     dval += d * gs->dom * gs->lum;
214     }
215     return(dval);
216     }
217    
218    
219     /*
220     * posindex - compute glare position index from:
221     *
222     * Source Direction
223     * View Direction
224     * View Up Direction
225     *
226     * All vectors are assumed to be normalized.
227     * This function is an implementation of the method proposed by
228     * Robert Levin in his 1975 JIES article.
229     * We return a value less than zero for improper positions.
230     */
231    
232     double
233     posindex(sd, vd, vu) /* compute position index */
234     FVECT sd, vd, vu;
235     {
236     double sigma, tau;
237     double d;
238    
239     d = DOT(sd,vd);
240     if (d <= 0.0)
241     return(-1.0);
242     sigma = acos(d) * (180./PI);
243     tau = acos(DOT(sd,vu)/sqrt(1.0-d*d)) * (180./PI);
244     return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
245     + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
246     ) );
247     }
248    
249    
250     double
251 greg 1.2 guth_dgr(gd) /* compute Guth discomfort glare rating */
252 greg 1.1 struct glare_dir *gd;
253     {
254     #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
255     register struct glare_src *gs;
256     double p;
257     double sum;
258 greg 1.3 double wtot, brsum;
259 greg 1.1 int n;
260    
261 greg 1.3 sum = wtot = brsum = 0.0; n = 0;
262 greg 1.1 for (gs = all_srcs; gs != NULL; gs = gs->next) {
263     p = posindex(gs->dir, midview.vdir, midview.vup);
264     if (p <= FTINY)
265     continue;
266     sum += gs->lum * q(gs->dom) / p;
267 greg 1.3 brsum += gs->lum * gs->dom;
268     wtot += gs->dom;
269 greg 1.1 n++;
270     }
271     if (n == 0)
272     return(0.0);
273     else
274     return( pow(
275 greg 1.3 .5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
276 greg 1.1 pow((double)n, -.0914) ) );
277     #undef q
278     }
279    
280    
281     extern double erf(), erfc();
282    
283     #ifndef M_SQRT2
284     #define M_SQRT2 1.41421356237309504880
285     #endif
286    
287     #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
288    
289    
290     double
291 greg 1.2 guth_vcp(gd) /* compute Guth visual comfort probability */
292 greg 1.1 struct glare_dir *gd;
293     {
294 greg 1.3 return(100.*norm_integral(6.374-1.3227*log(guth_dgr(gd))));
295 greg 1.1 }