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

Comparing ray/src/gen/genssky.c (file contents):
Revision 2.1 by greg, Fri Jul 5 18:04:36 2024 UTC vs.
Revision 2.2 by greg, Fri Jul 19 23:38:28 2024 UTC

# Line 1 | Line 1
1 < // Main function for generating spectral sky
2 < // Cloudy sky computed as weight average of clear and cie overcast sky
1 > #ifndef lint
2 > static const char RCSid[] = "$Id$";
3 > #endif
4 > /* Main function for generating spectral sky */
5 > /* Cloudy sky computed as weight average of clear and cie overcast sky */
6  
4 – #include "copyright.h"
7   #include "atmos.h"
8 + #include "copyright.h"
9   #include "resolu.h"
10 + #include "rtio.h"
11   #include "view.h"
12 + #include <ctype.h>
13 + #ifdef _WIN32
14 + #include <windows.h>
15 + #else
16 + #include <errno.h>
17 + #include <sys/stat.h>
18 + #include <sys/types.h>
19 + #endif
20  
9 –
21   char *progname;
22  
23   const double ARCTIC_LAT = 67.;
# Line 17 | Line 28 | const double GNORM = 0.777778;
28  
29   const double D65EFF = 203.; /* standard illuminant D65 */
30  
31 < // Mean normalized relative daylight spectra where CCT = 6415K for overcast;
31 > /* Mean normalized relative daylight spectra where CCT = 6415K for overcast; */
32   const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
33                                1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
34                                1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
35                                0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
36  
37 + /* European and North American zones */
38 + struct {
39 +  char zname[8]; /* time zone name (all caps) */
40 +  float zmer;    /* standard meridian */
41 + } tzone[] = {{"YST", 135},   {"YDT", 120},   {"PST", 120},  {"PDT", 105},
42 +             {"MST", 105},   {"MDT", 90},    {"CST", 90},   {"CDT", 75},
43 +             {"EST", 75},    {"EDT", 60},    {"AST", 60},   {"ADT", 45},
44 +             {"NST", 52.5},  {"NDT", 37.5},  {"GMT", 0},    {"BST", -15},
45 +             {"CET", -15},   {"CEST", -30},  {"EET", -30},  {"EEST", -45},
46 +             {"AST", -45},   {"ADT", -60},   {"GST", -60},  {"GDT", -75},
47 +             {"IST", -82.5}, {"IDT", -97.5}, {"JST", -135}, {"NDT", -150},
48 +             {"NZST", -180}, {"NZDT", -195}, {"", 0}};
49 +
50 + static int make_directory(const char *path) {
51 + #ifdef _WIN32
52 +  if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
53 +    return 1;
54 +  }
55 +  return 0;
56 + #else
57 +  if (mkdir(path, 0777) == 0 || errno == EEXIST) {
58 +    return 1;
59 +  }
60 +  return 0;
61 + #endif
62 + }
63 +
64 + static int cvthour(char *hs, int *tsolar, double *hour) {
65 +  char *cp = hs;
66 +  int i, j;
67 +
68 +  if ((*tsolar = *cp == '+'))
69 +    cp++; /* solar time? */
70 +  while (isdigit(*cp))
71 +    cp++;
72 +  if (*cp == ':')
73 +    *hour = atoi(hs) + atoi(++cp) / 60.0;
74 +  else {
75 +    *hour = atof(hs);
76 +    if (*cp == '.')
77 +      cp++;
78 +  }
79 +  while (isdigit(*cp))
80 +    cp++;
81 +  if (!*cp)
82 +    return (0);
83 +  if (*tsolar || !isalpha(*cp)) {
84 +    fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
85 +    exit(1);
86 +  }
87 +  i = 0;
88 +  do {
89 +    for (j = 0; cp[j]; j++)
90 +      if (toupper(cp[j]) != tzone[i].zname[j])
91 +        break;
92 +    if (!cp[j] && !tzone[i].zname[j]) {
93 +      s_meridian = tzone[i].zmer * (PI / 180);
94 +      return (1);
95 +    }
96 +  } while (tzone[i++].zname[0]);
97 +
98 +  fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
99 +  fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
100 +  for (i = 1; tzone[i].zname[0]; i++)
101 +    fprintf(stderr, " %s", tzone[i].zname);
102 +  putc('\n', stderr);
103 +  exit(1);
104 + }
105 +
106 + static void basename(const char *path, char *output, size_t outsize) {
107 +  const char *last_slash = strrchr(path, '/');
108 +  const char *last_backslash = strrchr(path, '\\');
109 +  const char *filename = path;
110 +  const char *last_dot;
111 +
112 +  if (last_slash && last_backslash) {
113 +    filename =
114 +        (last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
115 +  } else if (last_slash) {
116 +    filename = last_slash + 1;
117 +  } else if (last_backslash) {
118 +    filename = last_backslash + 1;
119 +  }
120 +
121 +  last_dot = strrchr(filename, '.');
122 +  if (last_dot) {
123 +    size_t length = last_dot - filename;
124 +    if (length < outsize) {
125 +      strncpy(output, filename, length);
126 +      output[length] = '\0';
127 +    } else {
128 +      strncpy(output, filename, outsize - 1);
129 +      output[outsize - 1] = '\0';
130 +    }
131 +  }
132 + }
133 +
134 + char *join_paths(const char *path1, const char *path2) {
135 +  size_t len1 = strlen(path1);
136 +  size_t len2 = strlen(path2);
137 +  int need_separator = (path1[len1 - 1] != DIRSEP);
138 +
139 +  char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
140 +  if (!result)
141 +    return NULL;
142 +
143 +  strcpy(result, path1);
144 +  if (need_separator) {
145 +    result[len1] = DIRSEP;
146 +    len1++;
147 +  }
148 +  strcpy(result + len1, path2);
149 +
150 +  return result;
151 + }
152 +
153   static inline double wmean2(const double a, const double b, const double x) {
154    return a * (1 - x) + b * x;
155   }
# Line 42 | Line 169 | static double get_zenith_brightness(const double sundi
169    return zenithbr;
170   }
171  
172 < // from gensky.c
172 > /* from gensky.c */
173   static double get_overcast_brightness(const double dz, const double zenithbr) {
174    double groundbr = zenithbr * GNORM;
175    return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
176                 pow(dz + 1.01, -10), groundbr);
177   }
178  
179 < static void write_rad_file(FILE *fp, const double *sun_radiance,
180 <                           const FVECT sundir, const char skyfile[PATH_MAX],
181 <                           const char grndfile[PATH_MAX]) {
179 > static void write_header(const int argc, char **argv, const double cloud_cover,
180 >                         const double grefl, const int res) {
181 >  printf("# ");
182 >  for (int i = 0; i < argc; i++) {
183 >    printf("%s ", argv[i]);
184 >  }
185 >  printf("\n");
186 >  printf("#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
187 >         cloud_cover, grefl, res);
188 > }
189 >
190 > static void write_rad(const double *sun_radiance, const FVECT sundir,
191 >                      const char skyfile[PATH_MAX],
192 >                      const char grndfile[PATH_MAX]) {
193    if (sundir[2] > 0) {
194 <    fprintf(fp, "void spectrum sunrad\n0\n0\n22 380 780 ");
194 >    printf("void spectrum sunrad\n0\n0\n22 380 780 ");
195 >    /* Normalize to one */
196 >    double sum = 0.0;
197      for (int i = 0; i < NSSAMP; ++i) {
198 <      fprintf(fp, "%.1f ", sun_radiance[i] * WVLSPAN);
198 >      sum += sun_radiance[i];
199      }
200 <    fprintf(fp, "\n\nsunrad light solar\n0\n0\n3 1 1 1\n\n");
201 <    fprintf(fp, "solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0],
202 <            sundir[1], sundir[2]);
200 >    double mean = sum / NSSAMP;
201 >    for (int i = 0; i < NSSAMP; ++i) {
202 >      printf("%.3f ", sun_radiance[i] / mean);
203 >    }
204 >    double intensity = mean * WVLSPAN;
205 >    printf("\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n", intensity,
206 >           intensity, intensity);
207 >    printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
208 >           sundir[2]);
209    }
210 <  fprintf(fp,
211 <          "void specpict skyfunc\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
212 <          "-mx\n0\n0\n\n",
67 <          skyfile);
68 <  fprintf(fp, "skyfunc glow sky_glow\n0\n0\n4 1 1 1 0\n\n");
69 <  fprintf(fp, "sky_glow source sky\n0\n0\n4 0 0 1 180\n\n");
210 >  printf("void specpict skymap\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
211 >         "-mx\n0\n0\n\n",
212 >         skyfile);
213  
214 <  fprintf(fp,
215 <          "void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
216 <          "-my\n0\n0\n\n",
217 <          grndfile);
75 <  fprintf(fp, "grndmap glow ground_glow\n0\n0\n4 1 1 1 0\n\n");
76 <  fprintf(fp, "ground_glow source ground_source\n0\n0\n4 0 0 -1 180\n\n");
214 >  printf("void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
215 >         "-my\n0\n0\n\n",
216 >         grndfile);
217 >  printf("void mixfunc skyfunc\n4 skymap grndmap if(Dz,1,0) .\n0\n0\n");
218   }
219  
220   static void write_hsr_header(FILE *fp, RESOLU *res) {
221    float wvsplit[4] = {380, 480, 588,
222 <                      780}; // RGB wavelength limits+partitions (nm)
222 >                      780}; /* RGB wavelength limits+partitions (nm) */
223    newheader("RADIANCE", fp);
224    fputncomp(NSSAMP, fp);
225    fputwlsplit(wvsplit, fp);
# Line 91 | Line 232 | int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_
232                    DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
233                    const double cloud_cover, const FVECT sundir,
234                    const double grefl, const int res, const char *outname) {
94 –
95 –  char radfile[PATH_MAX];
235    char skyfile[PATH_MAX];
236    char grndfile[PATH_MAX];
98 –  if (!snprintf(radfile, sizeof(radfile), "%s.rad", outname)) {
99 –    fprintf(stderr, "Error setting rad file name\n");
100 –    return 0;
101 –  };
237    if (!snprintf(skyfile, sizeof(skyfile), "%s_sky.hsr", outname)) {
238      fprintf(stderr, "Error setting sky file name\n");
239      return 0;
# Line 180 | Line 315 | int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_
315    fclose(skyfp);
316    fclose(grndfp);
317  
318 <  // Get solar radiance
318 >  /* Get solar radiance */
319    double sun_radiance[NSSAMP] = {0};
320    get_solar_radiance(tau_clear, scat_clear, scat1m_clear, sundir, radius,
321                       sun_ct, sun_radiance);
# Line 193 | Line 328 | int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_
328      }
329    }
330  
331 <  FILE *rfp = fopen(radfile, "w");
197 <  write_rad_file(rfp, sun_radiance, sundir, skyfile, grndfile);
198 <  fclose(rfp);
331 >  write_rad(sun_radiance, sundir, skyfile, grndfile);
332    return 1;
333   }
334  
335 < static DpPaths get_dppaths(const double aod, const char *tag) {
335 > static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
336 >                           const char *tag) {
337    DpPaths paths;
338  
339 <  snprintf(paths.tau, PATH_MAX, "tau_%s_%.2f.dat", tag, aod);
340 <  snprintf(paths.scat, PATH_MAX, "scat_%s_%.2f.dat", tag, aod);
341 <  snprintf(paths.scat1m, PATH_MAX, "scat1m_%s_%.2f.dat", tag, aod);
342 <  snprintf(paths.irrad, PATH_MAX, "irrad_%s_%.2f.dat", tag, aod);
339 >  snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
340 >           mname, aod);
341 >  snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
342 >           mname, aod);
343 >  snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
344 >           tag, mname, aod);
345 >  snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
346 >           mname, aod);
347  
348    return paths;
349   }
350  
351 < static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag, const int is_summer,
351 > static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
352 >                                         const int is_summer,
353                                           const double s_latitude) {
354 <  // Set rayleigh density profile
355 <  if (fabs(s_latitude*180.0 / PI) > ARCTIC_LAT) {
354 >  /* Set rayleigh density profile */
355 >  if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
356      tag[0] = 's';
357      if (is_summer) {
358        tag[1] = 's';
# Line 224 | Line 363 | static void set_rayleigh_density_profile(Atmosphere *a
363        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
364        atmos->beta_r0 = BR0_SW;
365      }
366 <  } else if (fabs(s_latitude*180.0/PI) > TROPIC_LAT) {
366 >  } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
367      tag[0] = 'm';
368      if (is_summer) {
369        tag[1] = 's';
# Line 245 | Line 384 | static void set_rayleigh_density_profile(Atmosphere *a
384   }
385  
386   static Atmosphere init_atmos(const double aod, const double grefl) {
387 <  Atmosphere atmos = {
388 <      .ozone_density = {.layers =
389 <                            {
390 <                                {.width = 25000.0,
391 <                                 .exp_term = 0.0,
392 <                                 .exp_scale = 0.0,
393 <                                 .linear_term = 1.0 / 15000.0,
394 <                                 .constant_term = -2.0 / 3.0},
395 <                                {.width = AH,
396 <                                 .exp_term = 0.0,
397 <                                 .exp_scale = 0.0,
398 <                                 .linear_term = -1.0 / 15000.0,
399 <                                 .constant_term = 8.0 / 3.0},
400 <                            }},
401 <      .rayleigh_density = {.layers =
402 <                               {
403 <                                   {.width = AH,
404 <                                    .exp_term = 1.0,
405 <                                    .exp_scale = -1.0 / HR_MS,
406 <                                    .linear_term = 0.0,
407 <                                    .constant_term = 0.0},
408 <                               }},
409 <      .beta_r0 = BR0_MS,
410 <      .beta_scale = aod / AOD0_CA,
411 <      .beta_m = NULL,
273 <      .grefl = grefl
274 <  };
387 >  Atmosphere atmos = {.ozone_density = {.layers =
388 >                                            {
389 >                                                {.width = 25000.0,
390 >                                                 .exp_term = 0.0,
391 >                                                 .exp_scale = 0.0,
392 >                                                 .linear_term = 1.0 / 15000.0,
393 >                                                 .constant_term = -2.0 / 3.0},
394 >                                                {.width = AH,
395 >                                                 .exp_term = 0.0,
396 >                                                 .exp_scale = 0.0,
397 >                                                 .linear_term = -1.0 / 15000.0,
398 >                                                 .constant_term = 8.0 / 3.0},
399 >                                            }},
400 >                      .rayleigh_density = {.layers =
401 >                                               {
402 >                                                   {.width = AH,
403 >                                                    .exp_term = 1.0,
404 >                                                    .exp_scale = -1.0 / HR_MS,
405 >                                                    .linear_term = 0.0,
406 >                                                    .constant_term = 0.0},
407 >                                               }},
408 >                      .beta_r0 = BR0_MS,
409 >                      .beta_scale = aod / AOD0_CA,
410 >                      .beta_m = NULL,
411 >                      .grefl = grefl};
412    return atmos;
413   }
414  
# Line 284 | Line 421 | int main(int argc, char *argv[]) {
421    int sorder = 4;
422    int year = 0;
423    int tsolar = 0;
424 +  int got_meridian = 0;
425    double grefl = 0.2;
426    double ccover = 0.0;
427    int res = 128;
428    double aod = AOD0_CA;
429    char *outname = "out";
430    char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
431 +  char mie_name[20] = "mie_ca";
432    char lstag[3];
433 +  char *ddir = ".";
434  
435 +  if (!strcmp(argv[1], "-defaults")) {
436 +    printf("-i %d\t\t\t\t#scattering order\n", sorder);
437 +    printf("-g %f\t\t\t#ground reflectance\n", grefl);
438 +    printf("-c %f\t\t\t#cloud cover\n", ccover);
439 +    printf("-r %d\t\t\t\t#image resolution\n", res);
440 +    printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
441 +    printf("-f %s\t\t\t\t#output name (-f)\n", outname);
442 +    printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
443 +    exit(1);
444 +  }
445 +
446    if (argc < 4) {
447 <    fprintf(stderr, "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r res -n nproc -c ccover -l mie -g grefl -f outpath\n",
447 >    fprintf(stderr,
448 >            "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r "
449 >            "res -n nproc -c ccover -l mie -g grefl -f outpath\n",
450              argv[0]);
451      return 0;
452    }
453  
454    month = atoi(argv[1]);
455 +  if (month < 1 || month > 12) {
456 +    fprintf(stderr, "bad month");
457 +    exit(1);
458 +  }
459    day = atoi(argv[2]);
460 <  hour = atof(argv[3]);
460 >  if (day < 1 || day > 31) {
461 >    fprintf(stderr, "bad month");
462 >    exit(1);
463 >  }
464 >  got_meridian = cvthour(argv[3], &tsolar, &hour);
465  
466    if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
467      fprintf(stderr, "Cannot compute solar angle\n");
# Line 313 | Line 474 | int main(int argc, char *argv[]) {
474        case 'a':
475          s_latitude = atof(argv[++i]) * (PI / 180.0);
476          break;
316 –      case 'g':
317 –        grefl = atof(argv[++i]);
318 –        break;
477        case 'c':
478          ccover = atof(argv[++i]);
479          break;
480        case 'd':
481          aod = atof(argv[++i]);
482          break;
483 +      case 'f':
484 +        outname = argv[++i];
485 +        break;
486 +      case 'g':
487 +        grefl = atof(argv[++i]);
488 +        break;
489        case 'i':
490          sorder = atoi(argv[++i]);
491          break;
492        case 'l':
493          mie_path = argv[++i];
494 +        basename(mie_path, mie_name, sizeof(mie_name));
495          break;
496        case 'm':
497 +        if (got_meridian) {
498 +          ++i;
499 +          break;
500 +        }
501          s_meridian = atof(argv[++i]) * (PI / 180.0);
502          break;
334 –      case 'o':
335 –        s_longitude = atof(argv[++i]) * (PI / 180.0);
336 –        break;
503        case 'n':
504          num_threads = atoi(argv[++i]);
505          break;
506 <      case 'y':
507 <        year = atoi(argv[++i]);
506 >      case 'o':
507 >        s_longitude = atof(argv[++i]) * (PI / 180.0);
508          break;
509 <      case 'f':
510 <        outname = argv[++i];
509 >      case 'p':
510 >        ddir = argv[++i];
511          break;
512        case 'r':
513          res = atoi(argv[++i]);
514          break;
515 +      case 'y':
516 +        year = atoi(argv[++i]);
517 +        break;
518        default:
519          fprintf(stderr, "Unknown option %s\n", argv[i]);
520          exit(1);
521        }
522      }
523    }
524 +  if (year && (year < 1950) | (year > 2050))
525 +    fprintf(stderr, "%s: warning - year should be in range 1950-2050\n",
526 +            progname);
527 +  if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180)
528 +    fprintf(stderr,
529 +            "%s: warning - %.1f hours btwn. standard meridian and longitude\n",
530 +            progname, (s_longitude - s_meridian) * 12 / PI);
531  
532    Atmosphere clear_atmos = init_atmos(aod, grefl);
533  
# Line 361 | Line 537 | int main(int argc, char *argv[]) {
537    }
538    set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
539  
540 <  // Load mie density data
540 >  /* Load mie density data */
541    DATARRAY *mie_dp = getdata(mie_path);
542    if (mie_dp == NULL) {
543      fprintf(stderr, "Error reading mie data\n");
# Line 369 | Line 545 | int main(int argc, char *argv[]) {
545    }
546    clear_atmos.beta_m = mie_dp;
547  
548 <  DpPaths clear_paths = get_dppaths(aod, lstag);
548 >  char gsdir[PATH_MAX];
549 >  size_t siz = strlen(ddir);
550 >  if (ISDIRSEP(ddir[siz-1]))
551 >    ddir[siz-1] = '\0';
552 >  snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
553 >  printf("gsdir: %s\n", gsdir);
554 >  if (!make_directory(gsdir)) {
555 >    fprintf(stderr, "Failed creating atmos_data directory");
556 >    exit(1);
557 >  }
558 >  DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
559  
560    if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
561        getpath(clear_paths.scat, ".", R_OK) == NULL ||
562        getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
563        getpath(clear_paths.irrad, ".", R_OK) == NULL) {
564 <    printf("# Precomputing...\n");
564 >    printf("# Pre-computing...\n");
565      if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
566 <      fprintf(stderr, "Precompute failed\n");
566 >      fprintf(stderr, "Pre-compute failed\n");
567        return 0;
568      }
569    }
# Line 386 | Line 572 | int main(int argc, char *argv[]) {
572    DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
573    DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
574    DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
575 +
576 +  write_header(argc, argv, ccover, grefl, res);
577  
578    if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
579                       irrad_clear_dp, ccover, sundir, grefl, res, outname)) {

Diff Legend

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