ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.14
Committed: Mon Jan 31 12:53:49 1994 UTC (30 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.13: +1 -3 lines
Log Message:
fixed bug in calculation of ground brightness for uniform and overcast skies

File Contents

# Content
1 /* Copyright (c) 1992 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * gensky.c - program to generate sky functions.
9 * Our zenith is along the Z-axis, the X-axis
10 * points east, and the Y-axis points north.
11 * Radiance is in watts/steradian/sq. meter.
12 *
13 * 3/26/86
14 */
15
16 #include <stdio.h>
17
18 #include <math.h>
19
20 #include "color.h"
21
22 extern char *strcpy(), *strcat(), *malloc();
23 extern double stadj(), sdec(), sazi(), salt();
24
25 #define PI 3.141592654
26
27 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
28
29 #define S_CLEAR 1
30 #define S_OVER 2
31 #define S_UNIF 3
32 #define S_INTER 4
33
34 #define overcast (skytype==S_OVER|skytype==S_UNIF)
35
36 double normsc();
37 /* sun calculation constants */
38 extern double s_latitude;
39 extern double s_longitude;
40 extern double s_meridian;
41 /* required values */
42 int month, day; /* date */
43 double hour; /* time */
44 int tsolar; /* 0=standard, 1=solar */
45 double altitude, azimuth; /* or solar angles */
46 /* default values */
47 int skytype = S_CLEAR; /* sky type */
48 int dosun = 1;
49 double zenithbr = 0.0;
50 int u_zenith = 0; /* -1=irradiance, 1=radiance */
51 double turbidity = 2.75;
52 double gprefl = 0.2;
53 /* computed values */
54 double sundir[3];
55 double groundbr;
56 double F2;
57 double solarbr = 0.0;
58 int u_solar = 0; /* -1=irradiance, 1=radiance */
59
60 char *progname;
61 char errmsg[128];
62
63
64 main(argc, argv)
65 int argc;
66 char *argv[];
67 {
68 int i;
69
70 progname = argv[0];
71 if (argc == 2 && !strcmp(argv[1], "-defaults")) {
72 printdefaults();
73 exit(0);
74 }
75 if (argc < 4)
76 userror("arg count");
77 if (!strcmp(argv[1], "-ang")) {
78 altitude = atof(argv[2]) * (PI/180);
79 azimuth = atof(argv[3]) * (PI/180);
80 month = 0;
81 } else {
82 month = atoi(argv[1]);
83 if (month < 1 || month > 12)
84 userror("bad month");
85 day = atoi(argv[2]);
86 if (day < 1 || day > 31)
87 userror("bad day");
88 hour = atof(argv[3]);
89 if (hour < 0 || hour >= 24)
90 userror("bad hour");
91 tsolar = argv[3][0] == '+';
92 }
93 for (i = 4; i < argc; i++)
94 if (argv[i][0] == '-' || argv[i][0] == '+')
95 switch (argv[i][1]) {
96 case 's':
97 skytype = S_CLEAR;
98 dosun = argv[i][0] == '+';
99 break;
100 case 'r':
101 case 'R':
102 u_solar = argv[i][1]=='R' ? -1 : 1;
103 solarbr = atof(argv[++i]);
104 break;
105 case 'c':
106 skytype = S_OVER;
107 break;
108 case 'u':
109 skytype = S_UNIF;
110 break;
111 case 'i':
112 skytype = S_INTER;
113 dosun = argv[i][0] == '+';
114 break;
115 case 't':
116 turbidity = atof(argv[++i]);
117 break;
118 case 'b':
119 case 'B':
120 u_zenith = argv[i][1]=='B' ? -1 : 1;
121 zenithbr = atof(argv[++i]);
122 break;
123 case 'g':
124 gprefl = atof(argv[++i]);
125 break;
126 case 'a':
127 s_latitude = atof(argv[++i]) * (PI/180);
128 break;
129 case 'o':
130 s_longitude = atof(argv[++i]) * (PI/180);
131 break;
132 case 'm':
133 s_meridian = atof(argv[++i]) * (PI/180);
134 break;
135 default:
136 sprintf(errmsg, "unknown option: %s", argv[i]);
137 userror(errmsg);
138 }
139 else
140 userror("bad option");
141
142 if (fabs(s_meridian-s_longitude) > 30*PI/180)
143 fprintf(stderr,
144 "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
145 progname, (s_longitude-s_meridian)*12/PI);
146
147 printhead(argc, argv);
148
149 computesky();
150 printsky();
151 }
152
153
154 computesky() /* compute sky parameters */
155 {
156 double normfactor;
157 /* compute solar direction */
158 if (month) { /* from date and time */
159 int jd;
160 double sd, st;
161
162 jd = jdate(month, day); /* Julian date */
163 sd = sdec(jd); /* solar declination */
164 if (tsolar) /* solar time */
165 st = hour;
166 else
167 st = hour + stadj(jd);
168 altitude = salt(sd, st);
169 azimuth = sazi(sd, st);
170 printf("# Solar altitude and azimuth: %.1f %.1f\n",
171 180./PI*altitude, 180./PI*azimuth);
172 }
173 if (!overcast && altitude > 87.*PI/180.) {
174 fprintf(stderr,
175 "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
176 progname);
177 printf(
178 "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
179 altitude = 87.*PI/180.;
180 }
181 sundir[0] = -sin(azimuth)*cos(altitude);
182 sundir[1] = -cos(azimuth)*cos(altitude);
183 sundir[2] = sin(altitude);
184
185 /* Compute normalization factor */
186 switch (skytype) {
187 case S_UNIF:
188 normfactor = 1.0;
189 break;
190 case S_OVER:
191 normfactor = 0.777778;
192 break;
193 case S_CLEAR:
194 F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
195 0.45*sundir[2]*sundir[2]);
196 normfactor = normsc()/F2/PI;
197 break;
198 case S_INTER:
199 F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) *
200 exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
201 normfactor = normsc()/F2/PI;
202 break;
203 }
204 /* Compute zenith brightness */
205 if (u_zenith == -1)
206 zenithbr /= normfactor*PI;
207 else if (u_zenith == 0) {
208 if (overcast)
209 zenithbr = 8.6*sundir[2] + .123;
210 else
211 zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
212 if (skytype == S_INTER)
213 zenithbr = (zenithbr + 8.6*sundir[2] + .123)/2.0;
214 if (zenithbr < 0.0)
215 zenithbr = 0.0;
216 else
217 zenithbr *= 1000.0/SKYEFFICACY;
218 }
219 /* Compute horizontal radiance */
220 groundbr = zenithbr*normfactor;
221 printf("# Ground ambient level: %.1f\n", groundbr);
222 if (!overcast && sundir[2] > 0.0 && (!u_solar || solarbr > 0.0)) {
223 if (u_solar == -1)
224 solarbr /= 6e-5*sundir[2];
225 else if (u_solar == 0) {
226 solarbr = 1.5e9/SUNEFFICACY *
227 (1.147 - .147/(sundir[2]>.16?sundir[2]:.16));
228 if (skytype == S_INTER)
229 solarbr *= 0.15; /* fudge factor! */
230 }
231 groundbr += 6e-5/PI*solarbr*sundir[2];
232 } else
233 dosun = 0;
234 groundbr *= gprefl;
235 }
236
237
238 printsky() /* print out sky */
239 {
240 if (dosun) {
241 printf("\nvoid light solar\n");
242 printf("0\n0\n");
243 printf("3 %.2e %.2e %.2e\n", solarbr, solarbr, solarbr);
244 printf("\nsolar source sun\n");
245 printf("0\n0\n");
246 printf("4 %f %f %f 0.5\n", sundir[0], sundir[1], sundir[2]);
247 }
248
249 printf("\nvoid brightfunc skyfunc\n");
250 printf("2 skybr skybright.cal\n");
251 printf("0\n");
252 if (overcast)
253 printf("3 %d %.2e %.2e\n", skytype, zenithbr, groundbr);
254 else
255 printf("7 %d %.2e %.2e %.2e %f %f %f\n",
256 skytype, zenithbr, groundbr, F2,
257 sundir[0], sundir[1], sundir[2]);
258 }
259
260
261 printdefaults() /* print default values */
262 {
263 switch (skytype) {
264 case S_OVER:
265 printf("-c\t\t\t\t# Cloudy sky\n");
266 break;
267 case S_UNIF:
268 printf("-u\t\t\t\t# Uniform cloudy sky\n");
269 break;
270 case S_INTER:
271 if (dosun)
272 printf("+i\t\t\t\t# Intermediate sky with sun\n");
273 else
274 printf("-i\t\t\t\t# Intermediate sky without sun\n");
275 break;
276 case S_CLEAR:
277 if (dosun)
278 printf("+s\t\t\t\t# Sunny sky with sun\n");
279 else
280 printf("-s\t\t\t\t# Sunny sky without sun\n");
281 break;
282 }
283 printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
284 if (zenithbr > 0.0)
285 printf("-b %f\t\t\t# Zenith radiance (watts/ster/m2\n", zenithbr);
286 else
287 printf("-t %f\t\t\t# Atmospheric turbidity\n", turbidity);
288 printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
289 printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
290 printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
291 }
292
293
294 userror(msg) /* print usage error and quit */
295 char *msg;
296 {
297 if (msg != NULL)
298 fprintf(stderr, "%s: Use error - %s\n", progname, msg);
299 fprintf(stderr, "Usage: %s month day hour [options]\n", progname);
300 fprintf(stderr, " Or: %s -ang altitude azimuth [options]\n", progname);
301 fprintf(stderr, " Or: %s -defaults\n", progname);
302 exit(1);
303 }
304
305
306 double
307 normsc() /* compute normalization factor (E0*F2/L0) */
308 {
309 static double nfc[2][5] = {
310 /* clear sky approx. */
311 {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
312 /* intermediate sky approx. */
313 {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
314 };
315 register double *nf;
316 double x, nsc;
317 register int i;
318 /* polynomial approximation */
319 nf = nfc[skytype==S_INTER];
320 x = (altitude - PI/4.0)/(PI/4.0);
321 nsc = nf[i=4];
322 while (i--)
323 nsc = nsc*x + nf[i];
324
325 return(nsc);
326 }
327
328
329 printhead(ac, av) /* print command header */
330 register int ac;
331 register char **av;
332 {
333 putchar('#');
334 while (ac--) {
335 putchar(' ');
336 fputs(*av++, stdout);
337 }
338 putchar('\n');
339 }