ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.20
Committed: Sat Feb 22 02:07:23 2003 UTC (21 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 2.19: +5 -5 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.20 static const char RCSid[] = "$Id$";
3 greg 1.1 #endif
4     /*
5     * gensky.c - program to generate sky functions.
6     * Our zenith is along the Z-axis, the X-axis
7     * points east, and the Y-axis points north.
8     * Radiance is in watts/steradian/sq. meter.
9     *
10     * 3/26/86
11     */
12    
13     #include <stdio.h>
14    
15 greg 2.20 #include <stdlib.h>
16    
17     #include <string.h>
18    
19 greg 1.1 #include <math.h>
20    
21 greg 2.17 #include <ctype.h>
22    
23 greg 1.6 #include "color.h"
24 greg 1.1
25 greg 2.17 extern double stadj(), sdec(), sazi(), salt(), tz2mer();
26 greg 1.1
27 greg 2.17 #ifndef PI
28     #define PI 3.14159265358979323846
29     #endif
30 greg 1.1
31     #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
32    
33 greg 2.13 #define S_CLEAR 1
34     #define S_OVER 2
35     #define S_UNIF 3
36     #define S_INTER 4
37    
38     #define overcast (skytype==S_OVER|skytype==S_UNIF)
39    
40 greg 1.1 double normsc();
41     /* sun calculation constants */
42     extern double s_latitude;
43     extern double s_longitude;
44     extern double s_meridian;
45 greg 2.17
46     #undef toupper
47     #define toupper(c) ((c) & ~0x20) /* ASCII trick to convert case */
48    
49     /* European and North American zones */
50     struct {
51     char zname[8]; /* time zone name (all caps) */
52     float zmer; /* standard meridian */
53     } tzone[] = {
54     "YST", 135, "YDT", 120,
55     "PST", 120, "PDT", 105,
56     "MST", 105, "MDT", 90,
57     "CST", 90, "CDT", 75,
58     "EST", 75, "EDT", 60,
59     "AST", 60, "ADT", 45,
60     "NST", 52.5, "NDT", 37.5,
61     "GMT", 0, "BST", -15,
62 gwlarson 2.19 "CET", -15, "CEST", -30,
63     "EET", -30, "EEST", -45,
64     "AST", -45, "ADT", -60,
65     "GST", -60, "GDT", -75,
66     "IST", -82.5, "IDT", -97.5,
67     "JST", -135, "NDT", -150,
68     "NZST", -180, "NZDT", -195,
69 greg 2.17 "", 0
70     };
71 greg 1.1 /* required values */
72 greg 2.3 int month, day; /* date */
73 greg 2.5 double hour; /* time */
74     int tsolar; /* 0=standard, 1=solar */
75 greg 2.3 double altitude, azimuth; /* or solar angles */
76 greg 1.1 /* default values */
77 greg 2.13 int skytype = S_CLEAR; /* sky type */
78 greg 1.1 int dosun = 1;
79 greg 2.12 double zenithbr = 0.0;
80     int u_zenith = 0; /* -1=irradiance, 1=radiance */
81 greg 1.1 double turbidity = 2.75;
82     double gprefl = 0.2;
83     /* computed values */
84     double sundir[3];
85     double groundbr;
86     double F2;
87 greg 2.12 double solarbr = 0.0;
88     int u_solar = 0; /* -1=irradiance, 1=radiance */
89 greg 1.1
90     char *progname;
91     char errmsg[128];
92    
93    
94     main(argc, argv)
95     int argc;
96     char *argv[];
97     {
98     int i;
99    
100     progname = argv[0];
101     if (argc == 2 && !strcmp(argv[1], "-defaults")) {
102     printdefaults();
103     exit(0);
104     }
105     if (argc < 4)
106     userror("arg count");
107 greg 2.3 if (!strcmp(argv[1], "-ang")) {
108     altitude = atof(argv[2]) * (PI/180);
109     azimuth = atof(argv[3]) * (PI/180);
110     month = 0;
111     } else {
112     month = atoi(argv[1]);
113 greg 2.6 if (month < 1 || month > 12)
114     userror("bad month");
115 greg 2.3 day = atoi(argv[2]);
116 greg 2.6 if (day < 1 || day > 31)
117     userror("bad day");
118 greg 2.16 cvthour(argv[3]);
119 greg 2.3 }
120 greg 1.1 for (i = 4; i < argc; i++)
121     if (argv[i][0] == '-' || argv[i][0] == '+')
122     switch (argv[i][1]) {
123     case 's':
124 greg 2.13 skytype = S_CLEAR;
125 greg 1.1 dosun = argv[i][0] == '+';
126     break;
127 greg 2.8 case 'r':
128 greg 2.12 case 'R':
129     u_solar = argv[i][1]=='R' ? -1 : 1;
130 greg 2.8 solarbr = atof(argv[++i]);
131     break;
132 greg 1.1 case 'c':
133 greg 2.13 skytype = S_OVER;
134 greg 1.1 break;
135 greg 2.13 case 'u':
136     skytype = S_UNIF;
137     break;
138     case 'i':
139     skytype = S_INTER;
140     dosun = argv[i][0] == '+';
141     break;
142 greg 1.1 case 't':
143     turbidity = atof(argv[++i]);
144     break;
145     case 'b':
146 greg 2.12 case 'B':
147     u_zenith = argv[i][1]=='B' ? -1 : 1;
148 greg 1.1 zenithbr = atof(argv[++i]);
149     break;
150     case 'g':
151     gprefl = atof(argv[++i]);
152     break;
153     case 'a':
154     s_latitude = atof(argv[++i]) * (PI/180);
155     break;
156     case 'o':
157     s_longitude = atof(argv[++i]) * (PI/180);
158     break;
159     case 'm':
160     s_meridian = atof(argv[++i]) * (PI/180);
161     break;
162     default:
163     sprintf(errmsg, "unknown option: %s", argv[i]);
164     userror(errmsg);
165     }
166     else
167     userror("bad option");
168 greg 1.5
169 greg 2.18 if (fabs(s_meridian-s_longitude) > 45*PI/180)
170 greg 1.5 fprintf(stderr,
171     "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
172     progname, (s_longitude-s_meridian)*12/PI);
173 greg 1.1
174     printhead(argc, argv);
175    
176     computesky();
177     printsky();
178 greg 2.15
179     exit(0);
180 greg 1.1 }
181    
182    
183     computesky() /* compute sky parameters */
184     {
185 greg 2.12 double normfactor;
186 greg 1.1 /* compute solar direction */
187 greg 2.3 if (month) { /* from date and time */
188     int jd;
189     double sd, st;
190    
191     jd = jdate(month, day); /* Julian date */
192     sd = sdec(jd); /* solar declination */
193 greg 2.5 if (tsolar) /* solar time */
194     st = hour;
195     else
196     st = hour + stadj(jd);
197 greg 2.3 altitude = salt(sd, st);
198     azimuth = sazi(sd, st);
199 greg 2.17 printf("# Local solar time: %.2f\n", st);
200 greg 2.13 printf("# Solar altitude and azimuth: %.1f %.1f\n",
201 greg 2.11 180./PI*altitude, 180./PI*azimuth);
202 greg 2.9 }
203 greg 2.13 if (!overcast && altitude > 87.*PI/180.) {
204 greg 2.9 fprintf(stderr,
205     "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
206     progname);
207     printf(
208     "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
209     altitude = 87.*PI/180.;
210 greg 2.3 }
211 greg 1.1 sundir[0] = -sin(azimuth)*cos(altitude);
212     sundir[1] = -cos(azimuth)*cos(altitude);
213     sundir[2] = sin(altitude);
214    
215 greg 2.12 /* Compute normalization factor */
216 greg 2.13 switch (skytype) {
217     case S_UNIF:
218 greg 2.12 normfactor = 1.0;
219 greg 2.13 break;
220     case S_OVER:
221 greg 2.12 normfactor = 0.777778;
222 greg 2.13 break;
223     case S_CLEAR:
224 greg 2.12 F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
225     0.45*sundir[2]*sundir[2]);
226 greg 2.13 normfactor = normsc()/F2/PI;
227     break;
228     case S_INTER:
229     F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) *
230     exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
231     normfactor = normsc()/F2/PI;
232     break;
233 greg 2.12 }
234 greg 1.1 /* Compute zenith brightness */
235 greg 2.12 if (u_zenith == -1)
236     zenithbr /= normfactor*PI;
237     else if (u_zenith == 0) {
238 greg 2.13 if (overcast)
239 greg 1.1 zenithbr = 8.6*sundir[2] + .123;
240 greg 2.12 else
241 greg 1.1 zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
242 greg 2.13 if (skytype == S_INTER)
243     zenithbr = (zenithbr + 8.6*sundir[2] + .123)/2.0;
244 greg 2.12 if (zenithbr < 0.0)
245     zenithbr = 0.0;
246     else
247 greg 1.6 zenithbr *= 1000.0/SKYEFFICACY;
248 greg 2.12 }
249 greg 1.1 /* Compute horizontal radiance */
250 greg 2.12 groundbr = zenithbr*normfactor;
251 greg 2.13 printf("# Ground ambient level: %.1f\n", groundbr);
252 greg 2.14 if (!overcast && sundir[2] > 0.0 && (!u_solar || solarbr > 0.0)) {
253 greg 2.12 if (u_solar == -1)
254     solarbr /= 6e-5*sundir[2];
255 greg 2.13 else if (u_solar == 0) {
256 greg 2.12 solarbr = 1.5e9/SUNEFFICACY *
257     (1.147 - .147/(sundir[2]>.16?sundir[2]:.16));
258 greg 2.13 if (skytype == S_INTER)
259     solarbr *= 0.15; /* fudge factor! */
260     }
261 greg 2.12 groundbr += 6e-5/PI*solarbr*sundir[2];
262     } else
263     dosun = 0;
264 greg 1.1 groundbr *= gprefl;
265     }
266    
267    
268     printsky() /* print out sky */
269     {
270     if (dosun) {
271     printf("\nvoid light solar\n");
272     printf("0\n0\n");
273 greg 1.6 printf("3 %.2e %.2e %.2e\n", solarbr, solarbr, solarbr);
274 greg 1.1 printf("\nsolar source sun\n");
275     printf("0\n0\n");
276     printf("4 %f %f %f 0.5\n", sundir[0], sundir[1], sundir[2]);
277     }
278    
279     printf("\nvoid brightfunc skyfunc\n");
280 greg 2.13 printf("2 skybr skybright.cal\n");
281 greg 1.1 printf("0\n");
282 greg 2.13 if (overcast)
283     printf("3 %d %.2e %.2e\n", skytype, zenithbr, groundbr);
284 greg 1.1 else
285 greg 2.13 printf("7 %d %.2e %.2e %.2e %f %f %f\n",
286     skytype, zenithbr, groundbr, F2,
287     sundir[0], sundir[1], sundir[2]);
288 greg 1.1 }
289    
290    
291     printdefaults() /* print default values */
292     {
293 greg 2.13 switch (skytype) {
294     case S_OVER:
295 greg 1.1 printf("-c\t\t\t\t# Cloudy sky\n");
296 greg 2.13 break;
297     case S_UNIF:
298     printf("-u\t\t\t\t# Uniform cloudy sky\n");
299     break;
300     case S_INTER:
301     if (dosun)
302     printf("+i\t\t\t\t# Intermediate sky with sun\n");
303     else
304     printf("-i\t\t\t\t# Intermediate sky without sun\n");
305     break;
306     case S_CLEAR:
307     if (dosun)
308     printf("+s\t\t\t\t# Sunny sky with sun\n");
309     else
310     printf("-s\t\t\t\t# Sunny sky without sun\n");
311     break;
312     }
313 greg 1.1 printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
314     if (zenithbr > 0.0)
315     printf("-b %f\t\t\t# Zenith radiance (watts/ster/m2\n", zenithbr);
316     else
317     printf("-t %f\t\t\t# Atmospheric turbidity\n", turbidity);
318     printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
319     printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
320     printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
321     }
322    
323    
324     userror(msg) /* print usage error and quit */
325     char *msg;
326     {
327     if (msg != NULL)
328     fprintf(stderr, "%s: Use error - %s\n", progname, msg);
329     fprintf(stderr, "Usage: %s month day hour [options]\n", progname);
330 greg 2.3 fprintf(stderr, " Or: %s -ang altitude azimuth [options]\n", progname);
331 greg 1.1 fprintf(stderr, " Or: %s -defaults\n", progname);
332     exit(1);
333     }
334    
335    
336     double
337 greg 2.13 normsc() /* compute normalization factor (E0*F2/L0) */
338 greg 1.1 {
339 greg 2.13 static double nfc[2][5] = {
340     /* clear sky approx. */
341     {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
342     /* intermediate sky approx. */
343     {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
344     };
345     register double *nf;
346 greg 1.1 double x, nsc;
347     register int i;
348     /* polynomial approximation */
349 greg 2.13 nf = nfc[skytype==S_INTER];
350     x = (altitude - PI/4.0)/(PI/4.0);
351     nsc = nf[i=4];
352     while (i--)
353 greg 1.1 nsc = nsc*x + nf[i];
354    
355     return(nsc);
356 greg 2.16 }
357    
358    
359     cvthour(hs) /* convert hour string */
360     char *hs;
361     {
362     register char *cp = hs;
363 greg 2.17 register int i, j;
364 greg 2.16
365 greg 2.17 if (tsolar = *cp == '+') cp++; /* solar time? */
366     while (isdigit(*cp)) cp++;
367     if (*cp == ':')
368     hour = atoi(hs) + atoi(++cp)/60.0;
369     else {
370 greg 2.16 hour = atof(hs);
371 greg 2.17 if (*cp == '.') cp++;
372     }
373     while (isdigit(*cp)) cp++;
374     if (!*cp)
375     return;
376     if (tsolar || !isalpha(*cp)) {
377     fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
378     exit(1);
379     }
380     i = 0;
381     do {
382     for (j = 0; cp[j]; j++)
383     if (toupper(cp[j]) != tzone[i].zname[j])
384     break;
385     if (!cp[j] && !tzone[i].zname[j]) {
386     s_meridian = tzone[i].zmer * (PI/180);
387     return;
388     }
389     } while (tzone[i++].zname[0]);
390    
391     fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
392     fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
393     for (i = 1; tzone[i].zname[0]; i++)
394     fprintf(stderr, " %s", tzone[i].zname);
395     putc('\n', stderr);
396     exit(1);
397 greg 1.1 }
398    
399    
400     printhead(ac, av) /* print command header */
401     register int ac;
402     register char **av;
403     {
404     putchar('#');
405     while (ac--) {
406     putchar(' ');
407     fputs(*av++, stdout);
408     }
409     putchar('\n');
410     }