ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 1.6
Committed: Mon Apr 22 10:11:21 1991 UTC (32 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.5: +1 -1 lines
Log Message:
ASCII -> ascii

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