ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.26
Committed: Wed Jul 30 17:30:27 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R2P1
Changes since 2.25: +13 -22 lines
Log Message:
Improved use of sun.h

File Contents

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