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.6 by greg, Wed Oct 9 17:22:42 2024 UTC vs.
Revision 2.9 by greg, Fri Jul 11 18:12:25 2025 UTC

# Line 1 | Line 1
1 – #include "color.h"
1   #ifndef lint
2 < static const char RCSid[] =
4 <    "$Id$";
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  
7   #include "atmos.h"
8   #include "copyright.h"
9 + #include "color.h"
10 + #include "paths.h"
11   #include "resolu.h"
12   #include "rtio.h"
13   #include <ctype.h>
# Line 19 | Line 19 | static const char RCSid[] =
19   #include <sys/types.h>
20   #endif
21  
22 – char *progname;
23 –
22   const double ARCTIC_LAT = 67.;
23   const double TROPIC_LAT = 23.;
24   const int SUMMER_START = 4;
# Line 31 | Line 29 | const double D65EFF = 203.; /* standard illuminant D65
29  
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};
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}};
37 > struct
38 > {
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) {
50 > static int
51 > make_directory
52 > (
53 >        const char *path
54 > )
55 > {
56   #ifdef _WIN32
57 <  if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
58 <    return 1;
59 <  }
60 <  return 0;
57 >        if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
58 >                return 1;
59 >        }
60 >        return 0;
61 >
62   #else
63 <  if (mkdir(path, 0777) == 0 || errno == EEXIST) {
64 <    return 1;
65 <  }
66 <  return 0;
63 >        if (mkdir(path, 0777) == 0 || errno == EEXIST) {
64 >                return 1;
65 >        }
66 >        return 0;
67 >
68   #endif
69   }
70  
71 < inline static float deg2rad(float deg) { return deg * (PI / 180.); }
71 > inline static float
72 > deg2rad
73 > (
74 >        float deg
75 > )
76 > {
77 >        return deg * (PI / 180.);
78 > }
79  
80 < static int cvthour(char *hs, int *tsolar, double *hour) {
81 <  char *cp = hs;
82 <  int i, j;
80 > static int
81 > cvthour
82 > (
83 >        char *hs,
84 >        int *tsolar,
85 >        double *hour
86 > )
87 > {
88 >        char *cp = hs;
89 >        int i, j;
90  
91 <  if ((*tsolar = *cp == '+'))
92 <    cp++; /* solar time? */
93 <  while (isdigit(*cp))
94 <    cp++;
95 <  if (*cp == ':')
96 <    *hour = atoi(hs) + atoi(++cp) / 60.0;
97 <  else {
98 <    *hour = atof(hs);
99 <    if (*cp == '.')
100 <      cp++;
101 <  }
102 <  while (isdigit(*cp))
103 <    cp++;
104 <  if (!*cp)
105 <    return (0);
106 <  if (*tsolar || !isalpha(*cp)) {
107 <    fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
108 <    exit(1);
109 <  }
110 <  i = 0;
111 <  do {
112 <    for (j = 0; cp[j]; j++)
113 <      if (toupper(cp[j]) != tzone[i].zname[j])
114 <        break;
115 <    if (!cp[j] && !tzone[i].zname[j]) {
116 <      s_meridian = tzone[i].zmer * (PI / 180);
117 <      return (1);
118 <    }
119 <  } while (tzone[i++].zname[0]);
91 >        if ((*tsolar = *cp == '+')) {
92 >                cp++;                                 /* solar time? */
93 >        }
94 >        while (isdigit(*cp)) {
95 >                cp++;
96 >        }
97 >        if (*cp == ':') {
98 >                *hour = atoi(hs) + atoi(++cp) / 60.0;
99 >        }else{
100 >                *hour = atof(hs);
101 >                if (*cp == '.') {
102 >                        cp++;
103 >                }
104 >        }
105 >        while (isdigit(*cp)) {
106 >                cp++;
107 >        }
108 >        if (!*cp) {
109 >                return (0);
110 >        }
111 >        if (*tsolar || !isalpha(*cp)) {
112 >                fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
113 >                exit(1);
114 >        }
115 >        i = 0;
116 >        do {
117 >                for (j = 0; cp[j]; j++) {
118 >                        if (toupper(cp[j]) != tzone[i].zname[j]) {
119 >                                break;
120 >                        }
121 >                }
122 >                if (!cp[j] && !tzone[i].zname[j]) {
123 >                        s_meridian = tzone[i].zmer * (PI / 180);
124 >                        return (1);
125 >                }
126 >        } while (tzone[i++].zname[0]);
127  
128 <  fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
129 <  fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
130 <  for (i = 1; tzone[i].zname[0]; i++)
131 <    fprintf(stderr, " %s", tzone[i].zname);
132 <  putc('\n', stderr);
133 <  exit(1);
128 >        fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
129 >        fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
130 >        for (i = 1; tzone[i].zname[0]; i++) {
131 >                fprintf(stderr, " %s", tzone[i].zname);
132 >        }
133 >        putc('\n', stderr);
134 >        exit(1);
135   }
136  
137 < static void basename(const char *path, char *output, size_t outsize) {
138 <  const char *last_slash = strrchr(path, '/');
139 <  const char *last_backslash = strrchr(path, '\\');
140 <  const char *filename = path;
141 <  const char *last_dot;
137 > static void
138 > basename
139 > (
140 >        const char *path,
141 >        char *output,
142 >        size_t outsize
143 > )
144 > {
145 >        const char *last_slash = strrchr(path, '/');
146 >        const char *last_backslash = strrchr(path, '\\');
147 >        const char *filename = path;
148 >        const char *last_dot;
149  
150 <  if (last_slash && last_backslash) {
151 <    filename =
152 <        (last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
153 <  } else if (last_slash) {
154 <    filename = last_slash + 1;
155 <  } else if (last_backslash) {
156 <    filename = last_backslash + 1;
157 <  }
150 >        if (last_slash && last_backslash) {
151 >                filename =
152 >                        (last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
153 >        } else if (last_slash) {
154 >                filename = last_slash + 1;
155 >        } else if (last_backslash) {
156 >                filename = last_backslash + 1;
157 >        }
158  
159 <  last_dot = strrchr(filename, '.');
160 <  if (last_dot) {
161 <    size_t length = last_dot - filename;
162 <    if (length < outsize) {
163 <      strncpy(output, filename, length);
164 <      output[length] = '\0';
165 <    } else {
166 <      strncpy(output, filename, outsize - 1);
167 <      output[outsize - 1] = '\0';
168 <    }
169 <  }
159 >        last_dot = strrchr(filename, '.');
160 >        if (last_dot) {
161 >                size_t length = last_dot - filename;
162 >                if (length < outsize) {
163 >                        strncpy(output, filename, length);
164 >                        output[length] = '\0';
165 >                } else {
166 >                        strncpy(output, filename, outsize - 1);
167 >                        output[outsize - 1] = '\0';
168 >                }
169 >        }
170   }
171  
172 < static char *join_paths(const char *path1, const char *path2) {
173 <  size_t len1 = strlen(path1);
174 <  size_t len2 = strlen(path2);
175 <  int need_separator = (path1[len1 - 1] != DIRSEP);
172 > static char *
173 > join_paths
174 > (
175 >        const char *path1,
176 >        const char *path2
177 > )
178 > {
179 >        size_t len1 = strlen(path1);
180 >        size_t len2 = strlen(path2);
181 >        int need_separator = (path1[len1 - 1] != DIRSEP);
182  
183 <  char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
184 <  if (!result)
185 <    return NULL;
183 >        char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
184 >        if (!result) {
185 >                return NULL;
186 >        }
187  
188 <  strcpy(result, path1);
189 <  if (need_separator) {
190 <    result[len1] = DIRSEP;
191 <    len1++;
192 <  }
193 <  strcpy(result + len1, path2);
188 >        strcpy(result, path1);
189 >        if (need_separator) {
190 >                result[len1] = DIRSEP;
191 >                len1++;
192 >        }
193 >        strcpy(result + len1, path2);
194  
195 <  return result;
195 >        return result;
196   }
197  
198 < static inline double wmean2(const double a, const double b, const double x) {
199 <  return a * (1 - x) + b * x;
198 > static inline double
199 > wmean2
200 > (
201 >        const double a,
202 >        const double b,
203 >        const double x
204 > )
205 > {
206 >        return a * (1 - x) + b * x;
207   }
208  
209 < static inline double wmean(const double a, const double x, const double b,
210 <                           const double y) {
211 <  return (a * x + b * y) / (a + b);
209 > static inline double
210 > wmean
211 > (
212 >        const double a,
213 >        const double x,
214 >        const double b,
215 >        const double y
216 > )
217 > {
218 >        return (a * x + b * y) / (a + b);
219   }
220  
221 < static double get_overcast_zenith_brightness(const double sundir[3]) {
222 <  double zenithbr;
223 <  if (sundir[2] < 0) {
224 <    zenithbr = 0;
225 <  } else {
226 <    zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
227 <  }
228 <  return zenithbr;
221 > static double
222 > get_overcast_zenith_brightness
223 > (
224 >        const double sundir[3]
225 > )
226 > {
227 >        double zenithbr;
228 >        if (sundir[2] < 0) {
229 >                zenithbr = 0;
230 >        } else {
231 >                zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
232 >        }
233 >        return zenithbr;
234   }
235  
236   /* from gensky.c */
237 < static double get_overcast_brightness(const double dz, const double zenithbr) {
238 <  double groundbr = zenithbr * GNORM;
239 <  return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
240 <               pow(dz + 1.01, -10), groundbr);
237 > static double
238 > get_overcast_brightness
239 > (
240 >        const double dz,
241 >        const double zenithbr
242 > )
243 > {
244 >        double groundbr = zenithbr * GNORM;
245 >        return wmean(
246 >                pow(dz + 1.01, 10),
247 >                zenithbr * (1 + 2 * dz) / 3,
248 >                pow(dz + 1.01, -10),
249 >                groundbr);
250   }
251  
252 < static void write_header(const int argc, char **argv, const double cloud_cover,
253 <                         const double grefl, const int res) {
254 <  int i;
255 <  printf("# ");
256 <  for (i = 0; i < argc; i++) {
257 <    printf("%s ", argv[i]);
258 <  }
259 <  printf("\n");
260 <  printf(
261 <      "#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
262 <      cloud_cover, grefl, res);
252 > static void
253 > write_header
254 > (
255 >        const int argc,
256 >        char **argv,
257 >        const double cloud_cover,
258 >        const double grefl,
259 >        const int res
260 > )
261 > {
262 >        int i;
263 >        printf("# ");
264 >        for (i = 0; i < argc; i++) {
265 >                printf("%s ", argv[i]);
266 >        }
267 >        printf("\n");
268 >        printf(
269 >                "#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: "
270 >                "%d\n\n",
271 >                cloud_cover,
272 >                grefl,
273 >                res);
274   }
275  
276 < static void write_rad(const double *sun_radiance, const double intensity,
277 <                      const FVECT sundir, const char *ddir,
278 <                      const char *skyfile) {
279 <  if (sundir[2] > 0) {
280 <    printf("void spectrum sunrad\n0\n0\n22 380 780 ");
281 <    int i;
282 <    for (i = 0; i < NSSAMP; ++i) {
283 <      printf("%.3f ", sun_radiance[i]);
284 <    }
285 <    printf("\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n", intensity,
286 <           intensity, intensity);
287 <    printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
288 <           sundir[2]);
289 <  }
290 <  printf("void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
291 <         "'1-Acos(Dz)/PI'\n0\n0\n\n",
292 <         skyfile);
276 > static void
277 > write_rad
278 > (
279 >        const double *sun_radiance,
280 >        const double intensity,
281 >        const FVECT sundir,
282 >        const char *ddir,
283 >        const char *skyfile
284 > )
285 > {
286 >        if (sundir[2] > 0) {
287 >                printf("void spectrum sunrad\n0\n0\n22 380 780 ");
288 >                int i;
289 >                for (i = 0; i < NSSAMP; ++i) {
290 >                        printf("%.3f ", sun_radiance[i]);
291 >                }
292 >                printf(
293 >                        "\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n",
294 >                        intensity,
295 >                        intensity,
296 >                        intensity);
297 >                printf(
298 >                        "solar source sun\n0\n0\n4 %f %f %f 0.533\n\n",
299 >                        sundir[0],
300 >                        sundir[1],
301 >                        sundir[2]);
302 >        }
303 >        printf(
304 >                "void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
305 >                "'1-Acos(Dz)/PI'\n0\n0\n\n",
306 >                skyfile);
307   }
308  
309 < static void write_hsr_header(FILE *fp, RESOLU *res) {
310 <  float wvsplit[4] = {380, 480, 588, 780};
311 <  newheader("RADIANCE", fp);
312 <  fputncomp(NSSAMP, fp);
313 <  fputwlsplit(wvsplit, fp);
314 <  fputformat(SPECFMT, fp);
315 <  fputc('\n', fp);
316 <  fputsresolu(res, fp);
309 > static void
310 > write_hsr_header
311 > (
312 >        FILE *fp,
313 >        RESOLU *res
314 > )
315 > {
316 >        newheader("RADIANCE", fp);
317 >        fputncomp(NSSAMP, fp);
318 >        fputwlsplit(WLPART, fp);
319 >        fputformat(SPECFMT, fp);
320 >        fputc('\n', fp);
321 >        fputsresolu(res, fp);
322   }
323  
324 < static inline float frac(float x) { return x - floor(x); }
324 > static void
325 > reverse_array_float
326 > (
327 >        float arr[],
328 >        int size
329 > )
330 > {
331 >        int start = 0;
332 >        int end = size - 1;
333  
334 < int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
335 <                  DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
336 <                  const double cloud_cover, const FVECT sundir,
337 <                  const double grefl, const int res, const char *outname,
338 <                  const char *ddir, const double dirnorm, const double difhor) {
339 <  char skyfile[PATH_MAX];
340 <  if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP,
341 <                outname)) {
234 <    fprintf(stderr, "Error setting sky file name\n");
235 <    return 0;
236 <  };
237 <  int xres = res;
238 <  int yres = xres / 2;
239 <  RESOLU rs = {PIXSTANDARD, xres, yres};
240 <  FILE *skyfp = fopen(skyfile, "w");
241 <  write_hsr_header(skyfp, &rs);
334 >        while (start < end) {
335 >                float temp = arr[start];
336 >                arr[start] = arr[end];
337 >                arr[end] = temp;
338 >                start++;
339 >                end--;
340 >        }
341 > }
342  
343 <  CNDX[3] = NSSAMP;
343 > int
344 > gen_spect_sky
345 > (
346 >        DATARRAY *tau_clear,
347 >        DATARRAY *scat_clear,
348 >        DATARRAY *scat1m_clear,
349 >        DATARRAY *irrad_clear,
350 >        const double cloud_cover,
351 >        const FVECT sundir,
352 >        const double grefl,
353 >        const int res,
354 >        const char *outname,
355 >        const char *ddir,
356 >        const double dirnorm,
357 >        const double difhor
358 > )
359 > {
360 >        char skyfile[PATH_MAX];
361 >        if (!snprintf(
362 >                        skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP, outname)) {
363 >                fprintf(stderr, "Error setting sky file name\n");
364 >                return 0;
365 >        }
366 >        ;
367 >        int xres = res;
368 >        int yres = xres / 2;
369 >        RESOLU rs = {PIXSTANDARD, xres, yres};
370 >        FILE *skyfp = fopen(skyfile, "w");
371 >        write_hsr_header(skyfp, &rs);
372  
373 <  FVECT view_point = {0, 0, ER + 10};
246 <  const double radius = VLEN(view_point);
247 <  const double sun_ct = fdot(view_point, sundir) / radius;
373 >        CNDX[3] = NSSAMP;
374  
375 <  double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
376 <  double overcast_grndbr = overcast_zenithbr * GNORM;
375 >        FVECT view_point = {0, 0, ER + 10};
376 >        const double radius = VLEN(view_point);
377 >        const double sun_ct = fdot(view_point, sundir) / radius;
378  
379 <  double dif_ratio = 1;
380 <  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};
379 >        double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
380 >        double overcast_grndbr = overcast_zenithbr * GNORM;
381  
382 <      float px = i / (xres - 1.0);
383 <      float py = j / (yres - 1.0);
384 <      float lambda = ((1 - py) * PI) - (PI / 2.0);
385 <      float phi = (px * 2.0 * PI) - PI;
382 >        double dif_ratio = 1;
383 >        if (difhor > 0) {
384 >                DATARRAY *indirect_irradiance_clear =
385 >                        get_indirect_irradiance(irrad_clear, radius, sun_ct);
386 >                double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
387 >                double diffuse_irradiance = 0;
388 >                int l;
389 >                for (l = 0; l < NSSAMP; ++l) {
390 >                        diffuse_irradiance +=
391 >                                indirect_irradiance_clear->arr.d[l] * 20;                                                                 /* 20nm interval */
392 >                }
393 >                free(indirect_irradiance_clear);
394 >                diffuse_irradiance =
395 >                        wmean2(diffuse_irradiance, overcast_ghi, cloud_cover);
396 >                if (diffuse_irradiance > 0) {
397 >                        dif_ratio =
398 >                                difhor / WHTEFFICACY / diffuse_irradiance / 1.15;                                                                 /* fudge */
399 >                }
400 >        }
401 >        int i, j, k;
402 >        for (j = 0; j < yres; ++j) {
403 >                for (i = 0; i < xres; ++i) {
404 >                        SCOLOR radiance = {0};
405 >                        SCOLR sky_sclr = {0};
406  
407 <      FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
408 <                    sin(lambda)};
407 >                        float px = i / (xres - 1.0);
408 >                        float py = j / (yres - 1.0);
409 >                        float lambda = ((1 - py) * PI) - (PI / 2.0);
410 >                        float phi = (px * 2.0 * PI) - PI;
411  
412 <      const double mu = fdot(view_point, rdir) / radius;
413 <      const double nu = fdot(rdir, sundir);
412 >                        FVECT rdir = {
413 >                                cos(lambda) * cos(phi), cos(lambda) * sin(phi), sin(lambda)
414 >                        };
415  
416 <      /* hit ground */
417 <      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 <      }
416 >                        const double mu = fdot(view_point, rdir) / radius;
417 >                        const double nu = fdot(rdir, sundir);
418  
419 <      for (k = 0; k < NSSAMP; ++k) {
420 <        radiance[k] *= WVLSPAN;
421 <      }
419 >                        /* hit ground */
420 >                        if (rdir[2] < 0) {
421 >                                get_ground_radiance(
422 >                                        tau_clear,
423 >                                        scat_clear,
424 >                                        scat1m_clear,
425 >                                        irrad_clear,
426 >                                        view_point,
427 >                                        rdir,
428 >                                        radius,
429 >                                        mu,
430 >                                        sun_ct,
431 >                                        nu,
432 >                                        grefl,
433 >                                        sundir,
434 >                                        radiance);
435 >                        } else {
436 >                                get_sky_radiance(
437 >                                        scat_clear, scat1m_clear, radius, mu, sun_ct, nu, radiance);
438 >                        }
439  
440 <      if (cloud_cover > 0) {
441 <        double skybr = get_overcast_brightness(rdir[2], overcast_zenithbr);
442 <        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 <      }
440 >                        for (k = 0; k < NSSAMP; ++k) {
441 >                                radiance[k] *= WVLSPAN;
442 >                        }
443  
444 <      for (k = 0; k < NSSAMP; ++k) {
445 <        radiance[k] *= dif_ratio;
446 <      }
444 >                        if (cloud_cover > 0) {
445 >                                double skybr =
446 >                                        get_overcast_brightness(rdir[2], overcast_zenithbr);
447 >                                if (rdir[2] < 0) {
448 >                                        for (k = 0; k < NSSAMP; ++k) {
449 >                                                radiance[k] = wmean2(
450 >                                                        radiance[k],
451 >                                                        overcast_grndbr * D6415[k],
452 >                                                        cloud_cover);
453 >                                        }
454 >                                } else {
455 >                                        for (k = 0; k < NSSAMP; ++k) {
456 >                                                radiance[k] =
457 >                                                        wmean2(radiance[k], skybr * D6415[k], cloud_cover);
458 >                                        }
459 >                                }
460 >                        }
461  
462 <      scolor2scolr(sky_sclr, radiance, NSSAMP);
463 <      putbinary(sky_sclr, LSCOLR, 1, skyfp);
464 <    }
316 <  }
317 <  fclose(skyfp);
462 >                        for (k = 0; k < NSSAMP; ++k) {
463 >                                radiance[k] *= dif_ratio;
464 >                        }
465  
466 <  /* 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 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 <  }
466 >                        reverse_array_float(radiance, NSSAMP);
467  
468 <  /* Normalize */
469 <  double sum = 0.0;
470 <  for (i = 0; i < NSSAMP; ++i) {
471 <    sum += sun_radiance[i];
472 <  }
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 <  }
468 >                        scolor2scolr(sky_sclr, radiance, NSSAMP);
469 >                        putbinary(sky_sclr, LSCOLR, 1, skyfp);
470 >                }
471 >        }
472 >        fclose(skyfp);
473  
474 <  write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
475 <  return 1;
474 >        /* Get solar radiance */
475 >        double sun_radiance[NSSAMP] = {0};
476 >        get_solar_radiance(
477 >                tau_clear,
478 >                scat_clear,
479 >                scat1m_clear,
480 >                sundir,
481 >                radius,
482 >                sun_ct,
483 >                sun_radiance);
484 >        if (cloud_cover > 0) {
485 >                double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
486 >                int i;
487 >                for (i = 0; i < NSSAMP; ++i) {
488 >                        sun_radiance[i] = wmean2(
489 >                                sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
490 >                }
491 >        }
492 >
493 >        /* Normalize */
494 >        double sum = 0.0;
495 >        for (i = 0; i < NSSAMP; ++i) {
496 >                sum += sun_radiance[i];
497 >        }
498 >        double mean = sum / NSSAMP;
499 >        for (i = 0; i < NSSAMP; ++i) {
500 >                sun_radiance[i] /= mean;
501 >        }
502 >        double intensity = mean * WVLSPAN;
503 >        if (dirnorm > 0) {
504 >                intensity = dirnorm / SOLOMG / WHTEFFICACY;
505 >        }
506 >
507 >        write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
508 >        return 1;
509   }
510  
511 < static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
512 <                           const char *tag) {
513 <  DpPaths paths;
511 > static DpPaths
512 > get_dppaths
513 > (
514 >        const char *dir,
515 >        const double aod,
516 >        const char *mname,
517 >        const char *tag
518 > )
519 > {
520 >        DpPaths paths;
521  
522 <  snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
523 <           mname, aod);
524 <  snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
525 <           mname, aod);
526 <  snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
527 <           tag, mname, aod);
528 <  snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
529 <           mname, aod);
522 >        snprintf(
523 >                paths.tau,
524 >                PATH_MAX,
525 >                "%s%ctau_%s_%s_%.2f.dat",
526 >                dir,
527 >                DIRSEP,
528 >                tag,
529 >                mname,
530 >                aod);
531 >        snprintf(
532 >                paths.scat,
533 >                PATH_MAX,
534 >                "%s%cscat_%s_%s_%.2f.dat",
535 >                dir,
536 >                DIRSEP,
537 >                tag,
538 >                mname,
539 >                aod);
540 >        snprintf(
541 >                paths.scat1m,
542 >                PATH_MAX,
543 >                "%s%cscat1m_%s_%s_%.2f.dat",
544 >                dir,
545 >                DIRSEP,
546 >                tag,
547 >                mname,
548 >                aod);
549 >        snprintf(
550 >                paths.irrad,
551 >                PATH_MAX,
552 >                "%s%cirrad_%s_%s_%.2f.dat",
553 >                dir,
554 >                DIRSEP,
555 >                tag,
556 >                mname,
557 >                aod);
558  
559 <  return paths;
559 >        return paths;
560   }
561  
562 < static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
563 <                                         const int is_summer,
564 <                                         const double s_latitude) {
565 <  if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
566 <    tag[0] = 's';
567 <    if (is_summer) {
568 <      tag[1] = 's';
569 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
570 <      atmos->beta_r0 = BR0_SS;
571 <    } else {
572 <      tag[1] = 'w';
573 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
574 <      atmos->beta_r0 = BR0_SW;
575 <    }
576 <  } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
577 <    tag[0] = 'm';
578 <    if (is_summer) {
579 <      tag[1] = 's';
580 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
581 <      atmos->beta_r0 = BR0_MS;
582 <    } else {
583 <      tag[1] = 'w';
584 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
585 <      atmos->beta_r0 = BR0_MW;
586 <    }
587 <  } else {
588 <    tag[0] = 't';
589 <    tag[1] = 'r';
590 <    atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
591 <    atmos->beta_r0 = BR0_T;
592 <  }
593 <  tag[2] = '\0';
562 > static void
563 > set_rayleigh_density_profile
564 > (
565 >        Atmosphere *atmos,
566 >        char *tag,
567 >        const int is_summer,
568 >        const double s_latitude
569 > )
570 > {
571 >        if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
572 >                tag[0] = 's';
573 >                if (is_summer) {
574 >                        tag[1] = 's';
575 >                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
576 >                        atmos->beta_r0 = BR0_SS;
577 >                } else {
578 >                        tag[1] = 'w';
579 >                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
580 >                        atmos->beta_r0 = BR0_SW;
581 >                }
582 >        } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
583 >                tag[0] = 'm';
584 >                if (is_summer) {
585 >                        tag[1] = 's';
586 >                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
587 >                        atmos->beta_r0 = BR0_MS;
588 >                } else {
589 >                        tag[1] = 'w';
590 >                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
591 >                        atmos->beta_r0 = BR0_MW;
592 >                }
593 >        } else {
594 >                tag[0] = 't';
595 >                tag[1] = 'r';
596 >                atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
597 >                atmos->beta_r0 = BR0_T;
598 >        }
599 >        tag[2] = '\0';
600   }
601  
602 < static Atmosphere init_atmos(const double aod, const double grefl) {
603 <  Atmosphere atmos = {.ozone_density = {.layers =
604 <                                            {
605 <                                                {.width = 25000.0,
606 <                                                 .exp_term = 0.0,
607 <                                                 .exp_scale = 0.0,
608 <                                                 .linear_term = 1.0 / 15000.0,
609 <                                                 .constant_term = -2.0 / 3.0},
610 <                                                {.width = AH,
611 <                                                 .exp_term = 0.0,
612 <                                                 .exp_scale = 0.0,
613 <                                                 .linear_term = -1.0 / 15000.0,
614 <                                                 .constant_term = 8.0 / 3.0},
615 <                                            }},
616 <                      .rayleigh_density = {.layers =
617 <                                               {
618 <                                                   {.width = AH,
619 <                                                    .exp_term = 1.0,
620 <                                                    .exp_scale = -1.0 / HR_MS,
621 <                                                    .linear_term = 0.0,
622 <                                                    .constant_term = 0.0},
623 <                                               }},
624 <                      .beta_r0 = BR0_MS,
625 <                      .beta_scale = aod / AOD0_CA,
626 <                      .beta_m = NULL,
627 <                      .grefl = grefl};
628 <  return atmos;
602 > static Atmosphere
603 > init_atmos
604 > (
605 >        const double aod,
606 >        const double grefl
607 > )
608 > {
609 >        Atmosphere atmos = {
610 >                .ozone_density =
611 >                {.layers =
612 >                 {
613 >                         {.width = 25000.0,
614 >                                          .exp_term = 0.0,
615 >                                          .exp_scale = 0.0,
616 >                                          .linear_term = 1.0 / 15000.0,
617 >                                          .constant_term = -2.0 / 3.0},
618 >                         {.width = AH,
619 >                                          .exp_term = 0.0,
620 >                                          .exp_scale = 0.0,
621 >                                          .linear_term = -1.0 / 15000.0,
622 >                                          .constant_term = 8.0 / 3.0},
623 >                 }},
624 >                .rayleigh_density =
625 >                {.layers =
626 >                 {
627 >                         {.width = AH,
628 >                                          .exp_term = 1.0,
629 >                                          .exp_scale = -1.0 / HR_MS,
630 >                                          .linear_term = 0.0,
631 >                                          .constant_term = 0.0},
632 >                 }},
633 >                .beta_r0 = BR0_MS,
634 >                .beta_scale = aod / AOD0_CA,
635 >                .beta_m = NULL,
636 >                .grefl = grefl
637 >        };
638 >        return atmos;
639   }
640  
641 < int main(int argc, char *argv[]) {
642 <  progname = argv[0];
643 <  int month, day;
644 <  double hour;
645 <  FVECT sundir;
646 <  int num_threads = 1;
647 <  int sorder = 4;
648 <  int year = 0;
649 <  int tsolar = 0;
650 <  int got_meridian = 0;
651 <  double grefl = 0.2;
652 <  double ccover = 0.0;
653 <  int res = 64;
654 <  double aod = AOD0_CA;
655 <  char *outname = "out";
656 <  char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
657 <  char mie_name[20] = "mie_ca";
658 <  char lstag[3];
659 <  char *ddir = ".";
660 <  int i;
661 <  double dirnorm = 0; /* direct normal illuminance */
662 <  double difhor = 0;  /* diffuse horizontal illuminance */
641 > int
642 > main
643 > (
644 >        int argc,
645 >        char *argv[]
646 > )
647 > {
648 >        int month, day;
649 >        double hour;
650 >        FVECT sundir;
651 >        int num_threads = 1;
652 >        int sorder = 4;
653 >        int year = 0;
654 >        int tsolar = 0;
655 >        int got_meridian = 0;
656 >        double grefl = 0.2;
657 >        double ccover = 0.0;
658 >        int res = 64;
659 >        double aod = AOD0_CA;
660 >        char *outname = "out";
661 >        char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
662 >        char mie_name[20] = "mie_ca";
663 >        char lstag[3];
664 >        char *ddir = ".";
665 >        int i;
666 >        double dirnorm = 0;                 /* direct normal illuminance */
667 >        double difhor = 0;                 /* diffuse horizontal illuminance */
668  
669 <  if (argc == 2 && !strcmp(argv[1], "-defaults")) {
670 <    printf("-i %d\t\t\t\t#scattering order\n", sorder);
671 <    printf("-g %f\t\t\t#ground reflectance\n", grefl);
672 <    printf("-c %f\t\t\t#cloud cover\n", ccover);
673 <    printf("-r %d\t\t\t\t#image resolution\n", res);
674 <    printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
675 <    printf("-f %s\t\t\t\t#output name (-f)\n", outname);
676 <    printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
677 <    exit(0);
678 <  }
669 >        fixargv0(argv[0]);
670 >        if (argc == 2 && !strcmp(argv[1], "-defaults")) {
671 >                printf("-i %d\t\t\t\t#scattering order\n", sorder);
672 >                printf("-g %f\t\t\t#ground reflectance\n", grefl);
673 >                printf("-c %f\t\t\t#cloud cover\n", ccover);
674 >                printf("-r %d\t\t\t\t#image resolution\n", res);
675 >                printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
676 >                printf("-f %s\t\t\t\t#output name (-f)\n", outname);
677 >                printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
678 >                exit(0);
679 >        }
680  
681 <  if (argc < 4) {
682 <    fprintf(stderr,
683 <            "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r "
684 <            "res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
685 <            "-g grefl -f outpath\n",
686 <            argv[0]);
687 <    return 0;
688 <  }
681 >        if (argc < 4) {
682 >                fprintf(
683 >                        stderr,
684 >                        "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod "
685 >                        "-r "
686 >                        "res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
687 >                        "-g grefl -f outpath\n",
688 >                        argv[0]);
689 >                return 0;
690 >        }
691  
692 <  month = atoi(argv[1]);
693 <  if (month < 1 || month > 12) {
694 <    fprintf(stderr, "bad month");
695 <    exit(1);
696 <  }
697 <  day = atoi(argv[2]);
698 <  if (day < 1 || day > 31) {
699 <    fprintf(stderr, "bad month");
700 <    exit(1);
701 <  }
702 <  got_meridian = cvthour(argv[3], &tsolar, &hour);
692 >        month = atoi(argv[1]);
693 >        if (month < 1 || month > 12) {
694 >                fprintf(stderr, "bad month");
695 >                exit(1);
696 >        }
697 >        day = atoi(argv[2]);
698 >        if (day < 1 || day > 31) {
699 >                fprintf(stderr, "bad month");
700 >                exit(1);
701 >        }
702 >        got_meridian = cvthour(argv[3], &tsolar, &hour);
703  
704 <  if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
705 <    fprintf(stderr, "Cannot compute solar angle\n");
706 <    exit(1);
707 <  }
704 >        if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
705 >                fprintf(stderr, "Cannot compute solar angle\n");
706 >                exit(1);
707 >        }
708  
709 <  for (i = 4; i < argc; i++) {
710 <    if (argv[i][0] == '-') {
711 <      switch (argv[i][1]) {
712 <      case 'a':
713 <        s_latitude = atof(argv[++i]) * (PI / 180.0);
714 <        break;
715 <      case 'c':
716 <        ccover = atof(argv[++i]);
717 <        break;
718 <      case 'd':
719 <        aod = atof(argv[++i]);
720 <        break;
721 <      case 'f':
722 <        outname = argv[++i];
723 <        break;
724 <      case 'g':
725 <        grefl = atof(argv[++i]);
726 <        break;
727 <      case 'i':
728 <        sorder = atoi(argv[++i]);
729 <        break;
730 <      case 'l':
731 <        mie_path = argv[++i];
732 <        basename(mie_path, mie_name, sizeof(mie_name));
733 <        break;
734 <      case 'm':
735 <        if (got_meridian) {
736 <          ++i;
737 <          break;
738 <        }
739 <        s_meridian = atof(argv[++i]) * (PI / 180.0);
740 <        break;
741 <      case 'n':
742 <        num_threads = atoi(argv[++i]);
743 <        break;
744 <      case 'o':
745 <        s_longitude = atof(argv[++i]) * (PI / 180.0);
746 <        break;
747 <      case 'L':
748 <        dirnorm = atof(argv[++i]);
749 <        difhor = atof(argv[++i]);
750 <        break;
751 <      case 'p':
752 <        ddir = argv[++i];
753 <        break;
754 <      case 'r':
755 <        res = atoi(argv[++i]);
756 <        break;
757 <      case 'y':
758 <        year = atoi(argv[++i]);
759 <        break;
760 <      default:
761 <        fprintf(stderr, "Unknown option %s\n", argv[i]);
762 <        exit(1);
763 <      }
764 <    }
765 <  }
766 <  if (year && (year < 1950) | (year > 2050))
767 <    fprintf(stderr, "%s: warning - year should be in range 1950-2050\n",
768 <            progname);
769 <  if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180)
770 <    fprintf(stderr,
771 <            "%s: warning - %.1f hours btwn. standard meridian and longitude\n",
772 <            progname, (s_longitude - s_meridian) * 12 / PI);
709 >        for (i = 4; i < argc; i++) {
710 >                if (argv[i][0] == '-') {
711 >                        switch (argv[i][1]) {
712 >                        case 'a':
713 >                                s_latitude = atof(argv[++i]) * (PI / 180.0);
714 >                                break;
715 >                        case 'c':
716 >                                ccover = atof(argv[++i]);
717 >                                break;
718 >                        case 'd':
719 >                                aod = atof(argv[++i]);
720 >                                break;
721 >                        case 'f':
722 >                                outname = argv[++i];
723 >                                break;
724 >                        case 'g':
725 >                                grefl = atof(argv[++i]);
726 >                                break;
727 >                        case 'i':
728 >                                sorder = atoi(argv[++i]);
729 >                                break;
730 >                        case 'l':
731 >                                mie_path = argv[++i];
732 >                                basename(mie_path, mie_name, sizeof(mie_name));
733 >                                break;
734 >                        case 'm':
735 >                                if (got_meridian) {
736 >                                        ++i;
737 >                                        break;
738 >                                }
739 >                                s_meridian = atof(argv[++i]) * (PI / 180.0);
740 >                                break;
741 >                        case 'n':
742 >                                num_threads = atoi(argv[++i]);
743 >                                break;
744 >                        case 'o':
745 >                                s_longitude = atof(argv[++i]) * (PI / 180.0);
746 >                                break;
747 >                        case 'L':
748 >                                dirnorm = atof(argv[++i]);
749 >                                difhor = atof(argv[++i]);
750 >                                break;
751 >                        case 'p':
752 >                                ddir = argv[++i];
753 >                                break;
754 >                        case 'r':
755 >                                res = atoi(argv[++i]);
756 >                                break;
757 >                        case 'y':
758 >                                year = atoi(argv[++i]);
759 >                                break;
760 >                        default:
761 >                                fprintf(stderr, "Unknown option %s\n", argv[i]);
762 >                                exit(1);
763 >                        }
764 >                }
765 >        }
766 >        if (year && (year < 1950) | (year > 2050)) {
767 >                fprintf(
768 >                        stderr,
769 >                        "%s: warning - year should be in range 1950-2050\n",
770 >                        progname);
771 >        }
772 >        if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180) {
773 >                fprintf(
774 >                        stderr,
775 >                        "%s: warning - %.1f hours btwn. standard meridian and "
776 >                        "longitude\n",
777 >                        progname,
778 >                        (s_longitude - s_meridian) * 12 / PI);
779 >        }
780  
781 <  Atmosphere clear_atmos = init_atmos(aod, grefl);
781 >        Atmosphere clear_atmos = init_atmos(aod, grefl);
782  
783 <  int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
784 <  if (s_latitude < 0) {
785 <    is_summer = !is_summer;
786 <  }
787 <  set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
783 >        int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
784 >        if (s_latitude < 0) {
785 >                is_summer = !is_summer;
786 >        }
787 >        set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
788  
789 <  /* Load mie density data */
790 <  DATARRAY *mie_dp = getdata(mie_path);
791 <  if (mie_dp == NULL) {
792 <    fprintf(stderr, "Error reading mie data\n");
793 <    return 0;
794 <  }
795 <  clear_atmos.beta_m = mie_dp;
789 >        /* Load mie density data */
790 >        DATARRAY *mie_dp = getdata(mie_path);
791 >        if (mie_dp == NULL) {
792 >                fprintf(stderr, "Error reading mie data\n");
793 >                return 0;
794 >        }
795 >        clear_atmos.beta_m = mie_dp;
796  
797 <  char gsdir[PATH_MAX];
798 <  size_t siz = strlen(ddir);
799 <  if (ISDIRSEP(ddir[siz - 1]))
800 <    ddir[siz - 1] = '\0';
801 <  snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
802 <  if (!make_directory(gsdir)) {
803 <    fprintf(stderr, "Failed creating atmos_data directory");
804 <    exit(1);
805 <  }
806 <  DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
797 >        char gsdir[PATH_MAX];
798 >        size_t siz = strlen(ddir);
799 >        if (ISDIRSEP(ddir[siz - 1])) {
800 >                ddir[siz - 1] = '\0';
801 >        }
802 >        snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
803 >        if (!make_directory(gsdir)) {
804 >                fprintf(stderr, "Failed creating atmos_data directory");
805 >                exit(1);
806 >        }
807 >        DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
808  
809 <  if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
810 <      getpath(clear_paths.scat, ".", R_OK) == NULL ||
811 <      getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
812 <      getpath(clear_paths.irrad, ".", R_OK) == NULL) {
813 <    printf("# Pre-computing...\n");
814 <    if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
815 <      fprintf(stderr, "Pre-compute failed\n");
816 <      return 0;
817 <    }
818 <  }
809 >        if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
810 >                getpath(clear_paths.scat, ".", R_OK) == NULL ||
811 >                getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
812 >                getpath(clear_paths.irrad, ".", R_OK) == NULL) {
813 >                printf("# Pre-computing...\n");
814 >                if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
815 >                        fprintf(stderr, "Pre-compute failed\n");
816 >                        return 0;
817 >                }
818 >        }
819  
820 <  DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
821 <  DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
822 <  DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
823 <  DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
820 >        DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
821 >        DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
822 >        DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
823 >        DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
824  
825 <  write_header(argc, argv, ccover, grefl, res);
825 >        write_header(argc, argv, ccover, grefl, res);
826  
827 <  if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
828 <                     irrad_clear_dp, ccover, sundir, grefl, res, outname, ddir,
829 <                     dirnorm, difhor)) {
830 <    fprintf(stderr, "gen_spect_sky failed\n");
831 <    exit(1);
832 <  }
827 >        if (!gen_spect_sky(
828 >                        tau_clear_dp,
829 >                        scat_clear_dp,
830 >                        scat1m_clear_dp,
831 >                        irrad_clear_dp,
832 >                        ccover,
833 >                        sundir,
834 >                        grefl,
835 >                        res,
836 >                        outname,
837 >                        ddir,
838 >                        dirnorm,
839 >                        difhor)) {
840 >                fprintf(stderr, "gen_spect_sky failed\n");
841 >                exit(1);
842 >        }
843  
844 <  freedata(mie_dp);
845 <  freedata(tau_clear_dp);
846 <  freedata(scat_clear_dp);
847 <  freedata(irrad_clear_dp);
848 <  freedata(scat1m_clear_dp);
844 >        freedata(mie_dp);
845 >        freedata(tau_clear_dp);
846 >        freedata(scat_clear_dp);
847 >        freedata(irrad_clear_dp);
848 >        freedata(scat1m_clear_dp);
849  
850 <  return 1;
850 >        return 1;
851   }

Diff Legend

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