ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.5
Committed: Fri Mar 20 15:29:10 1992 UTC (32 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +7 -2 lines
Log Message:
added ability to specify solar instead of standard time

File Contents

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