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.9 by greg, Wed Jul 2 01:56:35 2025 UTC vs.
Revision 1.10 by greg, Fri Jul 11 18:12:25 2025 UTC

# Line 22 | Line 22 | static const char RCSid[] = "$Id$";
22   #include "loadEPW.h"
23  
24  
25 < const   double  SUN_ANG_DEG             = 0.533;                 /* sun full-angle in degrees */
26 < const   double  ARCTIC_LAT              = 67.;
27 < const   double  TROPIC_LAT              = 23.;
28 < const   int             SUMMER_START    = 4;
29 < const   int             SUMMER_END              = 9;
30 < const   double  GNORM                   = 0.777778;
25 > const double SUN_ANG_DEG             = 0.533;                    /* sun full-angle in degrees */
26 > const double ARCTIC_LAT              = 67.;
27 > const double TROPIC_LAT              = 23.;
28 > const int SUMMER_START    = 4;
29 > const int SUMMER_END              = 9;
30 > const double GNORM                   = 0.777778;
31 > const double M_PER_KM        = 1e3;
32  
33   /* Mean normalized relative daylight spectra where CCT = 6415K for overcast */
34 < const   double  D6415[NSSAMP]   = {
34 > const double D6415[NSSAMP]   = {
35          0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
36          1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
37          1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
38 <        0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
38 >        0.82387, 0.87853, 0.82559, 0.75111, 0.78925
39 > };
40  
41   enum {
42 <        NSUNPATCH = 4      /* max. # patches to spread sun into */
42 >        NSUNPATCH = 4          /* max. # patches to spread sun into */
43   };
44  
45 < double  altitude;          /* Solar altitude (radians) */
46 < double  azimuth;                /* Solar azimuth (radians) */
47 < int             julian_date;    /* Julian date */
48 < double  sun_zenith;      /* Sun zenith angle (radians) */
49 < int             nskypatch;        /* number of Reinhart patches */
50 < float   *rh_palt;          /* sky patch altitudes (radians) */
51 < float   *rh_pazi;          /* sky patch azimuths (radians) */
52 < float   *rh_dom;                /* sky patch solid angle (sr) */
53 < FVECT   sundir;
54 < double  sun_ct;          /* cos(theta) of sun altitude angle */
45 > double altitude;           /* Solar altitude (radians) */
46 > double azimuth;                 /* Solar azimuth (radians) */
47 > int julian_date;                /* Julian date */
48 > double sun_zenith;       /* Sun zenith angle (radians) */
49 > int nskypatch;                    /* number of Reinhart patches */
50 > float   *rh_palt;          /* sky patch altitudes (radians) */
51 > float   *rh_pazi;          /* sky patch azimuths (radians) */
52 > float   *rh_dom;                /* sky patch solid angle (sr) */
53 > FVECT sundir;
54 > double sun_ct;           /* cos(theta) of sun altitude angle */
55  
56 < int             input                   = 0;                                    /* Input type */
57 < int             output                  = 0;                                    /* Output type */
58 < int             nsuns                   = NSUNPATCH;                    /* number of sun patches to use */
59 < double  fixed_sun_sa    = -1;                              /* fixed solid angle per sun? */
60 < int             verbose                 = 0;                                    /* progress reports to stderr? */
61 < int             outfmt                  = 'a';                            /* output format */
62 < int             rhsubdiv                = 1;                                    /* Reinhart sky subdivisions */
63 < COLOR   skycolor                = {.96, 1.004, 1.118};  /* sky coloration */
64 < COLOR   suncolor                = {1., 1., 1.};          /* sun color */
65 < double  grefl                   = .2;                              /* ground reflectance */
56 > int input                   = 0;                                                /* Input type */
57 > int output                  = 0;                                                /* Output type */
58 > int nsuns                   = NSUNPATCH;                                /* number of sun patches to use */
59 > double fixed_sun_sa    = -1;                               /* fixed solid angle per sun? */
60 > int verbose                 = 0;                                                /* progress reports to stderr? */
61 > int outfmt                  = 'a';                                        /* output format */
62 > int rhsubdiv                = 1;                                                /* Reinhart sky subdivisions */
63 > COLOR skycolor                = {.96, 1.004, 1.118};    /* sky coloration */
64 > COLOR suncolor                = {1., 1., 1.};            /* sun color */
65 > double grefl                   = .2;                               /* ground reflectance */
66  
67  
68 < static inline double deg_to_rad(double deg)
68 > static inline double
69 > deg_to_rad
70 > (
71 >        double deg
72 > )
73   {
74          return deg * (PI / 180.);
75   }
76  
77 < static inline double rad_to_deg(double rad)
77 > static inline double
78 > rad_to_deg
79 > (
80 >        double rad
81 > )
82   {
83          return rad * (180. / PI);
84   }
85  
86 < static inline void vectorize(double altitude, double azimuth, FVECT v)
86 > static inline void
87 > vectorize
88 > (
89 >        double altitude,
90 >        double azimuth,
91 >        FVECT v
92 > )
93   {
94          v[1] = cos(altitude);
95          v[0] = (v)[1] * sin(azimuth);
# Line 81 | Line 97 | static inline void vectorize(double altitude, double a
97          v[2] = sin(altitude);
98   }
99  
100 < static inline double wmean2(const double a, const double b, const double x)
100 > static inline double
101 > wmean2
102 > (
103 >        const double a,
104 >        const double b,
105 >        const double x
106 > )
107   {
108          return a * (1 - x) + b * x;
109   }
110  
111 < static inline double wmean(
112 <                const double a, const double x,
113 <                const double b, const double y)
111 > static inline double
112 > wmean
113 > (
114 >        const double a,
115 >        const double x,
116 >        const double b,
117 >        const double y
118 > )
119   {
120          return (a * x + b * y) / (a + b);
121   }
122  
123 < static int make_directory(const char *path)
123 > static int
124 > make_directory
125 > (
126 >        const char *path
127 > )
128   {
129   #ifdef _WIN32
130          if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
131                  return 1;
132          }
133          return 0;
134 +
135   #else
136          if (mkdir(path, 0777) == 0 || errno == EEXIST) {
137                  return 1;
138          }
139          return 0;
140 +
141   #endif
142   }
143  
144 < static const char *getfmtname(int fmt)
144 > static const char *
145 > getfmtname
146 > (
147 >        int fmt
148 > )
149   {
150          switch (fmt) {
151          case 'a':
152                  return ("ascii");
153 +
154          case 'f':
155                  return ("float");
156 +
157          case 'd':
158                  return ("double");
159          }
# Line 122 | Line 161 | static const char *getfmtname(int fmt)
161   }
162  
163  
164 < static double get_overcast_zenith_brightness(const double sundir[3])
164 > static double
165 > get_overcast_zenith_brightness
166 > (
167 >        const double sundir[3]
168 > )
169   {
170          double zenithbr;
171          if (sundir[2] < 0) {
# Line 135 | Line 178 | static double get_overcast_zenith_brightness(const dou
178  
179  
180   /* from gensky.c */
181 < static double get_overcast_brightness(const double dz, const double zenithbr)
181 > static double
182 > get_overcast_brightness
183 > (
184 >        const double dz,
185 >        const double zenithbr
186 > )
187   {
188 <  double groundbr = zenithbr * GNORM;
189 <  return wmean(pow(dz + 1.01, 10),
190 <                  zenithbr * (1 + 2 * dz) / 3,
191 <                  pow(dz + 1.01, -10), groundbr);
188 >        double groundbr = zenithbr * GNORM;
189 >        return wmean(pow(dz + 1.01, 10),
190 >                                 zenithbr * (1 + 2 * dz) / 3,
191 >                                 pow(dz + 1.01, -10), groundbr);
192   }
193  
194 < double solar_sunset(int month, int day)
194 > double
195 > solar_sunset
196 > (
197 >        int month,
198 >        int day
199 > )
200   {
201          float W;
202          W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
# Line 151 | Line 204 | double solar_sunset(int month, int day)
204   }
205  
206  
207 < double solar_sunrise(int month, int day)
207 > double
208 > solar_sunrise
209 > (
210 >        int month,
211 >        int day
212 > )
213   {
214          float W;
215          W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
216          return 12 - (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
217   }
218  
219 < int rh_init(void)
219 > int
220 > rh_init
221 > (
222 >        void
223 > )
224   {
225   #define NROW 7
226          static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6};
# Line 166 | Line 228 | int rh_init(void)
228          int p, i, j;
229          /* allocate patch angle arrays */
230          nskypatch = 0;
231 <        for (p = 0; p < NROW; p++)
231 >        for (p = 0; p < NROW; p++) {
232                  nskypatch += tnaz[p];
233 +        }
234          nskypatch *= rhsubdiv * rhsubdiv;
235          nskypatch += 2;
236          rh_palt = (float *)malloc(sizeof(float) * nskypatch);
# Line 177 | Line 240 | int rh_init(void)
240                  fprintf(stderr, "%s: out of memory in rh_init()\n", progname);
241                  exit(1);
242          }
243 <        rh_palt[0] = -PI / 2.; /* ground & zenith patches */
243 >        rh_palt[0] = -PI / 2.;     /* ground & zenith patches */
244          rh_pazi[0] = 0.;
245          rh_dom[0] = 2. * PI;
246          rh_palt[nskypatch - 1] = PI / 2.;
247          rh_pazi[nskypatch - 1] = 0.;
248          rh_dom[nskypatch - 1] = 2. * PI * (1. - cos(alpha * .5));
249 <        p = 1; /* "normal" patches */
249 >        p = 1;     /* "normal" patches */
250          for (i = 0; i < NROW * rhsubdiv; i++) {
251                  const float ralt = alpha * (i + .5);
252                  const int ninrow = tnaz[i / rhsubdiv] * rhsubdiv;
253                  const float dom =
254 <                2. * PI * (sin(alpha * (i + 1)) - sin(alpha * i)) / (double)ninrow;
254 >                        2. * PI * (sin(alpha * (i + 1)) - sin(alpha * i)) / (double)ninrow;
255                  for (j = 0; j < ninrow; j++) {
256                          rh_palt[p] = ralt;
257                          rh_pazi[p] = 2. * PI * j / (double)ninrow;
# Line 196 | Line 259 | int rh_init(void)
259                  }
260          }
261          return nskypatch;
262 +
263   #undef NROW
264   }
265  
266   /* Resize daylight matrix (GW) */
267 < float *resize_dmatrix(float *mtx_data, int nsteps, int npatch)
267 > float *
268 > resize_dmatrix
269 > (
270 >        float *mtx_data,
271 >        int nsteps,
272 >        int npatch
273 > )
274   {
275 <        if (mtx_data == NULL)
275 >        if (mtx_data == NULL) {
276                  mtx_data = (float * ) malloc(sizeof(float) * NSSAMP * nsteps * npatch);
277 <        else
277 >        }else{
278                  mtx_data = (float * ) realloc(mtx_data,
279 <                                sizeof(float) * NSSAMP * nsteps * npatch);
279 >                                                                          sizeof(float) * NSSAMP * nsteps * npatch);
280 >        }
281          if (mtx_data == NULL) {
282                  fprintf(stderr,
283                                  "%s: out of memory in resize_dmatrix(%d,%d)\n",
284                                  progname, nsteps, npatch);
285                  exit(1);
286 <  }
287 <  return mtx_data;
286 >        }
287 >        return mtx_data;
288   }
289  
290 < static Atmosphere init_atmos(const double aod, const double grefl)
290 > static Atmosphere
291 > init_atmos
292 > (
293 >        const double aod,
294 >        const double grefl
295 > )
296   {
297          Atmosphere atmos = {
298                  .ozone_density = {
# Line 256 | Line 332 | static Atmosphere init_atmos(const double aod, const d
332          return atmos;
333   }
334  
335 < static DpPaths get_dppaths(const char *dir, const double aod,
336 <                const char *mname, const char *tag)
335 > static DpPaths
336 > get_dppaths
337 > (
338 >        const char *dir,
339 >        const double aod,
340 >        const char *mname,
341 >        const char *tag
342 > )
343   {
344          DpPaths paths;
345  
346          snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat",
347 <                        dir, DIRSEP, tag, mname, aod);
347 >                         dir, DIRSEP, tag, mname, aod);
348          snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat",
349 <                        dir, DIRSEP, tag, mname, aod);
349 >                         dir, DIRSEP, tag, mname, aod);
350          snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat",
351 <                        dir, DIRSEP, tag, mname, aod);
351 >                         dir, DIRSEP, tag, mname, aod);
352          snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat",
353 <                        dir, DIRSEP, tag, mname, aod);
353 >                         dir, DIRSEP, tag, mname, aod);
354  
355          return paths;
356   }
357  
358  
359 < static void set_rayleigh_density_profile(Atmosphere *atmos,
360 <                char *tag, const int is_summer, const double s_latitude)
359 > static void
360 > set_rayleigh_density_profile
361 > (
362 >        Atmosphere *atmos,
363 >        char *tag,
364 >        const int is_summer,
365 >        const double s_latitude
366 > )
367   {
368          /* Set rayleigh density profile */
369          if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
# Line 311 | Line 399 | static void set_rayleigh_density_profile(Atmosphere *a
399  
400  
401   /* Add in solar direct to nearest sky patches (GW) */
402 < void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m, DATARRAY *irrad,
403 <                double ccover, double dirnorm, float *parr)
402 > void
403 > add_direct
404 > (
405 >        DATARRAY *tau,
406 >        DATARRAY *scat,
407 >        DATARRAY *scat1m,
408 >        DATARRAY *irrad,
409 >        double ccover,
410 >        double dirnorm,
411 >        float *parr
412 > )
413   {
414          FVECT svec;
415          double near_dprod[NSUNPATCH];
# Line 321 | Line 418 | void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRA
418          int i, j, p;
419  
420          /* identify nsuns closest patches */
421 <        for (i = nsuns; i--;)
421 >        for (i = nsuns; i--;) {
422                  near_dprod[i] = -1.;
423 +        }
424          vectorize(altitude, azimuth, svec);
425          for (p = 1; p < nskypatch; p++) {
426                  FVECT pvec;
427                  double dprod;
428                  vectorize(rh_palt[p], rh_pazi[p], pvec);
429                  dprod = DOT(pvec, svec);
430 <                for (i = 0; i < nsuns; i++)
430 >                for (i = 0; i < nsuns; i++) {
431                          if (dprod > near_dprod[i]) {
432                                  for (j = nsuns; --j > i;) {
433                                          near_dprod[j] = near_dprod[j - 1];
# Line 339 | Line 437 | void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRA
437                                  near_patch[i] = p;
438                                  break;
439                          }
440 +                }
441          }
442          /* Get solar radiance */
443          double sun_radiance[NSSAMP] = {0};
# Line 363 | Line 462 | void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRA
462                  intensity = dirnorm / SOLOMG / WHTEFFICACY;
463          }
464          double dir_ratio = 1.;
465 <        if (mean > 0)
465 >        if (mean > 0) {
466                  dir_ratio = intensity / mean;
467 +        }
468          for (i = 0; i < NSSAMP; ++i) {
469                  sun_radiance[i] *= dir_ratio;
470          }
471  
472          /* weight by proximity */
473          wtot = 0;
474 <        for (i = nsuns; i--;)
474 >        for (i = nsuns; i--;) {
475                  wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
476 +        }
477          /* add to nearest patch radiances */
478          for (i = nsuns; i--;) {
479                  float *pdest = parr + NSSAMP * near_patch[i];
# Line 384 | Line 485 | void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRA
485   }
486  
487  
488 < void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m,
489 <                DATARRAY *irrad_clear, double ccover, double dif_ratio,
490 <                double overcast_zenithbr, FVECT view_point, float *parr)
488 > void
489 > calc_sky_patch_radiance
490 > (
491 >        DATARRAY *scat,
492 >        DATARRAY *scat1m,
493 >        DATARRAY *irrad_clear,
494 >        double ccover,
495 >        double dif_ratio,
496 >        double overcast_zenithbr,
497 >        FVECT view_point,
498 >        float *parr
499 > )
500   {
501 <        double  mu_sky; /* Sun-sky point azimuthal angle */
502 <        double  sspa;   /* Sun-sky point angle */
503 <        int      i;
501 >        double mu_sky;      /* Sun-sky point azimuthal angle */
502 >        double sspa;        /* Sun-sky point angle */
503 >        int i;
504          for (i = 1; i < nskypatch; i++) {
505 <                FVECT   rdir_sky;
505 >                FVECT rdir_sky;
506                  vectorize(rh_palt[i], rh_pazi[i], rdir_sky);
507                  mu_sky = fdot(view_point, rdir_sky) / ER;
508                  sspa = fdot(rdir_sky, sundir);
509  
510 <                SCOLOR  sky_radiance = {0};
510 >                SCOLOR sky_radiance = {0};
511                  get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
512 <                int     k;
512 >                int k;
513                  for (k = 0; k < NSSAMP; ++k) {
514                          sky_radiance[k] *= WVLSPAN;
515                  }
# Line 424 | Line 534 | void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY
534  
535  
536   /* Compute sky patch radiance values (modified by GW) */
537 < void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
538 <                DATARRAY *irrad, double ccover, double difhor, FVECT view_point, float *parr)
537 > void
538 > compute_sky
539 > (
540 >        DATARRAY *tau,
541 >        DATARRAY *scat,
542 >        DATARRAY *scat1m,
543 >        DATARRAY *irrad,
544 >        double ccover,
545 >        double difhor,
546 >        FVECT view_point,
547 >        float *parr
548 > )
549   {
550          float sun_zenith;
551          SCOLOR sky_radiance = {0};
# Line 440 | Line 560 | void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARR
560  
561          /* Calculate sun zenith angle (don't let it dip below horizon) */
562          /* Also limit minimum angle to keep circumsolar off zenith */
563 <        if (altitude <= 0.0)
563 >        if (altitude <= 0.0) {
564                  sun_zenith = deg_to_rad(90.0);
565 <        else if (altitude >= deg_to_rad(87.0))
565 >        }else if (altitude >= deg_to_rad(87.0)) {
566                  sun_zenith = deg_to_rad(3.0);
567 <        else
567 >        }else{
568                  sun_zenith = deg_to_rad(90.0) - altitude;
569 +        }
570  
571          double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
572  
# Line 457 | Line 578 | void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARR
578                  double diffuse_irradiance = 0;
579                  int l;
580                  for (l = 0; l < NSSAMP; ++l) {
581 <                        diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20;  /* 20nm interval */
581 >                        diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20;              /* 20nm interval */
582                  }
583                  /* free(indirect_irradiance_clear); */
584                  diffuse_irradiance = wmean2(diffuse_irradiance, overcast_ghi, ccover);
585                  if (diffuse_irradiance > 0) {
586 <                        dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;      /* fudge */
586 >                        dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;                  /* fudge */
587                  }
588          }
589  
# Line 476 | Line 597 | void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARR
597          calc_sky_patch_radiance(scat, scat1m, irrad, ccover, dif_ratio, overcast_zenithbr, view_point, parr);
598   }
599  
600 < int main(int argc, char *argv[])
600 > int
601 > main
602 > (
603 >        int argc,
604 >        char *argv[]
605 > )
606   {
607 <        EPWheader       *epw = NULL;                    /* EPW/WEA input file */
608 <        EPWrecord       erec;                                   /* current EPW/WEA input record */
609 <        int                     doheader = 1;                   /* output header? */
610 <        double          rotation = 0.0;                 /* site rotation (degrees) */
611 <        double          elevation = 0;                  /* site elevation (meters) */
612 <        int                     leap_day = 0;                   /* add leap day? */
613 <        int                     sun_hours_only = 0;             /* only output sun hours? */
614 <        int                     dir_is_horiz;                   /* direct is meas. on horizontal? */
615 <        int                     ntsteps = 0;                    /* number of time steps */
616 <        int                     tstorage = 0;                   /* number of allocated time steps */
617 <        int                     nstored = 0;                    /* number of time steps in matrix */
618 <        int                     last_monthly = 0;               /* month of last report */
619 <        double          dni;                                    /* direct normal illuminance */
620 <        double          dhi;                                    /* diffuse horizontal illuminance */
607 >        EPWheader       *epw = NULL;                        /* EPW/WEA input file */
608 >        EPWrecord erec;                                             /* current EPW/WEA input record */
609 >        int doheader = 1;                                           /* output header? */
610 >        double rotation = 0.0;                              /* site rotation (degrees) */
611 >        double elevation = 0;                               /* site elevation (meters) */
612 >        int leap_day = 0;                                           /* add leap day? */
613 >        int sun_hours_only = 0;                                     /* only output sun hours? */
614 >        int dir_is_horiz;                                           /* direct is meas. on horizontal? */
615 >        int ntsteps = 0;                                            /* number of time steps */
616 >        int tstorage = 0;                                           /* number of allocated time steps */
617 >        int nstored = 0;                                            /* number of time steps in matrix */
618 >        int last_monthly = 0;                                       /* month of last report */
619 >        double dni;                                                 /* direct normal illuminance */
620 >        double dhi;                                                 /* diffuse horizontal illuminance */
621  
622 <        float           *mtx_data = NULL;
623 <        int                     mtx_offset = 0;
624 <        double          timeinterval = 0;
625 <        char            lstag[3];
626 <        char            *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
627 <        char            *ddir = ".";
628 <        char            mie_name[20] = "mie_ca";
629 <        int                     num_threads = 1;
630 <        int                     sorder = 4;
631 <        int                     solar_only = 0;
632 <        int                     sky_only = 0;
633 <        int                     i, j;
634 <        FVECT           view_point      = {0, 0, ER};
622 >        float           *mtx_data = NULL;
623 >        int mtx_offset = 0;
624 >        double timeinterval = 0;
625 >        char lstag[3];
626 >        char            *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
627 >        char            *ddir = ".";
628 >        char mie_name[20] = "mie_ca";
629 >        int num_threads = 1;
630 >        int sorder = 4;
631 >        int solar_only = 0;
632 >        int sky_only = 0;
633 >        int i, j, k;
634 >        FVECT view_point      = {0, 0, ER};
635  
636          fixargv0(argv[0]);
637  
638          for (i = 1; i < argc && argv[i][0] == '-'; i++) {
639                  switch (argv[i][1]) {
640 <                case 'd': /* solar (direct) only */
640 >                case 'd':         /* solar (direct) only */
641                          solar_only = 1;
642                          break;
643 <                case 's': /* sky only (no direct) */
643 >                case 's':         /* sky only (no direct) */
644                          sky_only = 1;
645                          break;
646                  case 'g':
# Line 526 | Line 652 | int main(int argc, char *argv[])
652                  case 'n':
653                          num_threads = atoi(argv[++i]);
654                          break;
655 <                case 'r': /* rotate distribution */
656 <                        if (argv[i][2] && argv[i][2] != 'z')
655 >                case 'r':         /* rotate distribution */
656 >                        if (argv[i][2] && argv[i][2] != 'z') {
657                                  goto userr;
658 +                        }
659                          rotation = atof(argv[++i]);
660                          break;
661 <                case 'u': /* solar hours only */
661 >                case 'u':         /* solar hours only */
662                          sun_hours_only = 1;
663                          break;
664                  case 'p':
665                          ddir = argv[++i];
666                          break;
667 <                case 'v': /* verbose progress reports */
667 >                case 'v':         /* verbose progress reports */
668                          verbose++;
669                          break;
670 <                case 'h': /* turn off header */
670 >                case 'h':         /* turn off header */
671                          doheader = 0;
672                          break;
673 <                case '5': /* 5-phase calculation */
673 >                case '5':         /* 5-phase calculation */
674                          nsuns = 1;
675                          fixed_sun_sa = PI / 360. * atof(argv[++i]);
676                          if (fixed_sun_sa <= 0) {
677                                  fprintf(
678 <                                        stderr,
678 >                                        stderr,
679                                          "%s: missing solar disk size argument for '-5' option\n",
680                                          progname);
681                                  exit(1);
# Line 558 | Line 685 | int main(int argc, char *argv[])
685                  case 'i':
686                          timeinterval = atof(argv[++i]);
687                          break;
688 <                case 'o': /* output format */
688 >                case 'o':         /* output format */
689                          switch (argv[i][2]) {
690                          case 'f':
691                          case 'd':
# Line 573 | Line 700 | int main(int argc, char *argv[])
700                          goto userr;
701                  }
702          }
703 <        if (i < argc - 1)
703 >        if (i < argc - 1) {
704                  goto userr;
705 +        }
706  
707          epw = EPWopen(argv[i]);
708 <        if (epw == NULL)
708 >        if (epw == NULL) {
709                  exit(1);
710 +        }
711          if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
712                  fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
713                  exit(1);
714          }
715          if (verbose) {
716 <                if (i == argc - 1)
716 >                if (i == argc - 1) {
717                          fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
718 <                else
718 >                }else{
719                          fprintf(stderr, "%s: reading weather tape from <stdin>\n", progname);
720 +                }
721          }
722          s_latitude = epw->loc.latitude;
723          s_longitude = -epw->loc.longitude;
724          s_meridian = -15.*epw->loc.timezone;
725          elevation = epw->loc.elevation;
726 <        switch (epw->isWEA) {           /* translate units */
726 >        switch (epw->isWEA) {               /* translate units */
727          case WEAnot:
728          case WEAradnorm:
729 <                input = 1;              /* radiometric quantities */
730 <                dir_is_horiz = 0;       /* direct is perpendicular meas. */
729 >                input = 1;                      /* radiometric quantities */
730 >                dir_is_horiz = 0;               /* direct is perpendicular meas. */
731                  break;
732          case WEAradhoriz:
733 <                input = 1;              /* radiometric quantities */
734 <                dir_is_horiz = 1;       /* solar measured horizontally */
733 >                input = 1;                      /* radiometric quantities */
734 >                dir_is_horiz = 1;               /* solar measured horizontally */
735                  break;
736          case WEAphotnorm:
737 <                input = 2;              /* photometric quantities */
738 <                dir_is_horiz = 0;       /* direct is perpendicular meas. */
737 >                input = 2;                      /* photometric quantities */
738 >                dir_is_horiz = 0;               /* direct is perpendicular meas. */
739                  break;
740          default:
741                  goto fmterr;
# Line 616 | Line 746 | int main(int argc, char *argv[])
746          if (verbose) {
747                  fprintf(stderr, "%s: location '%s %s'\n", progname, epw->loc.city, epw->loc.country);
748                  fprintf(
749 <                        stderr,
749 >                        stderr,
750                          "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
751                          progname, s_latitude, s_longitude);
752 <                if (rotation != 0)
752 >                if (rotation != 0) {
753                          fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
754 +                }
755          }
756  
757          s_latitude = deg_to_rad(s_latitude);
# Line 644 | Line 775 | int main(int argc, char *argv[])
775          }
776  
777          while ((j = EPWread(epw, &erec)) > 0) {
778 <                const int          mo = erec.date.month+1;
779 <                const int          da = erec.date.day;
780 <                const double    hr = erec.date.hour;
781 <                double aod = erec.optdepth * 1e3;
778 >                const int mo = erec.date.month+1;
779 >                const int da = erec.date.day;
780 >                const double hr = erec.date.hour;
781 >                double aod = erec.optdepth * M_PER_KM;
782                  if (aod >= 999.0) {
652                        fprintf(stderr, "aod not set, using default value %.3f\n", AOD0_CA);
783                          aod = AOD0_CA;
784 +                        fprintf(stderr, "aod is not set, using default value %.3f\n", AOD0_CA);
785                  }
786                  double cc = erec.skycover;
787                  if (cc >= 99.0) {
657                        fprintf(stderr, "skycover not set, using default value 0.0\n");
788                          cc = 0.0;
789 +                        fprintf(stderr, "skycover is not set, using default value %.3f\n", 0.0);
790                  }
791 <                double            sda, sta, st;
792 <                int                      sun_in_sky;
791 >                double sda, sta, st;
792 >                int sun_in_sky;
793  
794                  /* compute solar position */
795                  if ((mo == 2) & (da == 29)) {
796                          julian_date = 60;
797                          leap_day = 1;
798 <                } else
798 >                } else{
799                          julian_date = jdate(mo, da) + leap_day;
800 +                }
801                  sda = sdec(julian_date);
802                  sta = stadj(julian_date);
803                  st = hr + sta;
804                  if (timeinterval > 0) {
805 <                        if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120)
805 >                        if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120) {
806                                  st = (st + timeinterval/120 + solar_sunrise(mo, da))/2;
807 <                        else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120)
807 >                        }else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120) {
808                                  st = (st - timeinterval/120 + solar_sunset(mo, da))/2;
809 +                        }
810                  }
811                  altitude = salt(sda, st);
812                  sun_in_sky = (altitude > -deg_to_rad(SUN_ANG_DEG / 2.));
# Line 682 | Line 815 | int main(int argc, char *argv[])
815  
816                  vectorize(altitude, azimuth, sundir);
817                  if (sun_hours_only && !sun_in_sky) {
818 <                        continue; /* skipping nighttime points */
818 >                        continue;             /* skipping nighttime points */
819                  }
820                  sun_ct = fdot(view_point, sundir) / ER;
821  
# Line 697 | Line 830 | int main(int argc, char *argv[])
830                          tstorage += (tstorage >> 1) + nstored + 7;
831                          mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
832                  }
833 <                ntsteps++; /* keep count of time steps */
833 >                ntsteps++;         /* keep count of time steps */
834  
835                  /* compute sky patch values */
836                  Atmosphere clear_atmos = init_atmos(aod, grefl);
# Line 711 | Line 844 | int main(int argc, char *argv[])
844  
845                  char gsdir[PATH_MAX];
846                  size_t siz = strlen(ddir);
847 <                if (ISDIRSEP(ddir[siz - 1]))
847 >                if (ISDIRSEP(ddir[siz - 1])) {
848                          ddir[siz - 1] = '\0';
849 <                        snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
849 >                }
850 >                snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
851                  if (!make_directory(gsdir)) {
852                          fprintf(stderr, "Failed creating atmos_data directory");
853                          exit(1);
# Line 736 | Line 870 | int main(int argc, char *argv[])
870                  DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
871                  DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
872  
873 <                if (!solar_only)
873 >                if (!solar_only) {
874                          compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
875 <                                  cc, dhi, view_point, mtx_data + mtx_offset);
876 <                if (!sky_only)
875 >                                                cc, dhi, view_point, mtx_data + mtx_offset);
876 >                }
877 >                if (!sky_only) {
878                          add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
879 <                                 cc, dni, mtx_data + mtx_offset);
879 >                                           cc, dni, mtx_data + mtx_offset);
880 >                }
881                  /* monthly reporting */
882 <                if (verbose && mo != last_monthly)
882 >                if (verbose && mo != last_monthly) {
883                          fprintf(stderr, "%s: stepping through month %d...\n", progname,
884 <                                last_monthly = mo);
884 >                                        last_monthly = mo);
885 >                }
886          }
887          if (j != EOF) {
888                  fprintf(stderr, "%s: error on input\n", progname);
# Line 758 | Line 895 | int main(int argc, char *argv[])
895                  exit(1);
896          }
897          /* write out matrix */
898 <        if (outfmt != 'a')
898 >        if (outfmt != 'a') {
899                  SET_FILE_BINARY(stdout);
900 +        }
901   #ifdef getc_unlocked
902          flockfile(stdout);
903   #endif
904 <        if (verbose)
904 >        if (verbose) {
905                  fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
906 <                        outfmt == 'a' ? "" : "binary ", nstored);
906 >                                outfmt == 'a' ? "" : "binary ", nstored);
907 >        }
908          if (doheader) {
909                  newheader("RADIANCE", stdout);
910                  printargs(argc, argv, stdout);
911                  printf("LATLONG= %.8f %.8f\n", rad_to_deg(s_latitude),
912 <                        -rad_to_deg(s_longitude));
912 >                           -rad_to_deg(s_longitude));
913                  printf("NROWS=%d\n", nskypatch);
914                  printf("NCOLS=%d\n", nstored);
915                  printf("NCOMP=%d\n", NSSAMP);
916 <                float wvsplit[4] = {380, 480, 588, 780};
917 <                fputwlsplit(wvsplit, stdout);
779 <                if ((outfmt == 'f') | (outfmt == 'd'))
916 >                fputwlsplit(WLPART, stdout);
917 >                if ((outfmt == 'f') | (outfmt == 'd')) {
918                          fputendian(stdout);
919 +                }
920                  fputformat((char *)getfmtname(outfmt), stdout);
921                  putchar('\n');
922          }
# Line 787 | Line 926 | int main(int argc, char *argv[])
926                  switch (outfmt) {
927                  case 'a':
928                          for (j = 0; j < nstored; j++) {
929 <                                int k;
791 <                                for (k = 0; k < NSSAMP; k++) {
929 >                                for (k = NSSAMP - 1; k >= 0; k--) {
930                                          printf("%.3g ", mtx_data[mtx_offset + k]);
931                                  }
932                                  printf("\n");
933                                  mtx_offset += NSSAMP * nskypatch;
934                          }
935 <                        if (nstored > 1)
935 >                        if (nstored > 1) {
936                                  fputc('\n', stdout);
937 +                        }
938                          break;
939                  case 'f':
940                          for (j = 0; j < nstored; j++) {
941 <                                putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
941 >                                float ment[NSSAMP];
942 >                                for (k = NSSAMP - 1; k >= 0; k--) {
943 >                                        ment[k] = mtx_data[mtx_offset + k];
944 >                                }
945 >                                putbinary(ment, sizeof(float), NSSAMP, stdout);
946                                  mtx_offset += NSSAMP * nskypatch;
947                          }
948                          break;
949                  case 'd':
950                          for (j = 0; j < nstored; j++) {
951                                  double ment[NSSAMP];
952 <                                for (j = 0; j < NSSAMP; j++)
953 <                                        ment[j] = mtx_data[mtx_offset + j];
952 >                                for (k = NSSAMP - 1; k >= 0; k--) {
953 >                                        ment[j] = mtx_data[mtx_offset + k];
954 >                                }
955                                  putbinary(ment, sizeof(double), NSSAMP, stdout);
956                                  mtx_offset += NSSAMP * nskypatch;
957                          }
958                          break;
959                  }
960 <                if (ferror(stdout))
960 >                if (ferror(stdout)) {
961                          goto writerr;
962 +                }
963          }
964          return 0;
965 +
966   userr:
967          fprintf(stderr,
968 <                  "Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-r "
969 <                  "deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
970 <                  progname);
968 >                        "Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-r "
969 >                        "deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
970 >                        progname);
971          exit(1);
972   fmterr:
973          fprintf(stderr, "%s: weather tape format error in header\n", progname);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines