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.7 by greg, Thu Apr 10 23:30:58 2025 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 >    if (diffuse_irradiance > 0) {
264 >        dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;       /* fudge */
265 >    }
266 >  }
267 >  int i, j, k;
268 >  for (j = 0; j < yres; ++j) {
269 >    for (i = 0; i < xres; ++i) {
270 >      SCOLOR radiance = {0};
271        SCOLR sky_sclr = {0};
138 –      SCOLR ground_sclr = {0};
272  
273 <      pix2loc(loc, &rs, i, j);
274 <      viewray(rorg, rdir_sky, &skyview, loc[0], loc[1]);
275 <      viewray(rorg, rdir_grnd, &grndview, loc[0], loc[1]);
273 >      float px = i / (xres - 1.0);
274 >      float py = j / (yres - 1.0);
275 >      float lambda = ((1 - py) * PI) - (PI / 2.0);
276 >      float phi = (px * 2.0 * PI) - PI;
277  
278 <      const double mu_sky = fdot(view_point, rdir_sky) / radius;
279 <      const double nu_sky = fdot(rdir_sky, sundir);
278 >      FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
279 >                    sin(lambda)};
280  
281 <      const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
282 <      const double nu_grnd = fdot(rdir_grnd, sundir);
281 >      const double mu = fdot(view_point, rdir) / radius;
282 >      const double nu = fdot(rdir, sundir);
283  
284 <      get_sky_radiance(scat_clear, scat1m_clear, radius, mu_sky, sun_ct, nu_sky,
285 <                       sky_radiance);
286 <      get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
287 <                          view_point, rdir_grnd, radius, mu_grnd, sun_ct,
288 <                          nu_grnd, grefl, sundir, ground_radiance);
284 >      /* hit ground */
285 >      if (rdir[2] < 0) {
286 >        get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
287 >                            view_point, rdir, radius, mu, sun_ct, nu, grefl,
288 >                            sundir, radiance);
289 >      } else {
290 >        get_sky_radiance(scat_clear, scat1m_clear, radius, mu, sun_ct, nu,
291 >                         radiance);
292 >      }
293  
294 <      for (int k = 0; k < NSSAMP; ++k) {
295 <        sky_radiance[k] *= WVLSPAN;
158 <        ground_radiance[k] *= WVLSPAN;
294 >      for (k = 0; k < NSSAMP; ++k) {
295 >        radiance[k] *= WVLSPAN;
296        }
297  
298        if (cloud_cover > 0) {
299 <        double zenithbr = get_zenith_brightness(sundir);
300 <        double grndbr = zenithbr * GNORM;
301 <        double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
302 <        for (int k = 0; k < NSSAMP; ++k) {
303 <          sky_radiance[k] =
304 <              wmean2(sky_radiance[k], skybr * D6415[k], cloud_cover);
305 <          ground_radiance[k] =
306 <              wmean2(ground_radiance[k], grndbr * D6415[k], cloud_cover);
299 >        double skybr = get_overcast_brightness(rdir[2], overcast_zenithbr);
300 >        if (rdir[2] < 0) {
301 >          for (k = 0; k < NSSAMP; ++k) {
302 >            radiance[k] = wmean2(radiance[k], overcast_grndbr * D6415[k], cloud_cover);
303 >          }
304 >        } else {
305 >          for (k = 0; k < NSSAMP; ++k) {
306 >            radiance[k] = wmean2(radiance[k], skybr * D6415[k], cloud_cover);
307 >          }
308          }
309        }
310  
311 <      scolor2scolr(sky_sclr, sky_radiance, 20);
312 <      putbinary(sky_sclr, LSCOLR, 1, skyfp);
311 >      for (k = 0; k < NSSAMP; ++k) {
312 >        radiance[k] *= dif_ratio;
313 >      }
314  
315 <      scolor2scolr(ground_sclr, ground_radiance, 20);
316 <      putbinary(ground_sclr, LSCOLR, 1, grndfp);
315 >      scolor2scolr(sky_sclr, radiance, NSSAMP);
316 >      putbinary(sky_sclr, LSCOLR, 1, skyfp);
317      }
318    }
319    fclose(skyfp);
181 –  fclose(grndfp);
320  
321 <  // Get solar radiance
321 >  /* Get solar radiance */
322    double sun_radiance[NSSAMP] = {0};
323    get_solar_radiance(tau_clear, scat_clear, scat1m_clear, sundir, radius,
324                       sun_ct, sun_radiance);
325    if (cloud_cover > 0) {
326 <    double zenithbr = get_zenith_brightness(sundir);
327 <    double skybr = get_overcast_brightness(sundir[2], zenithbr);
328 <    for (int i = 0; i < NSSAMP; ++i) {
326 >    double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
327 >    int i;
328 >    for (i = 0; i < NSSAMP; ++i) {
329        sun_radiance[i] =
330            wmean2(sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
331      }
332    }
333  
334 <  FILE *rfp = fopen(radfile, "w");
335 <  write_rad_file(rfp, sun_radiance, sundir, skyfile, grndfile);
336 <  fclose(rfp);
334 >  /* Normalize */
335 >  double sum = 0.0;
336 >  for (i = 0; i < NSSAMP; ++i) {
337 >    sum += sun_radiance[i];
338 >  }
339 >  double mean = sum / NSSAMP;
340 >  for (i = 0; i < NSSAMP; ++i) {
341 >    sun_radiance[i] /= mean;
342 >  }
343 >  double intensity = mean * WVLSPAN;
344 >  if (dirnorm > 0) {
345 >    intensity = dirnorm / SOLOMG / WHTEFFICACY;
346 >  }
347 >
348 >  write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
349    return 1;
350   }
351  
352 < static DpPaths get_dppaths(const double aod, const char *tag) {
352 > static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
353 >                           const char *tag) {
354    DpPaths paths;
355  
356 <  snprintf(paths.tau, PATH_MAX, "tau_%s_%.2f.dat", tag, aod);
357 <  snprintf(paths.scat, PATH_MAX, "scat_%s_%.2f.dat", tag, aod);
358 <  snprintf(paths.scat1m, PATH_MAX, "scat1m_%s_%.2f.dat", tag, aod);
359 <  snprintf(paths.irrad, PATH_MAX, "irrad_%s_%.2f.dat", tag, aod);
356 >  snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
357 >           mname, aod);
358 >  snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
359 >           mname, aod);
360 >  snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
361 >           tag, mname, aod);
362 >  snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
363 >           mname, aod);
364  
365    return paths;
366   }
367  
368 < static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag, const int is_summer,
368 > static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
369 >                                         const int is_summer,
370                                           const double s_latitude) {
371 <  // Set rayleigh density profile
216 <  if (fabs(s_latitude*180.0 / PI) > ARCTIC_LAT) {
371 >  if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
372      tag[0] = 's';
373      if (is_summer) {
374        tag[1] = 's';
# Line 224 | Line 379 | static void set_rayleigh_density_profile(Atmosphere *a
379        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
380        atmos->beta_r0 = BR0_SW;
381      }
382 <  } else if (fabs(s_latitude*180.0/PI) > TROPIC_LAT) {
382 >  } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
383      tag[0] = 'm';
384      if (is_summer) {
385        tag[1] = 's';
# Line 245 | Line 400 | static void set_rayleigh_density_profile(Atmosphere *a
400   }
401  
402   static Atmosphere init_atmos(const double aod, const double grefl) {
403 <  Atmosphere atmos = {
404 <      .ozone_density = {.layers =
405 <                            {
406 <                                {.width = 25000.0,
407 <                                 .exp_term = 0.0,
408 <                                 .exp_scale = 0.0,
409 <                                 .linear_term = 1.0 / 15000.0,
410 <                                 .constant_term = -2.0 / 3.0},
411 <                                {.width = AH,
412 <                                 .exp_term = 0.0,
413 <                                 .exp_scale = 0.0,
414 <                                 .linear_term = -1.0 / 15000.0,
415 <                                 .constant_term = 8.0 / 3.0},
416 <                            }},
417 <      .rayleigh_density = {.layers =
418 <                               {
419 <                                   {.width = AH,
420 <                                    .exp_term = 1.0,
421 <                                    .exp_scale = -1.0 / HR_MS,
422 <                                    .linear_term = 0.0,
423 <                                    .constant_term = 0.0},
424 <                               }},
425 <      .beta_r0 = BR0_MS,
426 <      .beta_scale = aod / AOD0_CA,
427 <      .beta_m = NULL,
273 <      .grefl = grefl
274 <  };
403 >  Atmosphere atmos = {.ozone_density = {.layers =
404 >                                            {
405 >                                                {.width = 25000.0,
406 >                                                 .exp_term = 0.0,
407 >                                                 .exp_scale = 0.0,
408 >                                                 .linear_term = 1.0 / 15000.0,
409 >                                                 .constant_term = -2.0 / 3.0},
410 >                                                {.width = AH,
411 >                                                 .exp_term = 0.0,
412 >                                                 .exp_scale = 0.0,
413 >                                                 .linear_term = -1.0 / 15000.0,
414 >                                                 .constant_term = 8.0 / 3.0},
415 >                                            }},
416 >                      .rayleigh_density = {.layers =
417 >                                               {
418 >                                                   {.width = AH,
419 >                                                    .exp_term = 1.0,
420 >                                                    .exp_scale = -1.0 / HR_MS,
421 >                                                    .linear_term = 0.0,
422 >                                                    .constant_term = 0.0},
423 >                                               }},
424 >                      .beta_r0 = BR0_MS,
425 >                      .beta_scale = aod / AOD0_CA,
426 >                      .beta_m = NULL,
427 >                      .grefl = grefl};
428    return atmos;
429   }
430  
# Line 284 | Line 437 | int main(int argc, char *argv[]) {
437    int sorder = 4;
438    int year = 0;
439    int tsolar = 0;
440 +  int got_meridian = 0;
441    double grefl = 0.2;
442    double ccover = 0.0;
443 <  int res = 128;
443 >  int res = 64;
444    double aod = AOD0_CA;
445    char *outname = "out";
446    char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
447 +  char mie_name[20] = "mie_ca";
448    char lstag[3];
449 +  char *ddir = ".";
450 +  int i;
451 +  double dirnorm = 0; /* direct normal illuminance */
452 +  double difhor = 0;  /* diffuse horizontal illuminance */
453  
454 +  if (argc == 2 && !strcmp(argv[1], "-defaults")) {
455 +    printf("-i %d\t\t\t\t#scattering order\n", sorder);
456 +    printf("-g %f\t\t\t#ground reflectance\n", grefl);
457 +    printf("-c %f\t\t\t#cloud cover\n", ccover);
458 +    printf("-r %d\t\t\t\t#image resolution\n", res);
459 +    printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
460 +    printf("-f %s\t\t\t\t#output name (-f)\n", outname);
461 +    printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
462 +    exit(0);
463 +  }
464 +
465    if (argc < 4) {
466 <    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",
466 >    fprintf(stderr,
467 >            "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r "
468 >            "res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
469 >            "-g grefl -f outpath\n",
470              argv[0]);
471      return 0;
472    }
473  
474    month = atoi(argv[1]);
475 +  if (month < 1 || month > 12) {
476 +    fprintf(stderr, "bad month");
477 +    exit(1);
478 +  }
479    day = atoi(argv[2]);
480 <  hour = atof(argv[3]);
480 >  if (day < 1 || day > 31) {
481 >    fprintf(stderr, "bad month");
482 >    exit(1);
483 >  }
484 >  got_meridian = cvthour(argv[3], &tsolar, &hour);
485  
486    if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
487      fprintf(stderr, "Cannot compute solar angle\n");
488      exit(1);
489    }
490  
491 <  for (int i = 4; i < argc; i++) {
491 >  for (i = 4; i < argc; i++) {
492      if (argv[i][0] == '-') {
493        switch (argv[i][1]) {
494        case 'a':
495          s_latitude = atof(argv[++i]) * (PI / 180.0);
496          break;
316 –      case 'g':
317 –        grefl = atof(argv[++i]);
318 –        break;
497        case 'c':
498          ccover = atof(argv[++i]);
499          break;
500        case 'd':
501          aod = atof(argv[++i]);
502          break;
503 +      case 'f':
504 +        outname = argv[++i];
505 +        break;
506 +      case 'g':
507 +        grefl = atof(argv[++i]);
508 +        break;
509        case 'i':
510          sorder = atoi(argv[++i]);
511          break;
512        case 'l':
513          mie_path = argv[++i];
514 +        basename(mie_path, mie_name, sizeof(mie_name));
515          break;
516        case 'm':
517 +        if (got_meridian) {
518 +          ++i;
519 +          break;
520 +        }
521          s_meridian = atof(argv[++i]) * (PI / 180.0);
522          break;
334 –      case 'o':
335 –        s_longitude = atof(argv[++i]) * (PI / 180.0);
336 –        break;
523        case 'n':
524          num_threads = atoi(argv[++i]);
525          break;
526 <      case 'y':
527 <        year = atoi(argv[++i]);
526 >      case 'o':
527 >        s_longitude = atof(argv[++i]) * (PI / 180.0);
528          break;
529 <      case 'f':
530 <        outname = argv[++i];
529 >      case 'L':
530 >        dirnorm = atof(argv[++i]);
531 >        difhor = atof(argv[++i]);
532          break;
533 +      case 'p':
534 +        ddir = argv[++i];
535 +        break;
536        case 'r':
537          res = atoi(argv[++i]);
538          break;
539 +      case 'y':
540 +        year = atoi(argv[++i]);
541 +        break;
542        default:
543          fprintf(stderr, "Unknown option %s\n", argv[i]);
544          exit(1);
545        }
546      }
547    }
548 +  if (year && (year < 1950) | (year > 2050))
549 +    fprintf(stderr, "%s: warning - year should be in range 1950-2050\n",
550 +            progname);
551 +  if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180)
552 +    fprintf(stderr,
553 +            "%s: warning - %.1f hours btwn. standard meridian and longitude\n",
554 +            progname, (s_longitude - s_meridian) * 12 / PI);
555  
556    Atmosphere clear_atmos = init_atmos(aod, grefl);
557  
# Line 361 | Line 561 | int main(int argc, char *argv[]) {
561    }
562    set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
563  
564 <  // Load mie density data
564 >  /* Load mie density data */
565    DATARRAY *mie_dp = getdata(mie_path);
566    if (mie_dp == NULL) {
567      fprintf(stderr, "Error reading mie data\n");
# Line 369 | Line 569 | int main(int argc, char *argv[]) {
569    }
570    clear_atmos.beta_m = mie_dp;
571  
572 <  DpPaths clear_paths = get_dppaths(aod, lstag);
572 >  char gsdir[PATH_MAX];
573 >  size_t siz = strlen(ddir);
574 >  if (ISDIRSEP(ddir[siz - 1]))
575 >    ddir[siz - 1] = '\0';
576 >  snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
577 >  if (!make_directory(gsdir)) {
578 >    fprintf(stderr, "Failed creating atmos_data directory");
579 >    exit(1);
580 >  }
581 >  DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
582  
583    if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
584        getpath(clear_paths.scat, ".", R_OK) == NULL ||
585        getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
586        getpath(clear_paths.irrad, ".", R_OK) == NULL) {
587 <    printf("# Precomputing...\n");
587 >    printf("# Pre-computing...\n");
588      if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
589 <      fprintf(stderr, "Precompute failed\n");
589 >      fprintf(stderr, "Pre-compute failed\n");
590        return 0;
591      }
592    }
# Line 387 | Line 596 | int main(int argc, char *argv[]) {
596    DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
597    DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
598  
599 +  write_header(argc, argv, ccover, grefl, res);
600 +
601    if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
602 <                     irrad_clear_dp, ccover, sundir, grefl, res, outname)) {
602 >                     irrad_clear_dp, ccover, sundir, grefl, res, outname, ddir,
603 >                     dirnorm, difhor)) {
604      fprintf(stderr, "gen_spect_sky failed\n");
605      exit(1);
606    }

Diff Legend

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