ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 2.4
Committed: Fri Mar 20 12:23:38 1992 UTC (32 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +10 -5 lines
Log Message:
added +c option for uniform cloudy sky

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