ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
Revision: 1.5
Committed: Tue Aug 6 08:22:06 1991 UTC (32 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.4: +6 -1 lines
Log Message:
added warning about mismatched standard meridians and longitudes

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 greg 1.5 extern double atof(), fabs();
57 greg 1.1 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 greg 1.5
105     if (fabs(s_meridian-s_longitude) > 30*PI/180)
106     fprintf(stderr,
107     "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
108     progname, (s_longitude-s_meridian)*12/PI);
109 greg 1.1
110     printhead(argc, argv);
111    
112     computesky();
113     printsky();
114     }
115    
116    
117     computesky() /* compute sky parameters */
118     {
119     int jd;
120     double sd, st;
121     double altitude, azimuth;
122     /* compute solar direction */
123     jd = jdate(month, day); /* Julian date */
124     sd = sdec(jd); /* solar declination */
125     st = hour + stadj(jd); /* solar time */
126     altitude = salt(sd, st);
127     azimuth = sazi(sd, st);
128     sundir[0] = -sin(azimuth)*cos(altitude);
129     sundir[1] = -cos(azimuth)*cos(altitude);
130     sundir[2] = sin(altitude);
131    
132     /* Compute zenith brightness */
133     if (zenithbr <= 0.0)
134     if (cloudy) {
135     zenithbr = 8.6*sundir[2] + .123;
136 greg 1.4 zenithbr *= 1000.0*.0064/3.;
137 greg 1.1 } else {
138     zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
139 greg 1.4 zenithbr *= 1000.0*.0064/3.;
140 greg 1.1 }
141 greg 1.2 if (zenithbr < 0.0)
142     zenithbr = 0.0;
143 greg 1.1 /* Compute horizontal radiance */
144     if (cloudy) {
145     groundbr = zenithbr*0.777778;
146     printf("# Ground ambient level: %f\n", groundbr);
147     } else {
148     F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
149     0.45*sundir[2]*sundir[2]);
150     groundbr = zenithbr*normsc(PI/2.0-altitude)/F2/PI;
151     printf("# Ground ambient level: %f\n", groundbr);
152     if (sundir[2] > 0.0) {
153     if (sundir[2] > .16)
154     solarbr = 2.47e6 - 3.15e5/sundir[2];
155     else
156     solarbr = 5e5;
157     groundbr += solarbr*6e-5*sundir[2]/PI;
158     } else
159     dosun = 0;
160     }
161     groundbr *= gprefl;
162     }
163    
164    
165     printsky() /* print out sky */
166     {
167     if (dosun) {
168     printf("\nvoid light solar\n");
169     printf("0\n0\n");
170     printf("3 %f %f %f\n", solarbr, solarbr, solarbr);
171     printf("\nsolar source sun\n");
172     printf("0\n0\n");
173     printf("4 %f %f %f 0.5\n", sundir[0], sundir[1], sundir[2]);
174     }
175    
176     printf("\nvoid brightfunc skyfunc\n");
177     printf("2 skybright skybright.cal\n");
178     printf("0\n");
179     if (cloudy)
180     printf("3 1 %f %f\n", zenithbr, groundbr);
181     else
182     printf("7 -1 %f %f %f %f %f %f\n", zenithbr, groundbr, F2,
183     sundir[0], sundir[1], sundir[2]);
184     }
185    
186    
187     printdefaults() /* print default values */
188     {
189     if (cloudy)
190     printf("-c\t\t\t\t# Cloudy sky\n");
191     else if (dosun)
192     printf("+s\t\t\t\t# Sunny sky with sun\n");
193     else
194     printf("-s\t\t\t\t# Sunny sky without sun\n");
195     printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
196     if (zenithbr > 0.0)
197     printf("-b %f\t\t\t# Zenith radiance (watts/ster/m2\n", zenithbr);
198     else
199     printf("-t %f\t\t\t# Atmospheric turbidity\n", turbidity);
200     printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
201     printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
202     printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
203     }
204    
205    
206     userror(msg) /* print usage error and quit */
207     char *msg;
208     {
209     if (msg != NULL)
210     fprintf(stderr, "%s: Use error - %s\n", progname, msg);
211     fprintf(stderr, "Usage: %s month day hour [options]\n", progname);
212     fprintf(stderr, " Or: %s -defaults\n", progname);
213     exit(1);
214     }
215    
216    
217     double
218     normsc(theta) /* compute normalization factor (E0*F2/L0) */
219     double theta;
220     {
221     static double nf[5] = {2.766521, 0.547665,
222     -0.369832, 0.009237, 0.059229};
223     double x, nsc;
224     register int i;
225     /* polynomial approximation */
226     x = (theta - PI/4.0)/(PI/4.0);
227     nsc = nf[4];
228     for (i = 3; i >= 0; i--)
229     nsc = nsc*x + nf[i];
230    
231     return(nsc);
232     }
233    
234    
235     printhead(ac, av) /* print command header */
236     register int ac;
237     register char **av;
238     {
239     putchar('#');
240     while (ac--) {
241     putchar(' ');
242     fputs(*av++, stdout);
243     }
244     putchar('\n');
245     }