ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensdaymtx.c
(Generate patch)

Comparing ray/src/gen/gensdaymtx.c (file contents):
Revision 1.2 by greg, Fri Aug 2 19:01:13 2024 UTC vs.
Revision 1.5 by greg, Thu Apr 10 23:30:58 2025 UTC

# Line 1 | Line 1
1   #ifndef lint
2   static const char RCSid[] = "$Id$";
3   #endif
4 < #include "atmos.h"
5 < #include "copyright.h"
6 < #include "data.h"
7 < #include "platform.h"
8 < #include "rtio.h"
9 < #include <ctype.h>
4 >
5   #include <stdlib.h>
6 + #include <ctype.h>
7   #ifdef _WIN32
8   #include <windows.h>
9   #else
# Line 16 | Line 12 | static const char RCSid[] = "$Id$";
12   #include <sys/types.h>
13   #endif
14  
15 + #include "atmos.h"
16 + #include "copyright.h"
17 + #include "data.h"
18 + #include "platform.h"
19 + #include "rtio.h"
20 + #include "rtmath.h"
21 + #include "sun.h"
22 + #include "loadEPW.h"
23 +
24 + #ifndef M_PI
25 +    #define M_PI 3.14159265358579
26 + #endif
27 +
28 + #define vector(v, alt, azi)                                                    \
29 +  ((v)[1] = cos(alt), (v)[0] = (v)[1] * sin(azi), (v)[1] *= cos(azi),          \
30 +   (v)[2] = sin(alt))
31 +
32 + #define rh_vector(v, i) vector(v, rh_palt[i], rh_pazi[i])
33 +
34 + #define rh_cos(i) tsin(rh_palt[i])
35 +
36 + #define solar_minute(jd, hr) ((24 * 60) * ((jd) - 1) + (int)((hr) * 60. + .5))
37 +
38 +
39   char *progname;
40  
41   double altitude;   /* Solar altitude (radians) */
# Line 71 | Line 91 | float *rh_dom;  /* sky patch solid angle (sr) */
91  
92   double sun_ct;
93  
74 – #define vector(v, alt, azi)                                                    \
75 –  ((v)[1] = cos(alt), (v)[0] = (v)[1] * sin(azi), (v)[1] *= cos(azi),          \
76 –   (v)[2] = sin(alt))
77 –
78 – #define rh_vector(v, i) vector(v, rh_palt[i], rh_pazi[i])
79 –
80 – #define rh_cos(i) tsin(rh_palt[i])
81 –
82 – #define solar_minute(jd, hr) ((24 * 60) * ((jd) - 1) + (int)((hr) * 60. + .5))
83 –
94   inline void vectorize(double altitude, double azimuth, FVECT v) {
95    v[1] = cos(altitude);
96    v[0] = (v)[1] * sin(azimuth);
# Line 123 | Line 133 | static inline double wmean(const double a, const doubl
133    return (a * x + b * y) / (a + b);
134   }
135  
136 < static double get_zenith_brightness(const double sundir[3]) {
136 >
137 > static double get_overcast_zenith_brightness(const double sundir[3]) {
138    double zenithbr;
139    if (sundir[2] < 0) {
140      zenithbr = 0;
# Line 140 | Line 151 | static double get_overcast_brightness(const double dz,
151                 pow(dz + 1.01, -10), groundbr);
152   }
153  
154 + double
155 + solar_sunset(int month, int day)
156 + {
157 +        float W;
158 +        W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
159 +        return(12 + (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15));
160 + }
161 +
162 +
163 + double
164 + solar_sunrise(int month, int day)
165 + {
166 +        float W;
167 +        W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
168 +        return(12 - (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15));
169 + }
170 +
171   int rh_init(void) {
172   #define NROW 7
173    static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6};
# Line 196 | Line 224 | float *resize_dmatrix(float *mtx_data, int nsteps, int
224   }
225  
226   static Atmosphere init_atmos(const double aod, const double grefl) {
227 <  Atmosphere atmos = {.ozone_density = {.layers =
227 >    Atmosphere atmos = {.ozone_density = {.layers =
228                                              {
229                                                  {.width = 25000.0,
230                                                   .exp_term = 0.0,
# Line 221 | Line 249 | static Atmosphere init_atmos(const double aod, const d
249                        .beta_scale = aod / AOD0_CA,
250                        .beta_m = NULL,
251                        .grefl = grefl};
252 <  return atmos;
252 >    return atmos;
253   }
254  
255   static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
256                             const char *tag) {
257 <  DpPaths paths;
257 >    DpPaths paths;
258  
259 <  snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
259 >    snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
260             mname, aod);
261 <  snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
261 >    snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
262             mname, aod);
263 <  snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
263 >    snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
264             tag, mname, aod);
265 <  snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
265 >    snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
266             mname, aod);
267  
268 <  return paths;
268 >    return paths;
269   }
270 +
271   static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
272                                           const int is_summer,
273                                           const double s_latitude) {
274 <  /* Set rayleigh density profile */
275 <  if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
276 <    tag[0] = 's';
277 <    if (is_summer) {
278 <      tag[1] = 's';
279 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
280 <      atmos->beta_r0 = BR0_SS;
274 >    /* Set rayleigh density profile */
275 >    if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
276 >        tag[0] = 's';
277 >        if (is_summer) {
278 >            tag[1] = 's';
279 >            atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
280 >            atmos->beta_r0 = BR0_SS;
281 >        } else {
282 >            tag[1] = 'w';
283 >            atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
284 >            atmos->beta_r0 = BR0_SW;
285 >        }
286 >    } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
287 >        tag[0] = 'm';
288 >        if (is_summer) {
289 >            tag[1] = 's';
290 >            atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
291 >            atmos->beta_r0 = BR0_MS;
292 >        } else {
293 >            tag[1] = 'w';
294 >            atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
295 >            atmos->beta_r0 = BR0_MW;
296 >        }
297      } else {
298 <      tag[1] = 'w';
299 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
300 <      atmos->beta_r0 = BR0_SW;
298 >        tag[0] = 't';
299 >        tag[1] = 'r';
300 >        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
301 >        atmos->beta_r0 = BR0_T;
302      }
303 <  } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
258 <    tag[0] = 'm';
259 <    if (is_summer) {
260 <      tag[1] = 's';
261 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
262 <      atmos->beta_r0 = BR0_MS;
263 <    } else {
264 <      tag[1] = 'w';
265 <      atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
266 <      atmos->beta_r0 = BR0_MW;
267 <    }
268 <  } else {
269 <    tag[0] = 't';
270 <    tag[1] = 'r';
271 <    atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
272 <    atmos->beta_r0 = BR0_T;
273 <  }
274 <  tag[2] = '\0';
303 >    tag[2] = '\0';
304   }
305 +
306   /* Add in solar direct to nearest sky patches (GW) */
307   void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
308 <                DATARRAY *irrad, double ccover, float *parr) {
309 <  FVECT svec;
310 <  double near_dprod[NSUNPATCH];
311 <  int near_patch[NSUNPATCH];
312 <  double wta[NSUNPATCH], wtot;
313 <  int i, j, p;
308 >                DATARRAY *irrad, double ccover, double dirnorm, float *parr) {
309 >    FVECT svec;
310 >    double near_dprod[NSUNPATCH];
311 >    int near_patch[NSUNPATCH];
312 >    double wta[NSUNPATCH], wtot;
313 >    int i, j, p;
314  
315 <  /* identify nsuns closest patches */
316 <  for (i = nsuns; i--;)
317 <    near_dprod[i] = -1.;
318 <  vectorize(altitude, azimuth, svec);
319 <  for (p = 1; p < nskypatch; p++) {
320 <    FVECT pvec;
321 <    double dprod;
322 <    vectorize(rh_palt[p], rh_pazi[p], pvec);
323 <    dprod = DOT(pvec, svec);
324 <    for (i = 0; i < nsuns; i++)
325 <      if (dprod > near_dprod[i]) {
326 <        for (j = nsuns; --j > i;) {
327 <          near_dprod[j] = near_dprod[j - 1];
328 <          near_patch[j] = near_patch[j - 1];
315 >    /* identify nsuns closest patches */
316 >    for (i = nsuns; i--;)
317 >        near_dprod[i] = -1.;
318 >    vectorize(altitude, azimuth, svec);
319 >    for (p = 1; p < nskypatch; p++) {
320 >        FVECT pvec;
321 >        double dprod;
322 >        vectorize(rh_palt[p], rh_pazi[p], pvec);
323 >        dprod = DOT(pvec, svec);
324 >        for (i = 0; i < nsuns; i++)
325 >            if (dprod > near_dprod[i]) {
326 >                for (j = nsuns; --j > i;) {
327 >                    near_dprod[j] = near_dprod[j - 1];
328 >                    near_patch[j] = near_patch[j - 1];
329 >                }
330 >                near_dprod[i] = dprod;
331 >                near_patch[i] = p;
332 >                break;
333 >            }
334 >    }
335 >    /* Get solar radiance */
336 >    double sun_radiance[NSSAMP] = {0};
337 >    get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
338 >    if (ccover > 0) {
339 >        double zenithbr = get_overcast_zenith_brightness(sundir);
340 >        double skybr = get_overcast_brightness(sundir[2], zenithbr);
341 >        int l;
342 >        for (l = 0; l < NSSAMP; ++l) {
343 >            sun_radiance[l] = wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
344          }
300 –        near_dprod[i] = dprod;
301 –        near_patch[i] = p;
302 –        break;
303 –      }
304 –  }
305 –  /* Get solar radiance */
306 –  double sun_radiance[NSSAMP] = {0};
307 –  get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
308 –  if (ccover > 0) {
309 –    double zenithbr = get_zenith_brightness(sundir);
310 –    double skybr = get_overcast_brightness(sundir[2], zenithbr);
311 –    for (int l = 0; l < NSSAMP; ++l) {
312 –      sun_radiance[l] =
313 –          wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
345      }
346 <  }
347 <  /* weight by proximity */
348 <  wtot = 0;
349 <  for (i = nsuns; i--;)
319 <    wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
320 <  /* add to nearest patch radiances */
321 <  for (i = nsuns; i--;) {
322 <    float *pdest = parr + NSSAMP * near_patch[i];
323 <    for (int k = 0; k < NSSAMP; k++) {
324 <      *pdest++ = sun_radiance[k] * wta[i] / wtot;
346 >    /* Normalize */
347 >    double sum = 0.0;
348 >    for (i = 0; i < NSSAMP; ++i) {
349 >        sum += sun_radiance[i];
350      }
351 <  }
327 < }
351 >    double mean = sum / NSSAMP;
352  
353 < void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m, double ccover,
354 <                             float *parr) {
355 <  int i;
332 <  double mu_sky; /* Sun-sky point azimuthal angle */
333 <  double sspa;   /* Sun-sky point angle */
334 <  double zsa;    /* Zenithal sun angle */
335 <  FVECT view_point = {0, 0, ER};
336 <  for (i = 1; i < nskypatch; i++) {
337 <    FVECT rdir_sky;
338 <    vectorize(rh_palt[i], rh_pazi[i], rdir_sky);
339 <    mu_sky = fdot(view_point, rdir_sky) / ER;
340 <    sspa = fdot(rdir_sky, sundir);
341 <    SCOLOR sky_radiance = {0};
342 <
343 <    get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
344 <    for (int k = 0; k < NSSAMP; ++k) {
345 <      sky_radiance[k] *= WVLSPAN;
353 >    double intensity = mean * WVLSPAN;
354 >    if (dirnorm > 0) {
355 >        intensity = dirnorm / SOLOMG / WHTEFFICACY;
356      }
357 +    double dir_ratio = 1.;
358 +    if (mean > 0)
359 +        dir_ratio = intensity / mean;
360 +    for (i = 0; i < NSSAMP; ++i) {
361 +        sun_radiance[i] *= dir_ratio;
362 +    }
363  
364 <    if (ccover > 0) {
365 <      double zenithbr = get_zenith_brightness(sundir);
366 <      double grndbr = zenithbr * GNORM;
367 <      double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
368 <      for (int k = 0; k < NSSAMP; ++k) {
369 <        sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
370 <      }
364 >    /* weight by proximity */
365 >    wtot = 0;
366 >    for (i = nsuns; i--;)
367 >        wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
368 >    /* add to nearest patch radiances */
369 >    for (i = nsuns; i--;) {
370 >        float *pdest = parr + NSSAMP * near_patch[i];
371 >        int k;
372 >        for (k = 0; k < NSSAMP; k++) {
373 >            *pdest++ = sun_radiance[k] * wta[i] / wtot;
374 >        }
375      }
376 + }
377  
378 <    for (int k = 0; k < NSSAMP; ++k) {
379 <      parr[NSSAMP * i + k] = sky_radiance[k];
378 > void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m, DATARRAY *irrad_clear, double ccover, double dif_ratio,
379 >                             double overcast_zenithbr, float *parr) {
380 >    int i;
381 >    double mu_sky; /* Sun-sky point azimuthal angle */
382 >    double sspa;   /* Sun-sky point angle */
383 >    FVECT view_point = {0, 0, ER};
384 >    for (i = 1; i < nskypatch; i++) {
385 >        FVECT rdir_sky;
386 >        int k;
387 >        vectorize(rh_palt[i], rh_pazi[i], rdir_sky);
388 >        mu_sky = fdot(view_point, rdir_sky) / ER;
389 >        sspa = fdot(rdir_sky, sundir);
390 >        SCOLOR sky_radiance = {0};
391 >
392 >        get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
393 >        for (k = 0; k < NSSAMP; ++k) {
394 >            sky_radiance[k] *= WVLSPAN;
395 >        }
396 >
397 >        if (ccover > 0) {
398 >            double skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr);
399 >            int k;
400 >            for (k = 0; k < NSSAMP; ++k) {
401 >                sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
402 >            }
403 >        }
404 >
405 >        /* calibration */
406 >        for (k = 0; k < NSSAMP; ++k) {
407 >            sky_radiance[k] *= dif_ratio;
408 >        }
409 >
410 >        for (k = 0; k < NSSAMP; ++k) {
411 >            parr[NSSAMP * i + k] = sky_radiance[k];
412 >        }
413      }
360 –  }
414   }
415  
416   /* Return maximum of two doubles */
# Line 365 | Line 418 | static inline double dmax(double a, double b) { return
418  
419   /* Compute sky patch radiance values (modified by GW) */
420   void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
421 <                 DATARRAY *irrad, double ccover, float *parr) {
422 <  int index; /* Category index */
423 <  int i;
424 <  float sun_zenith;
425 <  SCOLOR sky_radiance = {0};
426 <  SCOLOR ground_radiance = {0};
427 <  SCOLR sky_sclr = {0};
428 <  SCOLR ground_sclr = {0};
429 <  FVECT view_point = {0, 0, ER};
430 <  const double radius = VLEN(view_point);
431 <  const double sun_ct = fdot(view_point, sundir) / radius;
432 <  const FVECT rdir_grnd = {0, 0, -1};
380 <  const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
381 <  const double nu_grnd = fdot(rdir_grnd, sundir);
421 >                 DATARRAY *irrad, double ccover, double difhor, float *parr) {
422 >    float sun_zenith;
423 >    SCOLOR sky_radiance = {0};
424 >    SCOLOR ground_radiance = {0};
425 >    SCOLR sky_sclr = {0};
426 >    SCOLR ground_sclr = {0};
427 >    FVECT view_point = {0, 0, ER};
428 >    const double radius = VLEN(view_point);
429 >    const double sun_ct = fdot(view_point, sundir) / radius;
430 >    const FVECT rdir_grnd = {0, 0, -1};
431 >    const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
432 >    const double nu_grnd = fdot(rdir_grnd, sundir);
433  
434 <  /* Calculate sun zenith angle (don't let it dip below horizon) */
435 <  /* Also limit minimum angle to keep circumsolar off zenith */
436 <  if (altitude <= 0.0)
437 <    sun_zenith = DegToRad(90.0);
438 <  else if (altitude >= DegToRad(87.0))
439 <    sun_zenith = DegToRad(3.0);
440 <  else
441 <    sun_zenith = DegToRad(90.0) - altitude;
434 >    /* Calculate sun zenith angle (don't let it dip below horizon) */
435 >    /* Also limit minimum angle to keep circumsolar off zenith */
436 >    if (altitude <= 0.0)
437 >        sun_zenith = DegToRad(90.0);
438 >    else if (altitude >= DegToRad(87.0))
439 >        sun_zenith = DegToRad(3.0);
440 >    else
441 >        sun_zenith = DegToRad(90.0) - altitude;
442  
443 <  /* Compute ground radiance (include solar contribution if any) */
444 <  get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
445 <                      mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
446 <  for (int j = 0; j < NSSAMP; j++) {
447 <    parr[j] *= WVLSPAN;
448 <  }
449 <  calc_sky_patch_radiance(scat, scat1m, ccover, parr);
443 >    double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
444 >
445 >    /* diffuse calibration factor */
446 >    double dif_ratio = 1;
447 >    if (difhor > 0) {
448 >        DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad, radius, sun_ct);
449 >        double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
450 >        double diffuse_irradiance = 0;
451 >        int l;
452 >        for (l = 0; l < NSSAMP; ++l) {
453 >            diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20;  /* 20nm interval */
454 >        }
455 >        /* free(indirect_irradiance_clear); */
456 >        diffuse_irradiance = wmean2(diffuse_irradiance, overcast_ghi, ccover);
457 >        if (diffuse_irradiance > 0) {
458 >            dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;       /* fudge */
459 >        }
460 >    }
461 >
462 >    /* Compute ground radiance (include solar contribution if any) */
463 >    get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
464 >                        mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
465 >    int j;
466 >    for (j = 0; j < NSSAMP; j++) {
467 >        parr[j] *= WVLSPAN;
468 >    }
469 >    calc_sky_patch_radiance(scat, scat1m, irrad, ccover, dif_ratio, overcast_zenithbr, parr);
470   }
471  
472   int main(int argc, char *argv[]) {
473  
474 <  char buf[256];
475 <  int doheader = 1; /* output header? */
476 <  double rotation = 0.0;
477 <  double elevation = 0;
478 <  int leap_day = 0;       /* add leap day? */
479 <  int sun_hours_only = 0; /* only output sun hours? */
480 <  float *mtx_data = NULL;
481 <  int ntsteps = 0;      /* number of time steps */
482 <  int tstorage = 0;     /* number of allocated time steps */
483 <  int nstored = 0;      /* number of time steps in matrix */
484 <  int last_monthly = 0; /* month of last report */
485 <  int mo, da;
486 <  double hr, aod, cc;
487 <  double dni, dhi;
488 <  int mtx_offset = 0;
489 <  int i, j;
490 <  char lstag[3];
491 <  char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
492 <  char *ddir = ".";
493 <  char mie_name[20] = "mie_ca";
494 <  int num_threads = 1;
495 <  int sorder = 4;
496 <  int solar_only = 0;
497 <  int sky_only = 0;
498 <  FVECT view_point = {0, 0, ER};
474 >    EPWheader   *epw = NULL;    /* EPW/WEA input file */
475 >    EPWrecord   erec;           /* current EPW/WEA input record */
476 >    int         doheader = 1; /* output header? */
477 >    double      rotation = 0.0;
478 >    double      elevation = 0;
479 >    int         leap_day = 0;       /* add leap day? */
480 >    int         sun_hours_only = 0; /* only output sun hours? */
481 >    int         dir_is_horiz;           /* direct is meas. on horizontal? */
482 >    float       *mtx_data = NULL;
483 >    int         ntsteps = 0;      /* number of time steps */
484 >    int         tstorage = 0;     /* number of allocated time steps */
485 >    int         nstored = 0;      /* number of time steps in matrix */
486 >    int         last_monthly = 0; /* month of last report */
487 >    double      dni, dhi;
488 >    int         mtx_offset = 0;
489 >        double      timeinterval = 0;
490 >    char        lstag[3];
491 >    char        *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
492 >    char        *ddir = ".";
493 >    char        mie_name[20] = "mie_ca";
494 >    int         num_threads = 1;
495 >    int         sorder = 4;
496 >    int         solar_only = 0;
497 >    int         sky_only = 0;
498 >    FVECT       view_point = {0, 0, ER};
499 >    int         i, j;
500  
501 <  progname = argv[0];
501 >    progname = argv[0];
502  
503 <  for (i = 1; i < argc && argv[i][0] == '-'; i++) {
504 <    switch (argv[i][1]) {
505 <    case 'd': /* solar (direct) only */
506 <      solar_only = 1;
507 <      break;
508 <    case 's': /* sky only (no direct) */
509 <      sky_only = 1;
510 <      break;
511 <    case 'g':
512 <      grefl = atof(argv[++i]);
513 <      break;
514 <    case 'm':
515 <      rhsubdiv = atoi(argv[++i]);
516 <      break;
517 <    case 'n':
518 <      num_threads = atoi(argv[++i]);
519 <      break;
520 <    case 'r': /* rotate distribution */
521 <      if (argv[i][2] && argv[i][2] != 'z')
503 >    for (i = 1; i < argc && argv[i][0] == '-'; i++) {
504 >        switch (argv[i][1]) {
505 >        case 'd': /* solar (direct) only */
506 >            solar_only = 1;
507 >            break;
508 >        case 's': /* sky only (no direct) */
509 >            sky_only = 1;
510 >            break;
511 >        case 'g':
512 >            grefl = atof(argv[++i]);
513 >            break;
514 >        case 'm':
515 >            rhsubdiv = atoi(argv[++i]);
516 >            break;
517 >        case 'n':
518 >            num_threads = atoi(argv[++i]);
519 >            break;
520 >        case 'r': /* rotate distribution */
521 >            if (argv[i][2] && argv[i][2] != 'z')
522 >                goto userr;
523 >            rotation = atof(argv[++i]);
524 >            break;
525 >        case 'u': /* solar hours only */
526 >            sun_hours_only = 1;
527 >            break;
528 >        case 'p':
529 >            ddir = argv[++i];
530 >            break;
531 >        case 'v': /* verbose progress reports */
532 >            verbose++;
533 >            break;
534 >        case 'h': /* turn off header */
535 >            doheader = 0;
536 >            break;
537 >        case '5': /* 5-phase calculation */
538 >            nsuns = 1;
539 >            fixed_sun_sa = PI / 360. * atof(argv[++i]);
540 >            if (fixed_sun_sa <= 0) {
541 >                fprintf(
542 >                    stderr,
543 >                    "%s: missing solar disk size argument for '-5' option\n",
544 >                    progname);
545 >                exit(1);
546 >            }
547 >            fixed_sun_sa *= fixed_sun_sa * PI;
548 >            break;
549 >        case 'i':
550 >            timeinterval = atof(argv[++i]);
551 >            break;
552 >        case 'o': /* output format */
553 >            switch (argv[i][2]) {
554 >            case 'f':
555 >            case 'd':
556 >            case 'a':
557 >                outfmt = argv[i][2];
558 >                break;
559 >            default:
560 >                goto userr;
561 >            }
562 >            break;
563 >        default:
564 >            goto userr;
565 >        }
566 >    }
567 >    if (i < argc - 1)
568          goto userr;
569 <      rotation = atof(argv[++i]);
570 <      break;
571 <    case 'u': /* solar hours only */
454 <      sun_hours_only = 1;
455 <      break;
456 <    case 'p':
457 <      ddir = argv[++i];
458 <      break;
459 <    case 'v': /* verbose progress reports */
460 <      verbose++;
461 <      break;
462 <    case 'h': /* turn off header */
463 <      doheader = 0;
464 <      break;
465 <    case '5': /* 5-phase calculation */
466 <      nsuns = 1;
467 <      fixed_sun_sa = PI / 360. * atof(argv[++i]);
468 <      if (fixed_sun_sa <= 0) {
469 <        fprintf(stderr,
470 <                "%s: missing solar disk size argument for '-5' option\n",
471 <                progname);
569 >
570 >    epw = EPWopen(argv[i]);
571 >    if (epw == NULL)
572          exit(1);
573 <      }
574 <      fixed_sun_sa *= fixed_sun_sa * PI;
575 <      break;
576 <    case 'o': /* output format */
577 <      switch (argv[i][2]) {
578 <      case 'f':
579 <      case 'd':
580 <      case 'a':
581 <        outfmt = argv[i][2];
573 >    if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
574 >        fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
575 >        exit(1);
576 >    }
577 >    if (verbose) {
578 >        if (i == argc - 1)
579 >            fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
580 >        else
581 >            fprintf(stderr, "%s: reading weather tape from <stdin>\n", progname);
582 >    }
583 >    s_latitude = epw->loc.latitude;
584 >    s_longitude = -epw->loc.longitude;
585 >    s_meridian = -15.*epw->loc.timezone;
586 >    elevation = epw->loc.elevation;
587 >    switch (epw->isWEA) {               /* translate units */
588 >    case WEAnot:
589 >    case WEAradnorm:
590 >        input = 1;              /* radiometric quantities */
591 >        dir_is_horiz = 0;    /* direct is perpendicular meas. */
592          break;
593 <      default:
594 <        goto userr;
595 <      }
596 <      break;
593 >    case WEAradhoriz:
594 >        input = 1;        /* radiometric quantities */
595 >        dir_is_horiz = 1;    /* solar measured horizontally */
596 >        break;
597 >    case WEAphotnorm:
598 >        input = 2;        /* photometric quantities */
599 >        dir_is_horiz = 0;    /* direct is perpendicular meas. */
600 >        break;
601      default:
602 <      goto userr;
602 >        goto fmterr;
603      }
490 –  }
491 –  if (i < argc - 1)
492 –    goto userr;
493 –  if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
494 –    fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
495 –    exit(1);
496 –  }
497 –  if (verbose) {
498 –    if (i == argc - 1)
499 –      fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
500 –    else
501 –      fprintf(stderr, "%s: reading weather tape from <stdin>\n", progname);
502 –  }
503 –  /* read weather tape header */
504 –  if (scanf("place %[^\r\n] ", buf) != 1)
505 –    goto fmterr;
506 –  if (scanf("latitude %lf\n", &s_latitude) != 1)
507 –    goto fmterr;
508 –  if (scanf("longitude %lf\n", &s_longitude) != 1)
509 –    goto fmterr;
510 –  if (scanf("time_zone %lf\n", &s_meridian) != 1)
511 –    goto fmterr;
512 –  if (scanf("site_elevation %lf\n", &elevation) != 1)
513 –    goto fmterr;
514 –  if (scanf("weather_data_file_units %d\n", &input) != 1)
515 –    goto fmterr;
604  
605 <  rh_init();
606 <  if (verbose) {
607 <    fprintf(stderr, "%s: location '%s'\n", progname, buf);
608 <    fprintf(stderr, "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
605 >    rh_init();
606 >
607 >    if (verbose) {
608 >        fprintf(stderr, "%s: location '%s %s'\n", progname, epw->loc.city, epw->loc.country);
609 >        fprintf(
610 >            stderr,
611 >            "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
612              progname, s_latitude, s_longitude);
613 <    if (rotation != 0)
614 <      fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
615 <  }
613 >        if (rotation != 0)
614 >            fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
615 >    }
616  
617 <  s_latitude = DegToRad(s_latitude);
618 <  s_longitude = DegToRad(s_longitude);
619 <  s_meridian = DegToRad(s_meridian);
620 <  /* initial allocation */
621 <  mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
617 >    s_latitude = DegToRad(s_latitude);
618 >    s_longitude = DegToRad(s_longitude);
619 >    s_meridian = DegToRad(s_meridian);
620 >    /* initial allocation */
621 >    mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
622  
623 <  /* Load mie density data */
624 <  DATARRAY *mie_dp = getdata(mie_path);
625 <  if (mie_dp == NULL) {
626 <    fprintf(stderr, "Error reading mie data\n");
627 <    return 0;
628 <  }
623 >    /* Load mie density data */
624 >    DATARRAY *mie_dp = getdata(mie_path);
625 >    if (mie_dp == NULL) {
626 >        fprintf(stderr, "Error reading mie data\n");
627 >        return 0;
628 >    }
629  
630 <  while (scanf("%d %d %lf %lf %lf %lf %lf\n", &mo, &da, &hr, &dni, &dhi, &aod,
631 <               &cc) == 7) {
632 <    double sda, sta;
633 <    int sun_in_sky;
634 <    /* compute solar position */
635 <    if ((mo == 2) & (da == 29)) {
545 <      julian_date = 60;
546 <      leap_day = 1;
547 <    } else
548 <      julian_date = jdate(mo, da) + leap_day;
549 <    sda = sdec(julian_date);
550 <    sta = stadj(julian_date);
551 <    altitude = salt(sda, hr + sta);
552 <    sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG / 2.));
630 >    if (epw->isWEA == WEAnot) {
631 >        fprintf(stderr, "EPW input\n");
632 >    } else if (epw->isWEA != WEAphotnorm) {
633 >        fprintf(stderr, "need WEA in photopic unit\n");
634 >        exit(1);
635 >    }
636  
637 <    azimuth = sazi(sda, hr + sta) + PI - DegToRad(rotation);
637 >    while ((j = EPWread(epw, &erec)) > 0) {
638 >        const int       mo = erec.date.month+1;
639 >        const int       da = erec.date.day;
640 >        const double    hr = erec.date.hour;
641 >        double aod = erec.optdepth;
642 >        double cc = erec.skycover;
643 >        double          sda, sta, st;
644 >        int             sun_in_sky;
645  
646 <    vectorize(altitude, azimuth, sundir);
647 <    if (sun_hours_only && sundir[2] <= 0.) {
648 <      continue; /* skipping nighttime points */
649 <    }
650 <    sun_ct = fdot(view_point, sundir) / ER;
646 >        if (aod == 0.0) {
647 >            aod = AOD0_CA;
648 >            fprintf(stderr, "aod is zero, using default value %.3f\n", AOD0_CA);
649 >        }
650 >        /* compute solar position */
651 >        if ((mo == 2) & (da == 29)) {
652 >            julian_date = 60;
653 >            leap_day = 1;
654 >        } else
655 >            julian_date = jdate(mo, da) + leap_day;
656 >        sda = sdec(julian_date);
657 >        sta = stadj(julian_date);
658 >        st = hr + sta;
659 >                if (timeinterval > 0) {
660 >                        if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120)
661 >                                st = (st + timeinterval/120 + solar_sunrise(mo, da))/2;
662 >                        else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120)
663 >                                st = (st - timeinterval/120 + solar_sunset(mo, da))/2;
664 >                }
665 >        altitude = salt(sda, st);
666 >        sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG / 2.));
667  
668 <    mtx_offset = NSSAMP * nskypatch * nstored;
563 <    nstored += 1;
564 <    /* make space for next row */
565 <    if (nstored > tstorage) {
566 <      tstorage += (tstorage >> 1) + nstored + 7;
567 <      mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
568 <    }
569 <    ntsteps++; /* keep count of time steps */
570 <               /* compute sky patch values */
571 <    Atmosphere clear_atmos = init_atmos(aod, grefl);
572 <    int is_summer = (mo >= SUMMER_START && mo <= SUMMER_END);
573 <    if (s_latitude < 0) {
574 <      is_summer = !is_summer;
575 <    }
576 <    set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
668 >        azimuth = sazi(sda, st) + PI - DegToRad(rotation);
669  
670 <    clear_atmos.beta_m = mie_dp;
670 >        vectorize(altitude, azimuth, sundir);
671 >        if (sun_hours_only && !sun_in_sky) {
672 >            continue; /* skipping nighttime points */
673 >        }
674 >        sun_ct = fdot(view_point, sundir) / ER;
675  
676 <    char gsdir[PATH_MAX];
677 <    size_t siz = strlen(ddir);
678 <    if (ISDIRSEP(ddir[siz - 1]))
583 <      ddir[siz - 1] = '\0';
584 <    snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
585 <    if (!make_directory(gsdir)) {
586 <      fprintf(stderr, "Failed creating atmos_data directory");
587 <      exit(1);
588 <    }
589 <    DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
676 >        dni = erec.dirillum;
677 >        dhi = erec.diffillum;
678 >        printf("%d %d %f %f %f %f %f\n", mo, da, hr, dni, dhi, aod, cc);
679  
680 <    if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
681 <        getpath(clear_paths.scat, ".", R_OK) == NULL ||
593 <        getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
594 <        getpath(clear_paths.irrad, ".", R_OK) == NULL) {
595 <      printf("# Pre-computing...\n");
596 <      if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
597 <        fprintf(stderr, "Pre-compute failed\n");
598 <        return 0;
599 <      }
600 <    }
680 >        mtx_offset = NSSAMP * nskypatch * nstored;
681 >        nstored += 1;
682  
683 <    DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
684 <    DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
685 <    DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
686 <    DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
683 >        /* make space for next row */
684 >        if (nstored > tstorage) {
685 >            tstorage += (tstorage >> 1) + nstored + 7;
686 >            mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
687 >        }
688 >        ntsteps++; /* keep count of time steps */
689  
690 <    if (!solar_only)
691 <      compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
692 <                  cc, mtx_data + mtx_offset);
693 <    if (!sky_only)
694 <      add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
695 <                 cc, mtx_data + mtx_offset);
696 <    /* update cumulative sky? */
614 <    for (i = NSSAMP * nskypatch * (ntsteps > 1); i--;)
615 <      mtx_data[i] += mtx_data[mtx_offset + i];
616 <    /* monthly reporting */
617 <    if (verbose && mo != last_monthly)
618 <      fprintf(stderr, "%s: stepping through month %d...\n", progname,
619 <              last_monthly = mo);
620 <    /* note whether leap-day was given */
690 >               /* compute sky patch values */
691 >        Atmosphere clear_atmos = init_atmos(aod, grefl);
692 >        int is_summer = (mo >= SUMMER_START && mo <= SUMMER_END);
693 >        if (s_latitude < 0) {
694 >            is_summer = !is_summer;
695 >        }
696 >        set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
697  
698 <    freedata(tau_clear_dp);
699 <    freedata(irrad_clear_dp);
700 <    freedata(scat_clear_dp);
701 <    freedata(scat1m_clear_dp);
702 <  }
703 <  freedata(mie_dp);
704 <  if (!ntsteps) {
705 <    fprintf(stderr, "%s: no valid time steps on input\n", progname);
706 <    exit(1);
707 <  }
708 <  /* check for junk at end */
709 <  while ((i = fgetc(stdin)) != EOF)
710 <    if (!isspace(i)) {
711 <      fprintf(stderr, "%s: warning - unexpected data past EOT: ", progname);
712 <      buf[0] = i;
713 <      buf[1] = '\0';
714 <      fgets(buf + 1, sizeof(buf) - 1, stdin);
715 <      fputs(buf, stderr);
716 <      fputc('\n', stderr);
717 <      break;
698 >        clear_atmos.beta_m = mie_dp;
699 >
700 >        char gsdir[PATH_MAX];
701 >        size_t siz = strlen(ddir);
702 >        if (ISDIRSEP(ddir[siz - 1]))
703 >            ddir[siz - 1] = '\0';
704 >            snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
705 >        if (!make_directory(gsdir)) {
706 >            fprintf(stderr, "Failed creating atmos_data directory");
707 >            exit(1);
708 >        }
709 >        DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
710 >
711 >        if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
712 >            getpath(clear_paths.scat, ".", R_OK) == NULL ||
713 >            getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
714 >            getpath(clear_paths.irrad, ".", R_OK) == NULL) {
715 >            fprintf(stderr, "# Pre-computing...\n");
716 >            if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
717 >                fprintf(stderr, "Pre-compute failed\n");
718 >                return 0;
719 >            }
720 >        }
721 >
722 >        DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
723 >        DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
724 >        DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
725 >        DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
726 >
727 >        if (!solar_only)
728 >            compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
729 >                  cc, dhi, mtx_data + mtx_offset);
730 >        if (!sky_only)
731 >            add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
732 >                 cc, dni, mtx_data + mtx_offset);
733 >        /* monthly reporting */
734 >        if (verbose && mo != last_monthly)
735 >            fprintf(stderr, "%s: stepping through month %d...\n", progname,
736 >                last_monthly = mo);
737      }
738 <  /* write out matrix */
739 <  if (outfmt != 'a')
740 <    SET_FILE_BINARY(stdout);
738 >    if (j != EOF) {
739 >        fprintf(stderr, "%s: error on input\n", progname);
740 >        exit(1);
741 >    }
742 >    EPWclose(epw); epw = NULL;
743 >    freedata(mie_dp);
744 >    if (!ntsteps) {
745 >        fprintf(stderr, "%s: no valid time steps on input\n", progname);
746 >        exit(1);
747 >    }
748 >    /* write out matrix */
749 >    if (outfmt != 'a')
750 >        SET_FILE_BINARY(stdout);
751   #ifdef getc_unlocked
752 <  flockfile(stdout);
752 >    flockfile(stdout);
753   #endif
754 <  if (verbose)
755 <    fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
754 >    if (verbose)
755 >        fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
756              outfmt == 'a' ? "" : "binary ", nstored);
757 <  if (doheader) {
758 <    newheader("RADIANCE", stdout);
759 <    printargs(argc, argv, stdout);
760 <    printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude),
761 <           -RadToDeg(s_longitude));
762 <    printf("NROWS=%d\n", nskypatch);
763 <    printf("NCOLS=%d\n", nstored);
764 <    printf("NCOMP=%d\n", NSSAMP);
765 <    if ((outfmt == 'f') | (outfmt == 'd'))
766 <      fputendian(stdout);
767 <    fputformat((char *)getfmtname(outfmt), stdout);
768 <    putchar('\n');
769 <  }
770 <  /* patches are rows (outer sort) */
771 <  for (i = 0; i < nskypatch; i++) {
772 <    mtx_offset = NSSAMP * i;
773 <    switch (outfmt) {
774 <    case 'a':
775 <      for (j = 0; j < nstored; j++) {
776 <        for (int k = 0; k < NSSAMP; k++) {
777 <          printf("%.3g \n", mtx_data[mtx_offset + k]);
757 >    if (doheader) {
758 >        newheader("RADIANCE", stdout);
759 >        printargs(argc, argv, stdout);
760 >        printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude),
761 >            -RadToDeg(s_longitude));
762 >        printf("NROWS=%d\n", nskypatch);
763 >        printf("NCOLS=%d\n", nstored);
764 >        printf("NCOMP=%d\n", NSSAMP);
765 >        if ((outfmt == 'f') | (outfmt == 'd'))
766 >            fputendian(stdout);
767 >        fputformat((char *)getfmtname(outfmt), stdout);
768 >        putchar('\n');
769 >    }
770 >    /* patches are rows (outer sort) */
771 >    for (i = 0; i < nskypatch; i++) {
772 >        mtx_offset = NSSAMP * i;
773 >        switch (outfmt) {
774 >        case 'a':
775 >            for (j = 0; j < nstored; j++) {
776 >                    int k;
777 >                for (k = 0; k < NSSAMP; k++) {
778 >                    printf("%.3g ", mtx_data[mtx_offset + k]);
779 >                }
780 >                printf("\n");
781 >                mtx_offset += NSSAMP * nskypatch;
782 >            }
783 >            if (nstored > 1)
784 >                fputc('\n', stdout);
785 >            break;
786 >        case 'f':
787 >            for (j = 0; j < nstored; j++) {
788 >                putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
789 >                mtx_offset += NSSAMP * nskypatch;
790 >            }
791 >            break;
792 >        case 'd':
793 >            for (j = 0; j < nstored; j++) {
794 >                double ment[NSSAMP];
795 >                for (j = 0; j < NSSAMP; j++)
796 >                    ment[j] = mtx_data[mtx_offset + j];
797 >                putbinary(ment, sizeof(double), NSSAMP, stdout);
798 >                mtx_offset += NSSAMP * nskypatch;
799 >            }
800 >            break;
801          }
802 <        printf("\n");
803 <        mtx_offset += NSSAMP * nskypatch;
676 <      }
677 <      if (nstored > 1)
678 <        fputc('\n', stdout);
679 <      break;
680 <    case 'f':
681 <      for (j = 0; j < nstored; j++) {
682 <        putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
683 <        mtx_offset += NSSAMP * nskypatch;
684 <      }
685 <      break;
686 <    case 'd':
687 <      for (j = 0; j < nstored; j++) {
688 <        double ment[NSSAMP];
689 <        for (j = 0; j < NSSAMP; j++)
690 <          ment[j] = mtx_data[mtx_offset + j];
691 <        putbinary(ment, sizeof(double), NSSAMP, stdout);
692 <        mtx_offset += NSSAMP * nskypatch;
693 <      }
694 <      break;
802 >        if (ferror(stdout))
803 >            goto writerr;
804      }
696 –    if (ferror(stdout))
697 –      goto writerr;
698 –  }
699 – alldone:
700 –  if (fflush(NULL) == EOF)
701 –    goto writerr;
702 –  if (verbose)
703 –    fprintf(stderr, "%s: done.\n", progname);
704 –  exit(0);
805   userr:
806 <  fprintf(stderr,
806 >    fprintf(stderr,
807            "Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-r "
808            "deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
809            progname);
810 <  exit(1);
810 >    exit(1);
811   fmterr:
812 <  fprintf(stderr, "%s: weather tape format error in header\n", progname);
813 <  exit(1);
812 >    fprintf(stderr, "%s: weather tape format error in header\n", progname);
813 >    exit(1);
814   writerr:
815 <  fprintf(stderr, "%s: write error on output\n", progname);
816 <  exit(1);
815 >    fprintf(stderr, "%s: write error on output\n", progname);
816 >    exit(1);
817   }

Diff Legend

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