ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensky.c
(Generate patch)

Comparing ray/src/gen/gensky.c (file contents):
Revision 1.6 by greg, Thu Oct 24 13:43:22 1991 UTC vs.
Revision 2.22 by schorsch, Sun Nov 16 10:29:38 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1986 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   *  gensky.c - program to generate sky functions.
6   *              Our zenith is along the Z-axis, the X-axis
# Line 14 | Line 11 | static char SCCSid[] = "$SunId$ LBL";
11   */
12  
13   #include  <stdio.h>
14 <
14 > #include  <stdlib.h>
15 > #include  <string.h>
16   #include  <math.h>
17 + #include  <ctype.h>
18  
19   #include  "color.h"
20  
21 < extern char  *strcpy(), *strcat(), *malloc();
22 < extern double  stadj(), sdec(), sazi(), salt();
21 > extern int jdate(int month, int day);
22 > extern double stadj(int  jd);
23 > extern double sdec(int  jd);
24 > extern double salt(double sd, double st);
25 > extern double sazi(double sd, double st);
26  
27 < #define  PI             3.141592654
27 > #ifndef  PI
28 > #define  PI             3.14159265358979323846
29 > #endif
30  
31   #define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
32  
33 + #define  S_CLEAR        1
34 + #define  S_OVER         2
35 + #define  S_UNIF         3
36 + #define  S_INTER        4
37 +
38 + #define  overcast       ((skytype==S_OVER)|(skytype==S_UNIF))
39 +
40   double  normsc();
41                                          /* sun calculation constants */
42   extern double  s_latitude;
43   extern double  s_longitude;
44   extern double  s_meridian;
45 +
46 + #undef  toupper
47 + #define  toupper(c)     ((c) & ~0x20)   /* ASCII trick to convert case */
48 +
49 +                                        /* European and North American zones */
50 + struct {
51 +        char    zname[8];       /* time zone name (all caps) */
52 +        float   zmer;           /* standard meridian */
53 + } tzone[] = {
54 +        {"YST", 135}, {"YDT", 120},
55 +        {"PST", 120}, {"PDT", 105},
56 +        {"MST", 105}, {"MDT", 90},
57 +        {"CST", 90}, {"CDT", 75},
58 +        {"EST", 75}, {"EDT", 60},
59 +        {"AST", 60}, {"ADT", 45},
60 +        {"NST", 52.5}, {"NDT", 37.5},
61 +        {"GMT", 0}, {"BST", -15},
62 +        {"CET", -15}, {"CEST", -30},
63 +        {"EET", -30}, {"EEST", -45},
64 +        {"AST", -45}, {"ADT", -60},
65 +        {"GST", -60}, {"GDT", -75},
66 +        {"IST", -82.5}, {"IDT", -97.5},
67 +        {"JST", -135}, {"NDT", -150},
68 +        {"NZST", -180}, {"NZDT", -195},
69 +        {"", 0}
70 + };
71                                          /* required values */
72 < int  month, day;
73 < double  hour;
72 > int  month, day;                                /* date */
73 > double  hour;                                   /* time */
74 > int  tsolar;                                    /* 0=standard, 1=solar */
75 > double  altitude, azimuth;                      /* or solar angles */
76                                          /* default values */
77 < int  cloudy = 0;
77 > int  skytype = S_CLEAR;                         /* sky type */
78   int  dosun = 1;
79 < double  zenithbr = -1.0;
79 > double  zenithbr = 0.0;
80 > int     u_zenith = 0;                           /* -1=irradiance, 1=radiance */
81   double  turbidity = 2.75;
82   double  gprefl = 0.2;
83                                          /* computed values */
84   double  sundir[3];
85   double  groundbr;
86   double  F2;
87 < double  solarbr;
87 > double  solarbr = 0.0;
88 > int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
89  
90   char  *progname;
91   char  errmsg[128];
92  
93 + void computesky(void);
94 + void printsky(void);
95 + void printdefaults(void);
96 + void userror(char  *msg);
97 + double normsc(void);
98 + void cvthour(char  *hs);
99 + void printhead(register int  ac, register char  **av);
100  
101 +
102 + int
103   main(argc, argv)
104   int  argc;
105   char  *argv[];
106   {
57 –        extern double  atof(), fabs();
107          int  i;
108  
109          progname = argv[0];
# Line 64 | Line 113 | char  *argv[];
113          }
114          if (argc < 4)
115                  userror("arg count");
116 <        month = atoi(argv[1]);
117 <        day = atoi(argv[2]);
118 <        hour = atof(argv[3]);
116 >        if (!strcmp(argv[1], "-ang")) {
117 >                altitude = atof(argv[2]) * (PI/180);
118 >                azimuth = atof(argv[3]) * (PI/180);
119 >                month = 0;
120 >        } else {
121 >                month = atoi(argv[1]);
122 >                if (month < 1 || month > 12)
123 >                        userror("bad month");
124 >                day = atoi(argv[2]);
125 >                if (day < 1 || day > 31)
126 >                        userror("bad day");
127 >                cvthour(argv[3]);
128 >        }
129          for (i = 4; i < argc; i++)
130                  if (argv[i][0] == '-' || argv[i][0] == '+')
131                          switch (argv[i][1]) {
132                          case 's':
133 <                                cloudy = 0;
133 >                                skytype = S_CLEAR;
134                                  dosun = argv[i][0] == '+';
135                                  break;
136 +                        case 'r':
137 +                        case 'R':
138 +                                u_solar = argv[i][1]=='R' ? -1 : 1;
139 +                                solarbr = atof(argv[++i]);
140 +                                break;
141                          case 'c':
142 <                                cloudy = 1;
79 <                                dosun = 0;
142 >                                skytype = S_OVER;
143                                  break;
144 +                        case 'u':
145 +                                skytype = S_UNIF;
146 +                                break;
147 +                        case 'i':
148 +                                skytype = S_INTER;
149 +                                dosun = argv[i][0] == '+';
150 +                                break;
151                          case 't':
152                                  turbidity = atof(argv[++i]);
153                                  break;
154                          case 'b':
155 +                        case 'B':
156 +                                u_zenith = argv[i][1]=='B' ? -1 : 1;
157                                  zenithbr = atof(argv[++i]);
158                                  break;
159                          case 'g':
# Line 103 | Line 175 | char  *argv[];
175                  else
176                          userror("bad option");
177  
178 <        if (fabs(s_meridian-s_longitude) > 30*PI/180)
178 >        if (fabs(s_meridian-s_longitude) > 45*PI/180)
179                  fprintf(stderr,
180          "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
181                          progname, (s_longitude-s_meridian)*12/PI);
# Line 112 | Line 184 | char  *argv[];
184  
185          computesky();
186          printsky();
187 +
188 +        exit(0);
189   }
190  
191  
192 < computesky()                    /* compute sky parameters */
192 > void
193 > computesky(void)                        /* compute sky parameters */
194   {
195 <        int  jd;
121 <        double  sd, st;
122 <        double  altitude, azimuth;
195 >        double  normfactor;
196                                          /* compute solar direction */
197 <        jd = jdate(month, day);                 /* Julian date */
198 <        sd = sdec(jd);                          /* solar declination */
199 <        st = hour + stadj(jd);                  /* solar time */
200 <        altitude = salt(sd, st);
201 <        azimuth = sazi(sd, st);
197 >        if (month) {                    /* from date and time */
198 >                int  jd;
199 >                double  sd, st;
200 >
201 >                jd = jdate(month, day);         /* Julian date */
202 >                sd = sdec(jd);                  /* solar declination */
203 >                if (tsolar)                     /* solar time */
204 >                        st = hour;
205 >                else
206 >                        st = hour + stadj(jd);
207 >                altitude = salt(sd, st);
208 >                azimuth = sazi(sd, st);
209 >                printf("# Local solar time: %.2f\n", st);
210 >                printf("# Solar altitude and azimuth: %.1f %.1f\n",
211 >                                180./PI*altitude, 180./PI*azimuth);
212 >        }
213 >        if (!overcast && altitude > 87.*PI/180.) {
214 >                fprintf(stderr,
215 > "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
216 >                                progname);
217 >                printf(
218 > "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
219 >                altitude = 87.*PI/180.;
220 >        }
221          sundir[0] = -sin(azimuth)*cos(altitude);
222          sundir[1] = -cos(azimuth)*cos(altitude);
223          sundir[2] = sin(altitude);
224  
225 +                                        /* Compute normalization factor */
226 +        switch (skytype) {
227 +        case S_UNIF:
228 +                normfactor = 1.0;
229 +                break;
230 +        case S_OVER:
231 +                normfactor = 0.777778;
232 +                break;
233 +        case S_CLEAR:
234 +                F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
235 +                                0.45*sundir[2]*sundir[2]);
236 +                normfactor = normsc()/F2/PI;
237 +                break;
238 +        case S_INTER:
239 +                F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) *
240 +                        exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
241 +                normfactor = normsc()/F2/PI;
242 +                break;
243 +        }
244                                          /* Compute zenith brightness */
245 <        if (zenithbr <= 0.0)
246 <                if (cloudy) {
245 >        if (u_zenith == -1)
246 >                zenithbr /= normfactor*PI;
247 >        else if (u_zenith == 0) {
248 >                if (overcast)
249                          zenithbr = 8.6*sundir[2] + .123;
250 <                        zenithbr *= 1000.0/SKYEFFICACY;
138 <                } else {
250 >                else
251                          zenithbr = (1.376*turbidity-1.81)*tan(altitude)+0.38;
252 +                if (skytype == S_INTER)
253 +                        zenithbr = (zenithbr + 8.6*sundir[2] + .123)/2.0;
254 +                if (zenithbr < 0.0)
255 +                        zenithbr = 0.0;
256 +                else
257                          zenithbr *= 1000.0/SKYEFFICACY;
141 –                }
142 –        if (zenithbr < 0.0)
143 –                zenithbr = 0.0;
144 –                                        /* Compute horizontal radiance */
145 –        if (cloudy) {
146 –                groundbr = zenithbr*0.777778;
147 –                printf("# Ground ambient level: %f\n", groundbr);
148 –        } else {
149 –                F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) +
150 –                                0.45*sundir[2]*sundir[2]);
151 –                groundbr = zenithbr*normsc(PI/2.0-altitude)/F2/PI;
152 –                printf("# Ground ambient level: %f\n", groundbr);
153 –                if (sundir[2] > 0.0) {
154 –                        if (sundir[2] > .16)
155 –                                solarbr = (1.5e9/SUNEFFICACY) *
156 –                                        (1.147 - .147/sundir[2]);
157 –                        else
158 –                                solarbr = 1.5e9/SUNEFFICACY*(1.147-.147/.16);
159 –                        groundbr += solarbr*6e-5*sundir[2]/PI;
160 –                } else
161 –                        dosun = 0;
258          }
259 +                                        /* Compute horizontal radiance */
260 +        groundbr = zenithbr*normfactor;
261 +        printf("# Ground ambient level: %.1f\n", groundbr);
262 +        if (!overcast && sundir[2] > 0.0 && (!u_solar || solarbr > 0.0)) {
263 +                if (u_solar == -1)
264 +                        solarbr /= 6e-5*sundir[2];
265 +                else if (u_solar == 0) {
266 +                        solarbr = 1.5e9/SUNEFFICACY *
267 +                        (1.147 - .147/(sundir[2]>.16?sundir[2]:.16));
268 +                        if (skytype == S_INTER)
269 +                                solarbr *= 0.15;        /* fudge factor! */
270 +                }
271 +                groundbr += 6e-5/PI*solarbr*sundir[2];
272 +        } else
273 +                dosun = 0;
274          groundbr *= gprefl;
275   }
276  
277  
278 < printsky()                      /* print out sky */
278 > void
279 > printsky(void)                  /* print out sky */
280   {
281          if (dosun) {
282                  printf("\nvoid light solar\n");
# Line 176 | Line 288 | printsky()                     /* print out sky */
288          }
289          
290          printf("\nvoid brightfunc skyfunc\n");
291 <        printf("2 skybright skybright.cal\n");
291 >        printf("2 skybr skybright.cal\n");
292          printf("0\n");
293 <        if (cloudy)
294 <                printf("3 1 %.2e %.2e\n", zenithbr, groundbr);
293 >        if (overcast)
294 >                printf("3 %d %.2e %.2e\n", skytype, zenithbr, groundbr);
295          else
296 <                printf("7 -1 %.2e %.2e %.2e %f %f %f\n", zenithbr, groundbr,
297 <                                F2, sundir[0], sundir[1], sundir[2]);
296 >                printf("7 %d %.2e %.2e %.2e %f %f %f\n",
297 >                                skytype, zenithbr, groundbr, F2,
298 >                                sundir[0], sundir[1], sundir[2]);
299   }
300  
301  
302 < printdefaults()                 /* print default values */
302 > void
303 > printdefaults(void)                     /* print default values */
304   {
305 <        if (cloudy)
305 >        switch (skytype) {
306 >        case S_OVER:
307                  printf("-c\t\t\t\t# Cloudy sky\n");
308 <        else if (dosun)
309 <                printf("+s\t\t\t\t# Sunny sky with sun\n");
310 <        else
311 <                printf("-s\t\t\t\t# Sunny sky without sun\n");
308 >                break;
309 >        case S_UNIF:
310 >                printf("-u\t\t\t\t# Uniform cloudy sky\n");
311 >                break;
312 >        case S_INTER:
313 >                if (dosun)
314 >                        printf("+i\t\t\t\t# Intermediate sky with sun\n");
315 >                else
316 >                        printf("-i\t\t\t\t# Intermediate sky without sun\n");
317 >                break;
318 >        case S_CLEAR:
319 >                if (dosun)
320 >                        printf("+s\t\t\t\t# Sunny sky with sun\n");
321 >                else
322 >                        printf("-s\t\t\t\t# Sunny sky without sun\n");
323 >                break;
324 >        }
325          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
326          if (zenithbr > 0.0)
327                  printf("-b %f\t\t\t# Zenith radiance (watts/ster/m2\n", zenithbr);
# Line 205 | Line 333 | printdefaults()                        /* print default values */
333   }
334  
335  
336 < userror(msg)                    /* print usage error and quit */
337 < char  *msg;
336 > void
337 > userror(                        /* print usage error and quit */
338 >        char  *msg
339 > )
340   {
341          if (msg != NULL)
342                  fprintf(stderr, "%s: Use error - %s\n", progname, msg);
343          fprintf(stderr, "Usage: %s month day hour [options]\n", progname);
344 +        fprintf(stderr, "   Or: %s -ang altitude azimuth [options]\n", progname);
345          fprintf(stderr, "   Or: %s -defaults\n", progname);
346          exit(1);
347   }
348  
349  
350   double
351 < normsc(theta)                   /* compute normalization factor (E0*F2/L0) */
221 < double  theta;
351 > normsc(void)                    /* compute normalization factor (E0*F2/L0) */
352   {
353 <        static double  nf[5] = {2.766521, 0.547665,
354 <                                -0.369832, 0.009237, 0.059229};
353 >        static double  nfc[2][5] = {
354 >                                /* clear sky approx. */
355 >                {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
356 >                                /* intermediate sky approx. */
357 >                {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
358 >        };
359 >        register double  *nf;
360          double  x, nsc;
361          register int  i;
362                                          /* polynomial approximation */
363 <        x = (theta - PI/4.0)/(PI/4.0);
364 <        nsc = nf[4];
365 <        for (i = 3; i >= 0; i--)
363 >        nf = nfc[skytype==S_INTER];
364 >        x = (altitude - PI/4.0)/(PI/4.0);
365 >        nsc = nf[i=4];
366 >        while (i--)
367                  nsc = nsc*x + nf[i];
368  
369          return(nsc);
370   }
371  
372  
373 < printhead(ac, av)               /* print command header */
374 < register int  ac;
375 < register char  **av;
373 > void
374 > cvthour(                        /* convert hour string */
375 >        char  *hs
376 > )
377 > {
378 >        register char  *cp = hs;
379 >        register int    i, j;
380 >
381 >        if ( (tsolar = *cp == '+') ) cp++;              /* solar time? */
382 >        while (isdigit(*cp)) cp++;
383 >        if (*cp == ':')
384 >                hour = atoi(hs) + atoi(++cp)/60.0;
385 >        else {
386 >                hour = atof(hs);
387 >                if (*cp == '.') cp++;
388 >        }
389 >        while (isdigit(*cp)) cp++;
390 >        if (!*cp)
391 >                return;
392 >        if (tsolar || !isalpha(*cp)) {
393 >                fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
394 >                exit(1);
395 >        }
396 >        i = 0;
397 >        do {
398 >                for (j = 0; cp[j]; j++)
399 >                        if (toupper(cp[j]) != tzone[i].zname[j])
400 >                                break;
401 >                if (!cp[j] && !tzone[i].zname[j]) {
402 >                        s_meridian = tzone[i].zmer * (PI/180);
403 >                        return;
404 >                }
405 >        } while (tzone[i++].zname[0]);
406 >
407 >        fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
408 >        fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
409 >        for (i = 1; tzone[i].zname[0]; i++)
410 >                fprintf(stderr, " %s", tzone[i].zname);
411 >        putc('\n', stderr);
412 >        exit(1);
413 > }
414 >
415 >
416 > void
417 > printhead(              /* print command header */
418 >        register int  ac,
419 >        register char  **av
420 > )
421   {
422          putchar('#');
423          while (ac--) {

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines