ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.24
Committed: Wed Jul 26 17:15:11 2006 UTC (16 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0, rad3R8, rad3R9
Changes since 2.23: +2 -2 lines
Log Message:
Changed default turbidity from 2.75 to 2.45

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: gensky.c,v 2.23 2004/09/10 18:19:24 greg Exp $";
3 #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 #include <stdlib.h>
15 #include <string.h>
16 #include <math.h>
17 #include <ctype.h>
18
19 #include "color.h"
20
21 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
27 #ifndef PI
28 #define PI 3.14159265358979323846
29 #endif
30
31 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
32
33 #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 double normsc();
41 /* sun calculation constants */
42 extern double s_latitude;
43 extern double s_longitude;
44 extern double s_meridian;
45
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 {"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 };
71 /* required values */
72 int month, day; /* date */
73 double hour; /* time */
74 int tsolar; /* 0=standard, 1=solar */
75 double altitude, azimuth; /* or solar angles */
76 /* default values */
77 int skytype = S_CLEAR; /* sky type */
78 int dosun = 1;
79 double zenithbr = 0.0;
80 int u_zenith = 0; /* -1=irradiance, 1=radiance */
81 double turbidity = 2.45;
82 double gprefl = 0.2;
83 /* computed values */
84 double sundir[3];
85 double groundbr;
86 double F2;
87 double solarbr = 0.0;
88 int u_solar = 0; /* -1=irradiance, 1=radiance */
89
90 char *progname;
91 char errmsg[128];
92
93 void computesky(void);
94 void printsky(void);
95 void printdefaults(void);
96 void userror(char *msg);
97 double normsc(void);
98 int cvthour(char *hs);
99 void printhead(register int ac, register char **av);
100
101
102 int
103 main(argc, argv)
104 int argc;
105 char *argv[];
106 {
107 int got_meridian = 0;
108 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 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 if (month < 1 || month > 12)
124 userror("bad month");
125 day = atoi(argv[2]);
126 if (day < 1 || day > 31)
127 userror("bad day");
128 got_meridian = cvthour(argv[3]);
129 }
130 for (i = 4; i < argc; i++)
131 if (argv[i][0] == '-' || argv[i][0] == '+')
132 switch (argv[i][1]) {
133 case 's':
134 skytype = S_CLEAR;
135 dosun = argv[i][0] == '+';
136 break;
137 case 'r':
138 case 'R':
139 u_solar = argv[i][1]=='R' ? -1 : 1;
140 solarbr = atof(argv[++i]);
141 break;
142 case 'c':
143 skytype = S_OVER;
144 break;
145 case 'u':
146 skytype = S_UNIF;
147 break;
148 case 'i':
149 skytype = S_INTER;
150 dosun = argv[i][0] == '+';
151 break;
152 case 't':
153 turbidity = atof(argv[++i]);
154 break;
155 case 'b':
156 case 'B':
157 u_zenith = argv[i][1]=='B' ? -1 : 1;
158 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 if (got_meridian) {
171 ++i;
172 break; /* time overrides */
173 }
174 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
183 if (fabs(s_meridian-s_longitude) > 45*PI/180)
184 fprintf(stderr,
185 "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
186 progname, (s_longitude-s_meridian)*12/PI);
187
188 printhead(argc, argv);
189
190 computesky();
191 printsky();
192
193 exit(0);
194 }
195
196
197 void
198 computesky(void) /* compute sky parameters */
199 {
200 double normfactor;
201 /* compute solar direction */
202 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 if (tsolar) /* solar time */
209 st = hour;
210 else
211 st = hour + stadj(jd);
212 altitude = salt(sd, st);
213 azimuth = sazi(sd, st);
214 printf("# Local solar time: %.2f\n", st);
215 printf("# Solar altitude and azimuth: %.1f %.1f\n",
216 180./PI*altitude, 180./PI*azimuth);
217 }
218 if (!overcast && altitude > 87.*PI/180.) {
219 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 }
226 sundir[0] = -sin(azimuth)*cos(altitude);
227 sundir[1] = -cos(azimuth)*cos(altitude);
228 sundir[2] = sin(altitude);
229
230 /* Compute normalization factor */
231 switch (skytype) {
232 case S_UNIF:
233 normfactor = 1.0;
234 break;
235 case S_OVER:
236 normfactor = 0.777778;
237 break;
238 case S_CLEAR:
239 F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
240 0.45*sundir[2]*sundir[2]);
241 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 }
249 /* Compute zenith brightness */
250 if (u_zenith == -1)
251 zenithbr /= normfactor*PI;
252 else if (u_zenith == 0) {
253 if (overcast)
254 zenithbr = 8.6*sundir[2] + .123;
255 else
256 zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
257 if (skytype == S_INTER)
258 zenithbr = (zenithbr + 8.6*sundir[2] + .123)/2.0;
259 if (zenithbr < 0.0)
260 zenithbr = 0.0;
261 else
262 zenithbr *= 1000.0/SKYEFFICACY;
263 }
264 /* Compute horizontal radiance */
265 groundbr = zenithbr*normfactor;
266 printf("# Ground ambient level: %.1f\n", groundbr);
267 if (!overcast && sundir[2] > 0.0 && (!u_solar || solarbr > 0.0)) {
268 if (u_solar == -1)
269 solarbr /= 6e-5*sundir[2];
270 else if (u_solar == 0) {
271 solarbr = 1.5e9/SUNEFFICACY *
272 (1.147 - .147/(sundir[2]>.16?sundir[2]:.16));
273 if (skytype == S_INTER)
274 solarbr *= 0.15; /* fudge factor! */
275 }
276 groundbr += 6e-5/PI*solarbr*sundir[2];
277 } else
278 dosun = 0;
279 groundbr *= gprefl;
280 }
281
282
283 void
284 printsky(void) /* print out sky */
285 {
286 if (dosun) {
287 printf("\nvoid light solar\n");
288 printf("0\n0\n");
289 printf("3 %.2e %.2e %.2e\n", solarbr, solarbr, solarbr);
290 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 printf("2 skybr skybright.cal\n");
297 printf("0\n");
298 if (overcast)
299 printf("3 %d %.2e %.2e\n", skytype, zenithbr, groundbr);
300 else
301 printf("7 %d %.2e %.2e %.2e %f %f %f\n",
302 skytype, zenithbr, groundbr, F2,
303 sundir[0], sundir[1], sundir[2]);
304 }
305
306
307 void
308 printdefaults(void) /* print default values */
309 {
310 switch (skytype) {
311 case S_OVER:
312 printf("-c\t\t\t\t# Cloudy sky\n");
313 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 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 void
342 userror( /* print usage error and quit */
343 char *msg
344 )
345 {
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 fprintf(stderr, " Or: %s -ang altitude azimuth [options]\n", progname);
350 fprintf(stderr, " Or: %s -defaults\n", progname);
351 exit(1);
352 }
353
354
355 double
356 normsc(void) /* compute normalization factor (E0*F2/L0) */
357 {
358 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 double x, nsc;
366 register int i;
367 /* polynomial approximation */
368 nf = nfc[skytype==S_INTER];
369 x = (altitude - PI/4.0)/(PI/4.0);
370 nsc = nf[i=4];
371 while (i--)
372 nsc = nsc*x + nf[i];
373
374 return(nsc);
375 }
376
377
378 int
379 cvthour( /* convert hour string */
380 char *hs
381 )
382 {
383 register char *cp = hs;
384 register int i, j;
385
386 if ( (tsolar = *cp == '+') ) cp++; /* solar time? */
387 while (isdigit(*cp)) cp++;
388 if (*cp == ':')
389 hour = atoi(hs) + atoi(++cp)/60.0;
390 else {
391 hour = atof(hs);
392 if (*cp == '.') cp++;
393 }
394 while (isdigit(*cp)) cp++;
395 if (!*cp)
396 return(0);
397 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 return(1);
409 }
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 }
419
420
421 void
422 printhead( /* print command header */
423 register int ac,
424 register char **av
425 )
426 {
427 putchar('#');
428 while (ac--) {
429 putchar(' ');
430 fputs(*av++, stdout);
431 }
432 putchar('\n');
433 }