ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 1.4
Committed: Tue Mar 19 12:09:27 1991 UTC (33 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.3: +2 -2 lines
Log Message:
changed calculation of zenith brightness to use white light

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