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

Diff Legend

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