ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.23
Committed: Fri Sep 10 18:19:24 2004 UTC (19 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad3R6, rad3R6P1
Changes since 2.22: +11 -6 lines
Log Message:
Made it so -m option is ignored in gensky if time zone is given

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.23 static const char RCSid[] = "$Id: gensky.c,v 2.22 2003/11/16 10:29:38 schorsch Exp $";
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 greg 2.20 #include <stdlib.h>
15     #include <string.h>
16 greg 1.1 #include <math.h>
17 greg 2.17 #include <ctype.h>
18    
19 greg 1.6 #include "color.h"
20 greg 1.1
21 schorsch 2.22 extern int jdate(int month, int day);
22     extern double stadj(int jd);
23     extern double sdec(int jd);
24     extern double salt(double sd, double st);
25     extern double sazi(double sd, double st);
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 schorsch 2.21 #define overcast ((skytype==S_OVER)|(skytype==S_UNIF))
39 greg 2.13
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 schorsch 2.22 {"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     {"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     {"", 0}
70 greg 2.17 };
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 schorsch 2.22 void computesky(void);
94     void printsky(void);
95     void printdefaults(void);
96     void userror(char *msg);
97     double normsc(void);
98 greg 2.23 int cvthour(char *hs);
99 schorsch 2.22 void printhead(register int ac, register char **av);
100    
101 greg 1.1
102 schorsch 2.22 int
103 greg 1.1 main(argc, argv)
104     int argc;
105     char *argv[];
106     {
107 greg 2.23 int got_meridian = 0;
108 greg 1.1 int i;
109    
110     progname = argv[0];
111     if (argc == 2 && !strcmp(argv[1], "-defaults")) {
112     printdefaults();
113     exit(0);
114     }
115     if (argc < 4)
116     userror("arg count");
117 greg 2.3 if (!strcmp(argv[1], "-ang")) {
118     altitude = atof(argv[2]) * (PI/180);
119     azimuth = atof(argv[3]) * (PI/180);
120     month = 0;
121     } else {
122     month = atoi(argv[1]);
123 greg 2.6 if (month < 1 || month > 12)
124     userror("bad month");
125 greg 2.3 day = atoi(argv[2]);
126 greg 2.6 if (day < 1 || day > 31)
127     userror("bad day");
128 greg 2.23 got_meridian = cvthour(argv[3]);
129 greg 2.3 }
130 greg 1.1 for (i = 4; i < argc; i++)
131     if (argv[i][0] == '-' || argv[i][0] == '+')
132     switch (argv[i][1]) {
133     case 's':
134 greg 2.13 skytype = S_CLEAR;
135 greg 1.1 dosun = argv[i][0] == '+';
136     break;
137 greg 2.8 case 'r':
138 greg 2.12 case 'R':
139     u_solar = argv[i][1]=='R' ? -1 : 1;
140 greg 2.8 solarbr = atof(argv[++i]);
141     break;
142 greg 1.1 case 'c':
143 greg 2.13 skytype = S_OVER;
144 greg 1.1 break;
145 greg 2.13 case 'u':
146     skytype = S_UNIF;
147     break;
148     case 'i':
149     skytype = S_INTER;
150     dosun = argv[i][0] == '+';
151     break;
152 greg 1.1 case 't':
153     turbidity = atof(argv[++i]);
154     break;
155     case 'b':
156 greg 2.12 case 'B':
157     u_zenith = argv[i][1]=='B' ? -1 : 1;
158 greg 1.1 zenithbr = atof(argv[++i]);
159     break;
160     case 'g':
161     gprefl = atof(argv[++i]);
162     break;
163     case 'a':
164     s_latitude = atof(argv[++i]) * (PI/180);
165     break;
166     case 'o':
167     s_longitude = atof(argv[++i]) * (PI/180);
168     break;
169     case 'm':
170 greg 2.23 if (got_meridian) {
171     ++i;
172     break; /* time overrides */
173     }
174 greg 1.1 s_meridian = atof(argv[++i]) * (PI/180);
175     break;
176     default:
177     sprintf(errmsg, "unknown option: %s", argv[i]);
178     userror(errmsg);
179     }
180     else
181     userror("bad option");
182 greg 1.5
183 greg 2.18 if (fabs(s_meridian-s_longitude) > 45*PI/180)
184 greg 1.5 fprintf(stderr,
185     "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
186     progname, (s_longitude-s_meridian)*12/PI);
187 greg 1.1
188     printhead(argc, argv);
189    
190     computesky();
191     printsky();
192 greg 2.15
193     exit(0);
194 greg 1.1 }
195    
196    
197 schorsch 2.22 void
198     computesky(void) /* compute sky parameters */
199 greg 1.1 {
200 greg 2.12 double normfactor;
201 greg 1.1 /* compute solar direction */
202 greg 2.3 if (month) { /* from date and time */
203     int jd;
204     double sd, st;
205    
206     jd = jdate(month, day); /* Julian date */
207     sd = sdec(jd); /* solar declination */
208 greg 2.5 if (tsolar) /* solar time */
209     st = hour;
210     else
211     st = hour + stadj(jd);
212 greg 2.3 altitude = salt(sd, st);
213     azimuth = sazi(sd, st);
214 greg 2.17 printf("# Local solar time: %.2f\n", st);
215 greg 2.13 printf("# Solar altitude and azimuth: %.1f %.1f\n",
216 greg 2.11 180./PI*altitude, 180./PI*azimuth);
217 greg 2.9 }
218 greg 2.13 if (!overcast && altitude > 87.*PI/180.) {
219 greg 2.9 fprintf(stderr,
220     "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
221     progname);
222     printf(
223     "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
224     altitude = 87.*PI/180.;
225 greg 2.3 }
226 greg 1.1 sundir[0] = -sin(azimuth)*cos(altitude);
227     sundir[1] = -cos(azimuth)*cos(altitude);
228     sundir[2] = sin(altitude);
229    
230 greg 2.12 /* Compute normalization factor */
231 greg 2.13 switch (skytype) {
232     case S_UNIF:
233 greg 2.12 normfactor = 1.0;
234 greg 2.13 break;
235     case S_OVER:
236 greg 2.12 normfactor = 0.777778;
237 greg 2.13 break;
238     case S_CLEAR:
239 greg 2.12 F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
240     0.45*sundir[2]*sundir[2]);
241 greg 2.13 normfactor = normsc()/F2/PI;
242     break;
243     case S_INTER:
244     F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) *
245     exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
246     normfactor = normsc()/F2/PI;
247     break;
248 greg 2.12 }
249 greg 1.1 /* Compute zenith brightness */
250 greg 2.12 if (u_zenith == -1)
251     zenithbr /= normfactor*PI;
252     else if (u_zenith == 0) {
253 greg 2.13 if (overcast)
254 greg 1.1 zenithbr = 8.6*sundir[2] + .123;
255 greg 2.12 else
256 greg 1.1 zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
257 greg 2.13 if (skytype == S_INTER)
258     zenithbr = (zenithbr + 8.6*sundir[2] + .123)/2.0;
259 greg 2.12 if (zenithbr < 0.0)
260     zenithbr = 0.0;
261     else
262 greg 1.6 zenithbr *= 1000.0/SKYEFFICACY;
263 greg 2.12 }
264 greg 1.1 /* Compute horizontal radiance */
265 greg 2.12 groundbr = zenithbr*normfactor;
266 greg 2.13 printf("# Ground ambient level: %.1f\n", groundbr);
267 greg 2.14 if (!overcast && sundir[2] > 0.0 && (!u_solar || solarbr > 0.0)) {
268 greg 2.12 if (u_solar == -1)
269     solarbr /= 6e-5*sundir[2];
270 greg 2.13 else if (u_solar == 0) {
271 greg 2.12 solarbr = 1.5e9/SUNEFFICACY *
272     (1.147 - .147/(sundir[2]>.16?sundir[2]:.16));
273 greg 2.13 if (skytype == S_INTER)
274     solarbr *= 0.15; /* fudge factor! */
275     }
276 greg 2.12 groundbr += 6e-5/PI*solarbr*sundir[2];
277     } else
278     dosun = 0;
279 greg 1.1 groundbr *= gprefl;
280     }
281    
282    
283 schorsch 2.22 void
284     printsky(void) /* print out sky */
285 greg 1.1 {
286     if (dosun) {
287     printf("\nvoid light solar\n");
288     printf("0\n0\n");
289 greg 1.6 printf("3 %.2e %.2e %.2e\n", solarbr, solarbr, solarbr);
290 greg 1.1 printf("\nsolar source sun\n");
291     printf("0\n0\n");
292     printf("4 %f %f %f 0.5\n", sundir[0], sundir[1], sundir[2]);
293     }
294    
295     printf("\nvoid brightfunc skyfunc\n");
296 greg 2.13 printf("2 skybr skybright.cal\n");
297 greg 1.1 printf("0\n");
298 greg 2.13 if (overcast)
299     printf("3 %d %.2e %.2e\n", skytype, zenithbr, groundbr);
300 greg 1.1 else
301 greg 2.13 printf("7 %d %.2e %.2e %.2e %f %f %f\n",
302     skytype, zenithbr, groundbr, F2,
303     sundir[0], sundir[1], sundir[2]);
304 greg 1.1 }
305    
306    
307 schorsch 2.22 void
308     printdefaults(void) /* print default values */
309 greg 1.1 {
310 greg 2.13 switch (skytype) {
311     case S_OVER:
312 greg 1.1 printf("-c\t\t\t\t# Cloudy sky\n");
313 greg 2.13 break;
314     case S_UNIF:
315     printf("-u\t\t\t\t# Uniform cloudy sky\n");
316     break;
317     case S_INTER:
318     if (dosun)
319     printf("+i\t\t\t\t# Intermediate sky with sun\n");
320     else
321     printf("-i\t\t\t\t# Intermediate sky without sun\n");
322     break;
323     case S_CLEAR:
324     if (dosun)
325     printf("+s\t\t\t\t# Sunny sky with sun\n");
326     else
327     printf("-s\t\t\t\t# Sunny sky without sun\n");
328     break;
329     }
330 greg 1.1 printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
331     if (zenithbr > 0.0)
332     printf("-b %f\t\t\t# Zenith radiance (watts/ster/m2\n", zenithbr);
333     else
334     printf("-t %f\t\t\t# Atmospheric turbidity\n", turbidity);
335     printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
336     printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
337     printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
338     }
339    
340    
341 schorsch 2.22 void
342     userror( /* print usage error and quit */
343     char *msg
344     )
345 greg 1.1 {
346     if (msg != NULL)
347     fprintf(stderr, "%s: Use error - %s\n", progname, msg);
348     fprintf(stderr, "Usage: %s month day hour [options]\n", progname);
349 greg 2.3 fprintf(stderr, " Or: %s -ang altitude azimuth [options]\n", progname);
350 greg 1.1 fprintf(stderr, " Or: %s -defaults\n", progname);
351     exit(1);
352     }
353    
354    
355     double
356 schorsch 2.22 normsc(void) /* compute normalization factor (E0*F2/L0) */
357 greg 1.1 {
358 greg 2.13 static double nfc[2][5] = {
359     /* clear sky approx. */
360     {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
361     /* intermediate sky approx. */
362     {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
363     };
364     register double *nf;
365 greg 1.1 double x, nsc;
366     register int i;
367     /* polynomial approximation */
368 greg 2.13 nf = nfc[skytype==S_INTER];
369     x = (altitude - PI/4.0)/(PI/4.0);
370     nsc = nf[i=4];
371     while (i--)
372 greg 1.1 nsc = nsc*x + nf[i];
373    
374     return(nsc);
375 greg 2.16 }
376    
377    
378 greg 2.23 int
379 schorsch 2.22 cvthour( /* convert hour string */
380     char *hs
381     )
382 greg 2.16 {
383     register char *cp = hs;
384 greg 2.17 register int i, j;
385 greg 2.16
386 schorsch 2.21 if ( (tsolar = *cp == '+') ) cp++; /* solar time? */
387 greg 2.17 while (isdigit(*cp)) cp++;
388     if (*cp == ':')
389     hour = atoi(hs) + atoi(++cp)/60.0;
390     else {
391 greg 2.16 hour = atof(hs);
392 greg 2.17 if (*cp == '.') cp++;
393     }
394     while (isdigit(*cp)) cp++;
395     if (!*cp)
396 greg 2.23 return(0);
397 greg 2.17 if (tsolar || !isalpha(*cp)) {
398     fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
399     exit(1);
400     }
401     i = 0;
402     do {
403     for (j = 0; cp[j]; j++)
404     if (toupper(cp[j]) != tzone[i].zname[j])
405     break;
406     if (!cp[j] && !tzone[i].zname[j]) {
407     s_meridian = tzone[i].zmer * (PI/180);
408 greg 2.23 return(1);
409 greg 2.17 }
410     } while (tzone[i++].zname[0]);
411    
412     fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
413     fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
414     for (i = 1; tzone[i].zname[0]; i++)
415     fprintf(stderr, " %s", tzone[i].zname);
416     putc('\n', stderr);
417     exit(1);
418 greg 1.1 }
419    
420    
421 schorsch 2.22 void
422     printhead( /* print command header */
423     register int ac,
424     register char **av
425     )
426 greg 1.1 {
427     putchar('#');
428     while (ac--) {
429     putchar(' ');
430     fputs(*av++, stdout);
431     }
432     putchar('\n');
433     }