ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.6
Committed: Tue Apr 14 10:29:24 1992 UTC (31 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +8 -2 lines
Log Message:
fixed bug in sky normalization calculation,
which resulted in incorrect ground ambient levels

File Contents

# Content
1 /* Copyright (c) 1986 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 #ifndef atof
23 extern double atof();
24 #endif
25 extern char *strcpy(), *strcat(), *malloc();
26 extern double stadj(), sdec(), sazi(), salt();
27
28 #define PI 3.141592654
29
30 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
31
32 double normsc();
33 /* sun calculation constants */
34 extern double s_latitude;
35 extern double s_longitude;
36 extern double s_meridian;
37 /* required values */
38 int month, day; /* date */
39 double hour; /* time */
40 int tsolar; /* 0=standard, 1=solar */
41 double altitude, azimuth; /* or solar angles */
42 /* default values */
43 int cloudy = 0; /* 1=standard, 2=uniform */
44 int dosun = 1;
45 double zenithbr = -1.0;
46 double turbidity = 2.75;
47 double gprefl = 0.2;
48 /* computed values */
49 double sundir[3];
50 double groundbr;
51 double F2;
52 double solarbr;
53
54 char *progname;
55 char errmsg[128];
56
57
58 main(argc, argv)
59 int argc;
60 char *argv[];
61 {
62 extern double fabs();
63 int i;
64
65 progname = argv[0];
66 if (argc == 2 && !strcmp(argv[1], "-defaults")) {
67 printdefaults();
68 exit(0);
69 }
70 if (argc < 4)
71 userror("arg count");
72 if (!strcmp(argv[1], "-ang")) {
73 altitude = atof(argv[2]) * (PI/180);
74 azimuth = atof(argv[3]) * (PI/180);
75 month = 0;
76 } else {
77 month = atoi(argv[1]);
78 if (month < 1 || month > 12)
79 userror("bad month");
80 day = atoi(argv[2]);
81 if (day < 1 || day > 31)
82 userror("bad day");
83 hour = atof(argv[3]);
84 if (hour < 0 || hour >= 24)
85 userror("bad hour");
86 tsolar = argv[3][0] == '+';
87 }
88 for (i = 4; i < argc; i++)
89 if (argv[i][0] == '-' || argv[i][0] == '+')
90 switch (argv[i][1]) {
91 case 's':
92 cloudy = 0;
93 dosun = argv[i][0] == '+';
94 break;
95 case 'c':
96 cloudy = argv[i][0] == '+' ? 2 : 1;
97 dosun = 0;
98 break;
99 case 't':
100 turbidity = atof(argv[++i]);
101 break;
102 case 'b':
103 zenithbr = atof(argv[++i]);
104 break;
105 case 'g':
106 gprefl = atof(argv[++i]);
107 break;
108 case 'a':
109 s_latitude = atof(argv[++i]) * (PI/180);
110 break;
111 case 'o':
112 s_longitude = atof(argv[++i]) * (PI/180);
113 break;
114 case 'm':
115 s_meridian = atof(argv[++i]) * (PI/180);
116 break;
117 default:
118 sprintf(errmsg, "unknown option: %s", argv[i]);
119 userror(errmsg);
120 }
121 else
122 userror("bad option");
123
124 if (fabs(s_meridian-s_longitude) > 30*PI/180)
125 fprintf(stderr,
126 "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
127 progname, (s_longitude-s_meridian)*12/PI);
128
129 printhead(argc, argv);
130
131 computesky();
132 printsky();
133 }
134
135
136 computesky() /* compute sky parameters */
137 {
138 /* compute solar direction */
139 if (month) { /* from date and time */
140 int jd;
141 double sd, st;
142
143 jd = jdate(month, day); /* Julian date */
144 sd = sdec(jd); /* solar declination */
145 if (tsolar) /* solar time */
146 st = hour;
147 else
148 st = hour + stadj(jd);
149 altitude = salt(sd, st);
150 azimuth = sazi(sd, st);
151 }
152 sundir[0] = -sin(azimuth)*cos(altitude);
153 sundir[1] = -cos(azimuth)*cos(altitude);
154 sundir[2] = sin(altitude);
155
156 /* Compute zenith brightness */
157 if (zenithbr <= 0.0)
158 if (cloudy) {
159 zenithbr = 8.6*sundir[2] + .123;
160 zenithbr *= 1000.0/WHTEFFICACY;
161 } else {
162 zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
163 zenithbr *= 1000.0/SKYEFFICACY;
164 }
165 if (zenithbr < 0.0)
166 zenithbr = 0.0;
167 /* Compute horizontal radiance */
168 if (cloudy) {
169 if (cloudy == 2)
170 groundbr = zenithbr;
171 else
172 groundbr = zenithbr*0.777778;
173 printf("# Ground ambient level: %f\n", groundbr);
174 } else {
175 F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
176 0.45*sundir[2]*sundir[2]);
177 groundbr = zenithbr*normsc(altitude)/F2/PI;
178 printf("# Ground ambient level: %f\n", groundbr);
179 if (sundir[2] > 0.0) {
180 if (sundir[2] > .16)
181 solarbr = (1.5e9/SUNEFFICACY) *
182 (1.147 - .147/sundir[2]);
183 else
184 solarbr = 1.5e9/SUNEFFICACY*(1.147-.147/.16);
185 groundbr += solarbr*6e-5*sundir[2]/PI;
186 } else
187 dosun = 0;
188 }
189 groundbr *= gprefl;
190 }
191
192
193 printsky() /* print out sky */
194 {
195 if (dosun) {
196 printf("\nvoid light solar\n");
197 printf("0\n0\n");
198 printf("3 %.2e %.2e %.2e\n", solarbr, solarbr, solarbr);
199 printf("\nsolar source sun\n");
200 printf("0\n0\n");
201 printf("4 %f %f %f 0.5\n", sundir[0], sundir[1], sundir[2]);
202 }
203
204 printf("\nvoid brightfunc skyfunc\n");
205 printf("2 skybright skybright.cal\n");
206 printf("0\n");
207 if (cloudy)
208 printf("3 %d %.2e %.2e\n", cloudy, zenithbr, groundbr);
209 else
210 printf("7 -1 %.2e %.2e %.2e %f %f %f\n", zenithbr, groundbr,
211 F2, sundir[0], sundir[1], sundir[2]);
212 }
213
214
215 printdefaults() /* print default values */
216 {
217 if (cloudy == 1)
218 printf("-c\t\t\t\t# Cloudy sky\n");
219 else if (cloudy == 2)
220 printf("+c\t\t\t\t# Uniform cloudy sky\n");
221 else if (dosun)
222 printf("+s\t\t\t\t# Sunny sky with sun\n");
223 else
224 printf("-s\t\t\t\t# Sunny sky without sun\n");
225 printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
226 if (zenithbr > 0.0)
227 printf("-b %f\t\t\t# Zenith radiance (watts/ster/m2\n", zenithbr);
228 else
229 printf("-t %f\t\t\t# Atmospheric turbidity\n", turbidity);
230 printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
231 printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
232 printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
233 }
234
235
236 userror(msg) /* print usage error and quit */
237 char *msg;
238 {
239 if (msg != NULL)
240 fprintf(stderr, "%s: Use error - %s\n", progname, msg);
241 fprintf(stderr, "Usage: %s month day hour [options]\n", progname);
242 fprintf(stderr, " Or: %s -ang altitude azimuth [options]\n", progname);
243 fprintf(stderr, " Or: %s -defaults\n", progname);
244 exit(1);
245 }
246
247
248 double
249 normsc(theta) /* compute normalization factor (E0*F2/L0) */
250 double theta;
251 {
252 static double nf[5] = {2.766521, 0.547665,
253 -0.369832, 0.009237, 0.059229};
254 double x, nsc;
255 register int i;
256 /* polynomial approximation */
257 x = (theta - PI/4.0)/(PI/4.0);
258 nsc = nf[4];
259 for (i = 3; i >= 0; i--)
260 nsc = nsc*x + nf[i];
261
262 return(nsc);
263 }
264
265
266 printhead(ac, av) /* print command header */
267 register int ac;
268 register char **av;
269 {
270 putchar('#');
271 while (ac--) {
272 putchar(' ');
273 fputs(*av++, stdout);
274 }
275 putchar('\n');
276 }