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.6 by greg, Wed Oct 9 17:22:42 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 > #include "color.h"
2 > #ifndef lint
3 > static const char RCSid[] =
4 >    "$Id$";
5 > #endif
6 > /* Main function for generating spectral sky */
7 > /* Cloudy sky computed as weight average of clear and cie overcast sky */
8  
4 – #include "copyright.h"
9   #include "atmos.h"
10 + #include "copyright.h"
11   #include "resolu.h"
12 < #include "view.h"
12 > #include "rtio.h"
13 > #include <ctype.h>
14 > #ifdef _WIN32
15 > #include <windows.h>
16 > #else
17 > #include <errno.h>
18 > #include <sys/stat.h>
19 > #include <sys/types.h>
20 > #endif
21  
9 –
22   char *progname;
23  
24   const double ARCTIC_LAT = 67.;
# Line 17 | Line 29 | const double GNORM = 0.777778;
29  
30   const double D65EFF = 203.; /* standard illuminant D65 */
31  
32 < // Mean normalized relative daylight spectra where CCT = 6415K for overcast;
32 > /* Mean normalized relative daylight spectra where CCT = 6415K for overcast; */
33   const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
34                                1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
35                                1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
36                                0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
37  
38 + /* European and North American zones */
39 + struct {
40 +  char zname[8]; /* time zone name (all caps) */
41 +  float zmer;    /* standard meridian */
42 + } tzone[] = {{"YST", 135},   {"YDT", 120},   {"PST", 120},  {"PDT", 105},
43 +             {"MST", 105},   {"MDT", 90},    {"CST", 90},   {"CDT", 75},
44 +             {"EST", 75},    {"EDT", 60},    {"AST", 60},   {"ADT", 45},
45 +             {"NST", 52.5},  {"NDT", 37.5},  {"GMT", 0},    {"BST", -15},
46 +             {"CET", -15},   {"CEST", -30},  {"EET", -30},  {"EEST", -45},
47 +             {"AST", -45},   {"ADT", -60},   {"GST", -60},  {"GDT", -75},
48 +             {"IST", -82.5}, {"IDT", -97.5}, {"JST", -135}, {"NDT", -150},
49 +             {"NZST", -180}, {"NZDT", -195}, {"", 0}};
50 +
51 + static int make_directory(const char *path) {
52 + #ifdef _WIN32
53 +  if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
54 +    return 1;
55 +  }
56 +  return 0;
57 + #else
58 +  if (mkdir(path, 0777) == 0 || errno == EEXIST) {
59 +    return 1;
60 +  }
61 +  return 0;
62 + #endif
63 + }
64 +
65 + inline static float deg2rad(float deg) { return deg * (PI / 180.); }
66 +
67 + static int cvthour(char *hs, int *tsolar, double *hour) {
68 +  char *cp = hs;
69 +  int i, j;
70 +
71 +  if ((*tsolar = *cp == '+'))
72 +    cp++; /* solar time? */
73 +  while (isdigit(*cp))
74 +    cp++;
75 +  if (*cp == ':')
76 +    *hour = atoi(hs) + atoi(++cp) / 60.0;
77 +  else {
78 +    *hour = atof(hs);
79 +    if (*cp == '.')
80 +      cp++;
81 +  }
82 +  while (isdigit(*cp))
83 +    cp++;
84 +  if (!*cp)
85 +    return (0);
86 +  if (*tsolar || !isalpha(*cp)) {
87 +    fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
88 +    exit(1);
89 +  }
90 +  i = 0;
91 +  do {
92 +    for (j = 0; cp[j]; j++)
93 +      if (toupper(cp[j]) != tzone[i].zname[j])
94 +        break;
95 +    if (!cp[j] && !tzone[i].zname[j]) {
96 +      s_meridian = tzone[i].zmer * (PI / 180);
97 +      return (1);
98 +    }
99 +  } while (tzone[i++].zname[0]);
100 +
101 +  fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
102 +  fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
103 +  for (i = 1; tzone[i].zname[0]; i++)
104 +    fprintf(stderr, " %s", tzone[i].zname);
105 +  putc('\n', stderr);
106 +  exit(1);
107 + }
108 +
109 + static void basename(const char *path, char *output, size_t outsize) {
110 +  const char *last_slash = strrchr(path, '/');
111 +  const char *last_backslash = strrchr(path, '\\');
112 +  const char *filename = path;
113 +  const char *last_dot;
114 +
115 +  if (last_slash && last_backslash) {
116 +    filename =
117 +        (last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
118 +  } else if (last_slash) {
119 +    filename = last_slash + 1;
120 +  } else if (last_backslash) {
121 +    filename = last_backslash + 1;
122 +  }
123 +
124 +  last_dot = strrchr(filename, '.');
125 +  if (last_dot) {
126 +    size_t length = last_dot - filename;
127 +    if (length < outsize) {
128 +      strncpy(output, filename, length);
129 +      output[length] = '\0';
130 +    } else {
131 +      strncpy(output, filename, outsize - 1);
132 +      output[outsize - 1] = '\0';
133 +    }
134 +  }
135 + }
136 +
137 + static char *join_paths(const char *path1, const char *path2) {
138 +  size_t len1 = strlen(path1);
139 +  size_t len2 = strlen(path2);
140 +  int need_separator = (path1[len1 - 1] != DIRSEP);
141 +
142 +  char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
143 +  if (!result)
144 +    return NULL;
145 +
146 +  strcpy(result, path1);
147 +  if (need_separator) {
148 +    result[len1] = DIRSEP;
149 +    len1++;
150 +  }
151 +  strcpy(result + len1, path2);
152 +
153 +  return result;
154 + }
155 +
156   static inline double wmean2(const double a, const double b, const double x) {
157    return a * (1 - x) + b * x;
158   }
# Line 32 | Line 162 | static inline double wmean(const double a, const doubl
162    return (a * x + b * y) / (a + b);
163   }
164  
165 < static double get_zenith_brightness(const double sundir[3]) {
165 > static double get_overcast_zenith_brightness(const double sundir[3]) {
166    double zenithbr;
167    if (sundir[2] < 0) {
168      zenithbr = 0;
# Line 42 | Line 172 | static double get_zenith_brightness(const double sundi
172    return zenithbr;
173   }
174  
175 < // from gensky.c
175 > /* from gensky.c */
176   static double get_overcast_brightness(const double dz, const double zenithbr) {
177    double groundbr = zenithbr * GNORM;
178    return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
179                 pow(dz + 1.01, -10), groundbr);
180   }
181  
182 < static void write_rad_file(FILE *fp, const double *sun_radiance,
183 <                           const FVECT sundir, const char skyfile[PATH_MAX],
184 <                           const char grndfile[PATH_MAX]) {
182 > static void write_header(const int argc, char **argv, const double cloud_cover,
183 >                         const double grefl, const int res) {
184 >  int i;
185 >  printf("# ");
186 >  for (i = 0; i < argc; i++) {
187 >    printf("%s ", argv[i]);
188 >  }
189 >  printf("\n");
190 >  printf(
191 >      "#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
192 >      cloud_cover, grefl, res);
193 > }
194 >
195 > static void write_rad(const double *sun_radiance, const double intensity,
196 >                      const FVECT sundir, const char *ddir,
197 >                      const char *skyfile) {
198    if (sundir[2] > 0) {
199 <    fprintf(fp, "void spectrum sunrad\n0\n0\n22 380 780 ");
200 <    for (int i = 0; i < NSSAMP; ++i) {
201 <      fprintf(fp, "%.1f ", sun_radiance[i] * WVLSPAN);
199 >    printf("void spectrum sunrad\n0\n0\n22 380 780 ");
200 >    int i;
201 >    for (i = 0; i < NSSAMP; ++i) {
202 >      printf("%.3f ", sun_radiance[i]);
203      }
204 <    fprintf(fp, "\n\nsunrad light solar\n0\n0\n3 1 1 1\n\n");
205 <    fprintf(fp, "solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0],
206 <            sundir[1], sundir[2]);
204 >    printf("\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n", intensity,
205 >           intensity, intensity);
206 >    printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
207 >           sundir[2]);
208    }
209 <  fprintf(fp,
210 <          "void specpict skyfunc\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
211 <          "-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");
70 <
71 <  fprintf(fp,
72 <          "void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
73 <          "-my\n0\n0\n\n",
74 <          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");
209 >  printf("void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
210 >         "'1-Acos(Dz)/PI'\n0\n0\n\n",
211 >         skyfile);
212   }
213  
214   static void write_hsr_header(FILE *fp, RESOLU *res) {
215 <  float wvsplit[4] = {380, 480, 588,
81 <                      780}; // RGB wavelength limits+partitions (nm)
215 >  float wvsplit[4] = {380, 480, 588, 780};
216    newheader("RADIANCE", fp);
217    fputncomp(NSSAMP, fp);
218    fputwlsplit(wvsplit, fp);
# Line 87 | Line 221 | static void write_hsr_header(FILE *fp, RESOLU *res) {
221    fputsresolu(res, fp);
222   }
223  
224 + static inline float frac(float x) { return x - floor(x); }
225 +
226   int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
227                    DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
228                    const double cloud_cover, const FVECT sundir,
229 <                  const double grefl, const int res, const char *outname) {
230 <
95 <  char radfile[PATH_MAX];
229 >                  const double grefl, const int res, const char *outname,
230 >                  const char *ddir, const double dirnorm, const double difhor) {
231    char skyfile[PATH_MAX];
232 <  char grndfile[PATH_MAX];
233 <  if (!snprintf(radfile, sizeof(radfile), "%s.rad", outname)) {
99 <    fprintf(stderr, "Error setting rad file name\n");
100 <    return 0;
101 <  };
102 <  if (!snprintf(skyfile, sizeof(skyfile), "%s_sky.hsr", outname)) {
232 >  if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP,
233 >                outname)) {
234      fprintf(stderr, "Error setting sky file name\n");
235      return 0;
236    };
237 <  if (!snprintf(grndfile, sizeof(grndfile), "%s_ground.hsr", outname)) {
238 <    fprintf(stderr, "Error setting ground file name\n");
239 <    return 0;
109 <  }
110 <  RESOLU rs = {PIXSTANDARD, res, res};
237 >  int xres = res;
238 >  int yres = xres / 2;
239 >  RESOLU rs = {PIXSTANDARD, xres, yres};
240    FILE *skyfp = fopen(skyfile, "w");
112 –  FILE *grndfp = fopen(grndfile, "w");
113 –  write_hsr_header(grndfp, &rs);
241    write_hsr_header(skyfp, &rs);
115 –  VIEW skyview = {VT_ANG, {0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, 1.,
116 –                  180.,   180.,         0.,           0.,           0.,
117 –                  0.,     {0., 0., 0.}, {0., 0., 0.}, 0.,           0.};
118 –  VIEW grndview = {
119 –      VT_ANG, {0., 0., 0.}, {0., 0., -1.}, {0., 1., 0.}, 1., 180., 180., 0., 0.,
120 –      0.,     0.,           {0., 0., 0.},  {0., 0., 0.}, 0., 0.};
121 –  setview(&skyview);
122 –  setview(&grndview);
242  
243    CNDX[3] = NSSAMP;
244  
245 <  FVECT view_point = {0, 0, ER};
245 >  FVECT view_point = {0, 0, ER + 10};
246    const double radius = VLEN(view_point);
247    const double sun_ct = fdot(view_point, sundir) / radius;
248 <  for (unsigned int j = 0; j < res; ++j) {
249 <    for (unsigned int i = 0; i < res; ++i) {
250 <      RREAL loc[2];
251 <      FVECT rorg = {0};
252 <      FVECT rdir_sky = {0};
253 <      FVECT rdir_grnd = {0};
254 <      SCOLOR sky_radiance = {0};
255 <      SCOLOR ground_radiance = {0};
248 >
249 >  double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
250 >  double overcast_grndbr = overcast_zenithbr * GNORM;
251 >
252 >  double dif_ratio = 1;
253 >  if (difhor > 0) {
254 >    DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad_clear, radius, sun_ct);
255 >    double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
256 >    double diffuse_irradiance = 0;
257 >    int l;
258 >    for (l = 0; l < NSSAMP; ++l) {
259 >      diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20;  /* 20nm interval */
260 >    }
261 >    free(indirect_irradiance_clear);
262 >    diffuse_irradiance = wmean2(diffuse_irradiance, overcast_ghi, cloud_cover);
263 >    dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;       /* fudge */
264 >  }
265 >  int i, j, k;
266 >  for (j = 0; j < yres; ++j) {
267 >    for (i = 0; i < xres; ++i) {
268 >      SCOLOR radiance = {0};
269        SCOLR sky_sclr = {0};
138 –      SCOLR ground_sclr = {0};
270  
271 <      pix2loc(loc, &rs, i, j);
272 <      viewray(rorg, rdir_sky, &skyview, loc[0], loc[1]);
273 <      viewray(rorg, rdir_grnd, &grndview, loc[0], loc[1]);
271 >      float px = i / (xres - 1.0);
272 >      float py = j / (yres - 1.0);
273 >      float lambda = ((1 - py) * PI) - (PI / 2.0);
274 >      float phi = (px * 2.0 * PI) - PI;
275  
276 <      const double mu_sky = fdot(view_point, rdir_sky) / radius;
277 <      const double nu_sky = fdot(rdir_sky, sundir);
276 >      FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
277 >                    sin(lambda)};
278  
279 <      const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
280 <      const double nu_grnd = fdot(rdir_grnd, sundir);
279 >      const double mu = fdot(view_point, rdir) / radius;
280 >      const double nu = fdot(rdir, sundir);
281  
282 <      get_sky_radiance(scat_clear, scat1m_clear, radius, mu_sky, sun_ct, nu_sky,
283 <                       sky_radiance);
284 <      get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
285 <                          view_point, rdir_grnd, radius, mu_grnd, sun_ct,
286 <                          nu_grnd, grefl, sundir, ground_radiance);
282 >      /* hit ground */
283 >      if (rdir[2] < 0) {
284 >        get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
285 >                            view_point, rdir, radius, mu, sun_ct, nu, grefl,
286 >                            sundir, radiance);
287 >      } else {
288 >        get_sky_radiance(scat_clear, scat1m_clear, radius, mu, sun_ct, nu,
289 >                         radiance);
290 >      }
291  
292 <      for (int k = 0; k < NSSAMP; ++k) {
293 <        sky_radiance[k] *= WVLSPAN;
158 <        ground_radiance[k] *= WVLSPAN;
292 >      for (k = 0; k < NSSAMP; ++k) {
293 >        radiance[k] *= WVLSPAN;
294        }
295  
296        if (cloud_cover > 0) {
297 <        double zenithbr = get_zenith_brightness(sundir);
298 <        double grndbr = zenithbr * GNORM;
299 <        double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
300 <        for (int k = 0; k < NSSAMP; ++k) {
301 <          sky_radiance[k] =
302 <              wmean2(sky_radiance[k], skybr * D6415[k], cloud_cover);
303 <          ground_radiance[k] =
304 <              wmean2(ground_radiance[k], grndbr * D6415[k], cloud_cover);
297 >        double skybr = get_overcast_brightness(rdir[2], overcast_zenithbr);
298 >        if (rdir[2] < 0) {
299 >          for (k = 0; k < NSSAMP; ++k) {
300 >            radiance[k] = wmean2(radiance[k], overcast_grndbr * D6415[k], cloud_cover);
301 >          }
302 >        } else {
303 >          for (k = 0; k < NSSAMP; ++k) {
304 >            radiance[k] = wmean2(radiance[k], skybr * D6415[k], cloud_cover);
305 >          }
306          }
307        }
308  
309 <      scolor2scolr(sky_sclr, sky_radiance, 20);
310 <      putbinary(sky_sclr, LSCOLR, 1, skyfp);
309 >      for (k = 0; k < NSSAMP; ++k) {
310 >        radiance[k] *= dif_ratio;
311 >      }
312  
313 <      scolor2scolr(ground_sclr, ground_radiance, 20);
314 <      putbinary(ground_sclr, LSCOLR, 1, grndfp);
313 >      scolor2scolr(sky_sclr, radiance, NSSAMP);
314 >      putbinary(sky_sclr, LSCOLR, 1, skyfp);
315      }
316    }
317    fclose(skyfp);
181 –  fclose(grndfp);
318  
319 <  // Get solar radiance
319 >  /* Get solar radiance */
320    double sun_radiance[NSSAMP] = {0};
321    get_solar_radiance(tau_clear, scat_clear, scat1m_clear, sundir, radius,
322                       sun_ct, sun_radiance);
323    if (cloud_cover > 0) {
324 <    double zenithbr = get_zenith_brightness(sundir);
325 <    double skybr = get_overcast_brightness(sundir[2], zenithbr);
326 <    for (int i = 0; i < NSSAMP; ++i) {
324 >    double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
325 >    int i;
326 >    for (i = 0; i < NSSAMP; ++i) {
327        sun_radiance[i] =
328            wmean2(sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
329      }
330    }
331  
332 <  FILE *rfp = fopen(radfile, "w");
333 <  write_rad_file(rfp, sun_radiance, sundir, skyfile, grndfile);
334 <  fclose(rfp);
332 >  /* Normalize */
333 >  double sum = 0.0;
334 >  for (i = 0; i < NSSAMP; ++i) {
335 >    sum += sun_radiance[i];
336 >  }
337 >  double mean = sum / NSSAMP;
338 >  for (i = 0; i < NSSAMP; ++i) {
339 >    sun_radiance[i] /= mean;
340 >  }
341 >  double intensity = mean * WVLSPAN;
342 >  if (dirnorm > 0) {
343 >    intensity = dirnorm / SOLOMG / WHTEFFICACY;
344 >  }
345 >
346 >  write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
347    return 1;
348   }
349  
350 < static DpPaths get_dppaths(const double aod, const char *tag) {
350 > static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
351 >                           const char *tag) {
352    DpPaths paths;
353  
354 <  snprintf(paths.tau, PATH_MAX, "tau_%s_%.2f.dat", tag, aod);
355 <  snprintf(paths.scat, PATH_MAX, "scat_%s_%.2f.dat", tag, aod);
356 <  snprintf(paths.scat1m, PATH_MAX, "scat1m_%s_%.2f.dat", tag, aod);
357 <  snprintf(paths.irrad, PATH_MAX, "irrad_%s_%.2f.dat", tag, aod);
354 >  snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
355 >           mname, aod);
356 >  snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
357 >           mname, aod);
358 >  snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
359 >           tag, mname, aod);
360 >  snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
361 >           mname, aod);
362  
363    return paths;
364   }
365  
366 < static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag, const int is_summer,
366 > static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
367 >                                         const int is_summer,
368                                           const double s_latitude) {
369 <  // Set rayleigh density profile
216 <  if (fabs(s_latitude*180.0 / PI) > ARCTIC_LAT) {
369 >  if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
370      tag[0] = 's';
371      if (is_summer) {
372        tag[1] = 's';
# Line 224 | Line 377 | static void set_rayleigh_density_profile(Atmosphere *a
377        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
378        atmos->beta_r0 = BR0_SW;
379      }
380 <  } else if (fabs(s_latitude*180.0/PI) > TROPIC_LAT) {
380 >  } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
381      tag[0] = 'm';
382      if (is_summer) {
383        tag[1] = 's';
# Line 245 | Line 398 | static void set_rayleigh_density_profile(Atmosphere *a
398   }
399  
400   static Atmosphere init_atmos(const double aod, const double grefl) {
401 <  Atmosphere atmos = {
402 <      .ozone_density = {.layers =
403 <                            {
404 <                                {.width = 25000.0,
405 <                                 .exp_term = 0.0,
406 <                                 .exp_scale = 0.0,
407 <                                 .linear_term = 1.0 / 15000.0,
408 <                                 .constant_term = -2.0 / 3.0},
409 <                                {.width = AH,
410 <                                 .exp_term = 0.0,
411 <                                 .exp_scale = 0.0,
412 <                                 .linear_term = -1.0 / 15000.0,
413 <                                 .constant_term = 8.0 / 3.0},
414 <                            }},
415 <      .rayleigh_density = {.layers =
416 <                               {
417 <                                   {.width = AH,
418 <                                    .exp_term = 1.0,
419 <                                    .exp_scale = -1.0 / HR_MS,
420 <                                    .linear_term = 0.0,
421 <                                    .constant_term = 0.0},
422 <                               }},
423 <      .beta_r0 = BR0_MS,
424 <      .beta_scale = aod / AOD0_CA,
425 <      .beta_m = NULL,
273 <      .grefl = grefl
274 <  };
401 >  Atmosphere atmos = {.ozone_density = {.layers =
402 >                                            {
403 >                                                {.width = 25000.0,
404 >                                                 .exp_term = 0.0,
405 >                                                 .exp_scale = 0.0,
406 >                                                 .linear_term = 1.0 / 15000.0,
407 >                                                 .constant_term = -2.0 / 3.0},
408 >                                                {.width = AH,
409 >                                                 .exp_term = 0.0,
410 >                                                 .exp_scale = 0.0,
411 >                                                 .linear_term = -1.0 / 15000.0,
412 >                                                 .constant_term = 8.0 / 3.0},
413 >                                            }},
414 >                      .rayleigh_density = {.layers =
415 >                                               {
416 >                                                   {.width = AH,
417 >                                                    .exp_term = 1.0,
418 >                                                    .exp_scale = -1.0 / HR_MS,
419 >                                                    .linear_term = 0.0,
420 >                                                    .constant_term = 0.0},
421 >                                               }},
422 >                      .beta_r0 = BR0_MS,
423 >                      .beta_scale = aod / AOD0_CA,
424 >                      .beta_m = NULL,
425 >                      .grefl = grefl};
426    return atmos;
427   }
428  
# Line 284 | Line 435 | int main(int argc, char *argv[]) {
435    int sorder = 4;
436    int year = 0;
437    int tsolar = 0;
438 +  int got_meridian = 0;
439    double grefl = 0.2;
440    double ccover = 0.0;
441 <  int res = 128;
441 >  int res = 64;
442    double aod = AOD0_CA;
443    char *outname = "out";
444    char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
445 +  char mie_name[20] = "mie_ca";
446    char lstag[3];
447 +  char *ddir = ".";
448 +  int i;
449 +  double dirnorm = 0; /* direct normal illuminance */
450 +  double difhor = 0;  /* diffuse horizontal illuminance */
451  
452 +  if (argc == 2 && !strcmp(argv[1], "-defaults")) {
453 +    printf("-i %d\t\t\t\t#scattering order\n", sorder);
454 +    printf("-g %f\t\t\t#ground reflectance\n", grefl);
455 +    printf("-c %f\t\t\t#cloud cover\n", ccover);
456 +    printf("-r %d\t\t\t\t#image resolution\n", res);
457 +    printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
458 +    printf("-f %s\t\t\t\t#output name (-f)\n", outname);
459 +    printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
460 +    exit(0);
461 +  }
462 +
463    if (argc < 4) {
464 <    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",
464 >    fprintf(stderr,
465 >            "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r "
466 >            "res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
467 >            "-g grefl -f outpath\n",
468              argv[0]);
469      return 0;
470    }
471  
472    month = atoi(argv[1]);
473 +  if (month < 1 || month > 12) {
474 +    fprintf(stderr, "bad month");
475 +    exit(1);
476 +  }
477    day = atoi(argv[2]);
478 <  hour = atof(argv[3]);
478 >  if (day < 1 || day > 31) {
479 >    fprintf(stderr, "bad month");
480 >    exit(1);
481 >  }
482 >  got_meridian = cvthour(argv[3], &tsolar, &hour);
483  
484    if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
485      fprintf(stderr, "Cannot compute solar angle\n");
486      exit(1);
487    }
488  
489 <  for (int i = 4; i < argc; i++) {
489 >  for (i = 4; i < argc; i++) {
490      if (argv[i][0] == '-') {
491        switch (argv[i][1]) {
492        case 'a':
493          s_latitude = atof(argv[++i]) * (PI / 180.0);
494          break;
316 –      case 'g':
317 –        grefl = atof(argv[++i]);
318 –        break;
495        case 'c':
496          ccover = atof(argv[++i]);
497          break;
498        case 'd':
499          aod = atof(argv[++i]);
500          break;
501 +      case 'f':
502 +        outname = argv[++i];
503 +        break;
504 +      case 'g':
505 +        grefl = atof(argv[++i]);
506 +        break;
507        case 'i':
508          sorder = atoi(argv[++i]);
509          break;
510        case 'l':
511          mie_path = argv[++i];
512 +        basename(mie_path, mie_name, sizeof(mie_name));
513          break;
514        case 'm':
515 +        if (got_meridian) {
516 +          ++i;
517 +          break;
518 +        }
519          s_meridian = atof(argv[++i]) * (PI / 180.0);
520          break;
334 –      case 'o':
335 –        s_longitude = atof(argv[++i]) * (PI / 180.0);
336 –        break;
521        case 'n':
522          num_threads = atoi(argv[++i]);
523          break;
524 <      case 'y':
525 <        year = atoi(argv[++i]);
524 >      case 'o':
525 >        s_longitude = atof(argv[++i]) * (PI / 180.0);
526          break;
527 <      case 'f':
528 <        outname = argv[++i];
527 >      case 'L':
528 >        dirnorm = atof(argv[++i]);
529 >        difhor = atof(argv[++i]);
530          break;
531 +      case 'p':
532 +        ddir = argv[++i];
533 +        break;
534        case 'r':
535          res = atoi(argv[++i]);
536          break;
537 +      case 'y':
538 +        year = atoi(argv[++i]);
539 +        break;
540        default:
541          fprintf(stderr, "Unknown option %s\n", argv[i]);
542          exit(1);
543        }
544      }
545    }
546 +  if (year && (year < 1950) | (year > 2050))
547 +    fprintf(stderr, "%s: warning - year should be in range 1950-2050\n",
548 +            progname);
549 +  if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180)
550 +    fprintf(stderr,
551 +            "%s: warning - %.1f hours btwn. standard meridian and longitude\n",
552 +            progname, (s_longitude - s_meridian) * 12 / PI);
553  
554    Atmosphere clear_atmos = init_atmos(aod, grefl);
555  
# Line 361 | Line 559 | int main(int argc, char *argv[]) {
559    }
560    set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
561  
562 <  // Load mie density data
562 >  /* Load mie density data */
563    DATARRAY *mie_dp = getdata(mie_path);
564    if (mie_dp == NULL) {
565      fprintf(stderr, "Error reading mie data\n");
# Line 369 | Line 567 | int main(int argc, char *argv[]) {
567    }
568    clear_atmos.beta_m = mie_dp;
569  
570 <  DpPaths clear_paths = get_dppaths(aod, lstag);
570 >  char gsdir[PATH_MAX];
571 >  size_t siz = strlen(ddir);
572 >  if (ISDIRSEP(ddir[siz - 1]))
573 >    ddir[siz - 1] = '\0';
574 >  snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
575 >  if (!make_directory(gsdir)) {
576 >    fprintf(stderr, "Failed creating atmos_data directory");
577 >    exit(1);
578 >  }
579 >  DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
580  
581    if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
582        getpath(clear_paths.scat, ".", R_OK) == NULL ||
583        getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
584        getpath(clear_paths.irrad, ".", R_OK) == NULL) {
585 <    printf("# Precomputing...\n");
585 >    printf("# Pre-computing...\n");
586      if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
587 <      fprintf(stderr, "Precompute failed\n");
587 >      fprintf(stderr, "Pre-compute failed\n");
588        return 0;
589      }
590    }
# Line 387 | Line 594 | int main(int argc, char *argv[]) {
594    DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
595    DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
596  
597 +  write_header(argc, argv, ccover, grefl, res);
598 +
599    if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
600 <                     irrad_clear_dp, ccover, sundir, grefl, res, outname)) {
600 >                     irrad_clear_dp, ccover, sundir, grefl, res, outname, ddir,
601 >                     dirnorm, difhor)) {
602      fprintf(stderr, "gen_spect_sky failed\n");
603      exit(1);
604    }

Diff Legend

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