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

Diff Legend

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