ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/radiance/ray/src/gen/gensdaymtx.c

Comparing ray/src/gen/gensdaymtx.c (file contents):
Revision 1.2 by greg, Fri Aug 2 19:01:13 2024 UTC vs.
Revision 1.6 by greg, Tue Apr 15 20:15:50 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	<	char *progname;
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
21	–	double altitude;   /* Solar altitude (radians) */
22	–	double azimuth;    /* Solar azimuth (radians) */
23	–	int julian_date;   /* Julian date */
24	–	double sun_zenith; /* Sun zenith angle (radians) */
25	–	int input = 0;     /* Input type */
26	–	int output = 0;    /* Output type */
27	–	FVECT sundir;
24
25	<	const double ARCTIC_LAT = 67.;
26	<	const double TROPIC_LAT = 23.;
27	<	const int SUMMER_START = 4;
28	<	const int SUMMER_END = 9;
29	<	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
35	–	const double D65EFF = 203.; /* standard illuminant D65 */
36	–
32		/* Mean normalized relative daylight spectra where CCT = 6415K for overcast */
33	<	const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
34	<	1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
35	<	1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
36	<	0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
37	<	/* Degrees into radians */
43	<	#define DegToRad(deg) ((deg) * (PI / 180.))
33	>	const   double  D6415[NSSAMP]   = {
34	>	0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
35	>	1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
36	>	1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
37	>	0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
38
39	<	/* Radiuans into degrees */
40	<	#define RadToDeg(rad) ((rad) * (180. / PI))
39	>	enum {
40	>	NSUNPATCH = 4      /* max. # patches to spread sun into */
41	>	};
42
43	<	#ifndef NSUNPATCH
44	<	#define NSUNPATCH 4 /* max. # patches to spread sun into */
45	<	#endif
43	>	char    *progname;
44	>	double  altitude;          /* Solar altitude (radians) */
45	>	double  azimuth;                /* Solar azimuth (radians) */
46	>	int             julian_date;    /* Julian date */
47	>	double  sun_zenith;      /* Sun zenith angle (radians) */
48	>	int             nskypatch;        /* number of Reinhart patches */
49	>	float   *rh_palt;          /* sky patch altitudes (radians) */
50	>	float   *rh_pazi;          /* sky patch azimuths (radians) */
51	>	float   *rh_dom;                /* sky patch solid angle (sr) */
52	>	FVECT   sundir;
53	>	double  sun_ct;          /* cos(theta) of sun altitude angle */
54
55	<	#define SUN_ANG_DEG 0.533 /* sun full-angle in degrees */
55	>	int             input                   = 0;                                    /* Input type */
56	>	int             output                  = 0;                                    /* Output type */
57	>	int             nsuns                   = NSUNPATCH;                    /* number of sun patches to use */
58	>	double  fixed_sun_sa    = -1;                              /* fixed solid angle per sun? */
59	>	int             verbose                 = 0;                                    /* progress reports to stderr? */
60	>	int             outfmt                  = 'a';                            /* output format */
61	>	int             rhsubdiv                = 1;                                    /* Reinhart sky subdivisions */
62	>	COLOR   skycolor                = {.96, 1.004, 1.118};  /* sky coloration */
63	>	COLOR   suncolor                = {1., 1., 1.};          /* sun color */
64	>	double  grefl                   = .2;                              /* ground reflectance */
65
54	–	int nsuns = NSUNPATCH;    /* number of sun patches to use */
55	–	double fixed_sun_sa = -1; /* fixed solid angle per sun? */
66
67	<	int verbose = 0; /* progress reports to stderr? */
67	>	static inline double deg_to_rad(double deg)
68	>	{
69	>	return deg * (PI / 180.);
70	>	}
71
72	<	int outfmt = 'a'; /* output format */
72	>	static inline double rad_to_deg(double rad)
73	>	{
74	>	return rad * (180. / PI);
75	>	}
76
77	<	int rhsubdiv = 1; /* Reinhart sky subdivisions */
77	>	static inline void vectorize(double altitude, double azimuth, FVECT v)
78	>	{
79	>	v[1] = cos(altitude);
80	>	v[0] = (v)[1] * sin(azimuth);
81	>	v[1] *= cos(azimuth);
82	>	v[2] = sin(altitude);
83	>	}
84
85	<	COLOR skycolor = {.96, 1.004, 1.118}; /* sky coloration */
86	<	COLOR suncolor = {1., 1., 1.};        /* sun color */
87	<	double grefl = .2;                    /* ground reflectance */
85	>	static inline double wmean2(const double a, const double b, const double x)
86	>	{
87	>	return a * (1 - x) + b * x;
88	>	}
89
90	<	int nskypatch;  /* number of Reinhart patches */
91	<	float *rh_palt; /* sky patch altitudes (radians) */
92	<	float *rh_pazi; /* sky patch azimuths (radians) */
93	<	float *rh_dom;  /* sky patch solid angle (sr) */
94	<
72	<	double sun_ct;
73	<
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	<
84	<	inline void vectorize(double altitude, double azimuth, FVECT v) {
85	<	v[1] = cos(altitude);
86	<	v[0] = (v)[1] * sin(azimuth);
87	<	v[1] *= cos(azimuth);
88	<	v[2] = sin(altitude);
90	>	static inline double wmean(
91	>	const double a, const double x,
92	>	const double b, const double y)
93	>	{
94	>	return (a * x + b * y) / (a + b);
95		}
96
97	<	static int make_directory(const char *path) {
97	>	static int make_directory(const char *path)
98	>	{
99		#ifdef _WIN32
100	<	if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
101	<	return 1;
102	<	}
103	<	return 0;
100	>	if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
101	>	return 1;
102	>	}
103	>	return 0;
104		#else
105	<	if (mkdir(path, 0777) == 0 || errno == EEXIST) {
106	<	return 1;
107	<	}
108	<	return 0;
105	>	if (mkdir(path, 0777) == 0 || errno == EEXIST) {
106	>	return 1;
107	>	}
108	>	return 0;
109		#endif
110		}
111
112	<	static const char *getfmtname(int fmt) {
113	<	switch (fmt) {
114	<	case 'a':
115	<	return ("ascii");
116	<	case 'f':
117	<	return ("float");
118	<	case 'd':
119	<	return ("double");
120	<	}
121	<	return ("unknown");
112	>	static const char *getfmtname(int fmt)
113	>	{
114	>	switch (fmt) {
115	>	case 'a':
116	>	return ("ascii");
117	>	case 'f':
118	>	return ("float");
119	>	case 'd':
120	>	return ("double");
121	>	}
122	>	return ("unknown");
123		}
124
117	–	static inline double wmean2(const double a, const double b, const double x) {
118	–	return a * (1 - x) + b * x;
119	–	}
125
126	<	static inline double wmean(const double a, const double x, const double b,
127	<	const double y) {
128	<	return (a * x + b * y) / (a + b);
126	>	static double get_overcast_zenith_brightness(const double sundir[3])
127	>	{
128	>	double zenithbr;
129	>	if (sundir[2] < 0) {
130	>	zenithbr = 0;
131	>	} else {
132	>	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFFICACY;
133	>	}
134	>	return zenithbr;
135		}
136
126	–	static double get_zenith_brightness(const double sundir[3]) {
127	–	double zenithbr;
128	–	if (sundir[2] < 0) {
129	–	zenithbr = 0;
130	–	} else {
131	–	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
132	–	}
133	–	return zenithbr;
134	–	}
137
138		/* from gensky.c */
139	<	static double get_overcast_brightness(const double dz, const double zenithbr) {
139	>	static double get_overcast_brightness(const double dz, const double zenithbr)
140	>	{
141		double groundbr = zenithbr * GNORM;
142	<	return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
143	<	pow(dz + 1.01, -10), groundbr);
142	>	return wmean(pow(dz + 1.01, 10),
143	>	zenithbr * (1 + 2 * dz) / 3,
144	>	pow(dz + 1.01, -10), groundbr);
145		}
146
147	<	int rh_init(void) {
147	>	double solar_sunset(int month, int day)
148	>	{
149	>	float W;
150	>	W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
151	>	return 12 + (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
152	>	}
153	>
154	>
155	>	double solar_sunrise(int month, int day)
156	>	{
157	>	float W;
158	>	W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
159	>	return 12 - (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
160	>	}
161	>
162	>	int rh_init(void)
163	>	{
164		#define NROW 7
165	<	static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6};
166	<	const double alpha = (PI / 2.) / (NROW * rhsubdiv + .5);
167	<	int p, i, j;
168	<	/* allocate patch angle arrays */
169	<	nskypatch = 0;
170	<	for (p = 0; p < NROW; p++)
171	<	nskypatch += tnaz[p];
172	<	nskypatch *= rhsubdiv * rhsubdiv;
173	<	nskypatch += 2;
174	<	rh_palt = (float *)malloc(sizeof(float) * nskypatch);
175	<	rh_pazi = (float *)malloc(sizeof(float) * nskypatch);
176	<	rh_dom = (float *)malloc(sizeof(float) * nskypatch);
177	<	if ((rh_palt == NULL) | (rh_pazi == NULL) | (rh_dom == NULL)) {
178	<	fprintf(stderr, "%s: out of memory in rh_init()\n", progname);
179	<	exit(1);
180	<	}
181	<	rh_palt[0] = -PI / 2.; /* ground & zenith patches */
182	<	rh_pazi[0] = 0.;
183	<	rh_dom[0] = 2. * PI;
184	<	rh_palt[nskypatch - 1] = PI / 2.;
185	<	rh_pazi[nskypatch - 1] = 0.;
186	<	rh_dom[nskypatch - 1] = 2. * PI * (1. - cos(alpha * .5));
187	<	p = 1; /* "normal" patches */
188	<	for (i = 0; i < NROW * rhsubdiv; i++) {
189	<	const float ralt = alpha * (i + .5);
190	<	const int ninrow = tnaz[i / rhsubdiv] * rhsubdiv;
191	<	const float dom =
192	<	2. * PI * (sin(alpha * (i + 1)) - sin(alpha * i)) / (double)ninrow;
193	<	for (j = 0; j < ninrow; j++) {
194	<	rh_palt[p] = ralt;
195	<	rh_pazi[p] = 2. * PI * j / (double)ninrow;
196	<	rh_dom[p++] = dom;
197	<	}
198	<	}
199	<	return nskypatch;
165	>	static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6};
166	>	const double alpha = (PI / 2.) / (NROW * rhsubdiv + .5);
167	>	int p, i, j;
168	>	/* allocate patch angle arrays */
169	>	nskypatch = 0;
170	>	for (p = 0; p < NROW; p++)
171	>	nskypatch += tnaz[p];
172	>	nskypatch *= rhsubdiv * rhsubdiv;
173	>	nskypatch += 2;
174	>	rh_palt = (float *)malloc(sizeof(float) * nskypatch);
175	>	rh_pazi = (float *)malloc(sizeof(float) * nskypatch);
176	>	rh_dom = (float *)malloc(sizeof(float) * nskypatch);
177	>	if ((rh_palt == NULL) | (rh_pazi == NULL) | (rh_dom == NULL)) {
178	>	fprintf(stderr, "%s: out of memory in rh_init()\n", progname);
179	>	exit(1);
180	>	}
181	>	rh_palt[0] = -PI / 2.; /* ground & zenith patches */
182	>	rh_pazi[0] = 0.;
183	>	rh_dom[0] = 2. * PI;
184	>	rh_palt[nskypatch - 1] = PI / 2.;
185	>	rh_pazi[nskypatch - 1] = 0.;
186	>	rh_dom[nskypatch - 1] = 2. * PI * (1. - cos(alpha * .5));
187	>	p = 1; /* "normal" patches */
188	>	for (i = 0; i < NROW * rhsubdiv; i++) {
189	>	const float ralt = alpha * (i + .5);
190	>	const int ninrow = tnaz[i / rhsubdiv] * rhsubdiv;
191	>	const float dom =
192	>	2. * PI * (sin(alpha * (i + 1)) - sin(alpha * i)) / (double)ninrow;
193	>	for (j = 0; j < ninrow; j++) {
194	>	rh_palt[p] = ralt;
195	>	rh_pazi[p] = 2. * PI * j / (double)ninrow;
196	>	rh_dom[p++] = dom;
197	>	}
198	>	}
199	>	return nskypatch;
200		#undef NROW
201		}
202
203		/* Resize daylight matrix (GW) */
204	<	float *resize_dmatrix(float *mtx_data, int nsteps, int npatch) {
205	<	if (mtx_data == NULL)
206	<	mtx_data = (float *)malloc(sizeof(float) * NSSAMP * nsteps * npatch);
207	<	else
208	<	mtx_data =
209	<	(float *)realloc(mtx_data, sizeof(float) * NSSAMP * nsteps * npatch);
210	<	if (mtx_data == NULL) {
211	<	fprintf(stderr, "%s: out of memory in resize_dmatrix(%d,%d)\n", progname,
212	<	nsteps, npatch);
213	<	exit(1);
204	>	float *resize_dmatrix(float *mtx_data, int nsteps, int npatch)
205	>	{
206	>	if (mtx_data == NULL)
207	>	mtx_data = (float * ) malloc(sizeof(float) * NSSAMP * nsteps * npatch);
208	>	else
209	>	mtx_data = (float * ) realloc(mtx_data,
210	>	sizeof(float) * NSSAMP * nsteps * npatch);
211	>	if (mtx_data == NULL) {
212	>	fprintf(stderr,
213	>	"%s: out of memory in resize_dmatrix(%d,%d)\n",
214	>	progname, nsteps, npatch);
215	>	exit(1);
216		}
217	<	return (mtx_data);
217	>	return mtx_data;
218		}
219
220	<	static Atmosphere init_atmos(const double aod, const double grefl) {
221	<	Atmosphere atmos = {.ozone_density = {.layers =
222	<	{
223	<	{.width = 25000.0,
224	<	.exp_term = 0.0,
225	<	.exp_scale = 0.0,
226	<	.linear_term = 1.0 / 15000.0,
227	<	.constant_term = -2.0 / 3.0},
228	<	{.width = AH,
229	<	.exp_term = 0.0,
230	<	.exp_scale = 0.0,
231	<	.linear_term = -1.0 / 15000.0,
232	<	.constant_term = 8.0 / 3.0},
233	<	}},
234	<	.rayleigh_density = {.layers =
235	<	{
236	<	{.width = AH,
237	<	.exp_term = 1.0,
238	<	.exp_scale = -1.0 / HR_MS,
239	<	.linear_term = 0.0,
240	<	.constant_term = 0.0},
241	<	}},
242	<	.beta_r0 = BR0_MS,
243	<	.beta_scale = aod / AOD0_CA,
244	<	.beta_m = NULL,
245	<	.grefl = grefl};
246	<	return atmos;
220	>	static Atmosphere init_atmos(const double aod, const double grefl)
221	>	{
222	>	Atmosphere atmos = {
223	>	.ozone_density = {
224	>	.layers = {
225	>	{
226	>	.width = 25000.0,
227	>	.exp_term = 0.0,
228	>	.exp_scale = 0.0,
229	>	.linear_term = 1.0 / 15000.0,
230	>	.constant_term = -2.0 / 3.0
231	>	},
232	>	{
233	>	.width = AH,
234	>	.exp_term = 0.0,
235	>	.exp_scale = 0.0,
236	>	.linear_term = -1.0 / 15000.0,
237	>	.constant_term = 8.0 / 3.0
238	>	},
239	>	}
240	>	},
241	>	.rayleigh_density = {
242	>	.layers = {
243	>	{
244	>	.width = AH,
245	>	.exp_term = 1.0,
246	>	.exp_scale = -1.0 / HR_MS,
247	>	.linear_term = 0.0,
248	>	.constant_term = 0.0
249	>	},
250	>	}
251	>	},
252	>	.beta_r0 = BR0_MS,
253	>	.beta_scale = aod / AOD0_CA,
254	>	.beta_m = NULL,
255	>	.grefl = grefl
256	>	};
257	>	return atmos;
258		}
259
260	<	static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
261	<	const char *tag) {
262	<	DpPaths paths;
260	>	static DpPaths get_dppaths(const char *dir, const double aod,
261	>	const char *mname, const char *tag)
262	>	{
263	>	DpPaths paths;
264
265	<	snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
266	<	mname, aod);
267	<	snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
268	<	mname, aod);
269	<	snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
270	<	tag, mname, aod);
271	<	snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
272	<	mname, aod);
265	>	snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat",
266	>	dir, DIRSEP, tag, mname, aod);
267	>	snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat",
268	>	dir, DIRSEP, tag, mname, aod);
269	>	snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat",
270	>	dir, DIRSEP, tag, mname, aod);
271	>	snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat",
272	>	dir, DIRSEP, tag, mname, aod);
273
274	<	return paths;
274	>	return paths;
275		}
276	<	static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
277	<	const int is_summer,
278	<	const double s_latitude) {
279	<	/* Set rayleigh density profile */
280	<	if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
281	<	tag[0] = 's';
282	<	if (is_summer) {
283	<	tag[1] = 's';
284	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
285	<	atmos->beta_r0 = BR0_SS;
286	<	} else {
287	<	tag[1] = 'w';
288	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
289	<	atmos->beta_r0 = BR0_SW;
290	<	}
291	<	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
292	<	tag[0] = 'm';
293	<	if (is_summer) {
294	<	tag[1] = 's';
295	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
296	<	atmos->beta_r0 = BR0_MS;
297	<	} else {
298	<	tag[1] = 'w';
299	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
300	<	atmos->beta_r0 = BR0_MW;
301	<	}
302	<	} else {
303	<	tag[0] = 't';
304	<	tag[1] = 'r';
305	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
306	<	atmos->beta_r0 = BR0_T;
307	<	}
308	<	tag[2] = '\0';
276	>
277	>
278	>	static void set_rayleigh_density_profile(Atmosphere *atmos,
279	>	char *tag, const int is_summer, const double s_latitude)
280	>	{
281	>	/* Set rayleigh density profile */
282	>	if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
283	>	tag[0] = 's';
284	>	if (is_summer) {
285	>	tag[1] = 's';
286	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
287	>	atmos->beta_r0 = BR0_SS;
288	>	} else {
289	>	tag[1] = 'w';
290	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
291	>	atmos->beta_r0 = BR0_SW;
292	>	}
293	>	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
294	>	tag[0] = 'm';
295	>	if (is_summer) {
296	>	tag[1] = 's';
297	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
298	>	atmos->beta_r0 = BR0_MS;
299	>	} else {
300	>	tag[1] = 'w';
301	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
302	>	atmos->beta_r0 = BR0_MW;
303	>	}
304	>	} else {
305	>	tag[0] = 't';
306	>	tag[1] = 'r';
307	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
308	>	atmos->beta_r0 = BR0_T;
309	>	}
310	>	tag[2] = '\0';
311		}
312	+
313	+
314		/* Add in solar direct to nearest sky patches (GW) */
315	<	void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
316	<	DATARRAY *irrad, double ccover, float *parr) {
317	<	FVECT svec;
318	<	double near_dprod[NSUNPATCH];
319	<	int near_patch[NSUNPATCH];
320	<	double wta[NSUNPATCH], wtot;
321	<	int i, j, p;
315	>	void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m, DATARRAY *irrad,
316	>	double ccover, double dirnorm, float *parr)
317	>	{
318	>	FVECT svec;
319	>	double near_dprod[NSUNPATCH];
320	>	int near_patch[NSUNPATCH];
321	>	double wta[NSUNPATCH], wtot;
322	>	int i, j, p;
323
324	<	/* identify nsuns closest patches */
325	<	for (i = nsuns; i--;)
326	<	near_dprod[i] = -1.;
327	<	vectorize(altitude, azimuth, svec);
328	<	for (p = 1; p < nskypatch; p++) {
329	<	FVECT pvec;
330	<	double dprod;
331	<	vectorize(rh_palt[p], rh_pazi[p], pvec);
332	<	dprod = DOT(pvec, svec);
333	<	for (i = 0; i < nsuns; i++)
334	<	if (dprod > near_dprod[i]) {
335	<	for (j = nsuns; --j > i;) {
336	<	near_dprod[j] = near_dprod[j - 1];
337	<	near_patch[j] = near_patch[j - 1];
338	<	}
339	<	near_dprod[i] = dprod;
340	<	near_patch[i] = p;
341	<	break;
342	<	}
343	<	}
344	<	/* Get solar radiance */
345	<	double sun_radiance[NSSAMP] = {0};
346	<	get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
347	<	if (ccover > 0) {
348	<	double zenithbr = get_zenith_brightness(sundir);
349	<	double skybr = get_overcast_brightness(sundir[2], zenithbr);
350	<	for (int l = 0; l < NSSAMP; ++l) {
351	<	sun_radiance[l] =
352	<	wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
353	<	}
354	<	}
355	<	/* weight by proximity */
356	<	wtot = 0;
357	<	for (i = nsuns; i--;)
358	<	wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
359	<	/* add to nearest patch radiances */
360	<	for (i = nsuns; i--;) {
361	<	float *pdest = parr + NSSAMP * near_patch[i];
362	<	for (int k = 0; k < NSSAMP; k++) {
363	<	*pdest++ = sun_radiance[k] * wta[i] / wtot;
364	<	}
365	<	}
324	>	/* identify nsuns closest patches */
325	>	for (i = nsuns; i--;)
326	>	near_dprod[i] = -1.;
327	>	vectorize(altitude, azimuth, svec);
328	>	for (p = 1; p < nskypatch; p++) {
329	>	FVECT pvec;
330	>	double dprod;
331	>	vectorize(rh_palt[p], rh_pazi[p], pvec);
332	>	dprod = DOT(pvec, svec);
333	>	for (i = 0; i < nsuns; i++)
334	>	if (dprod > near_dprod[i]) {
335	>	for (j = nsuns; --j > i;) {
336	>	near_dprod[j] = near_dprod[j - 1];
337	>	near_patch[j] = near_patch[j - 1];
338	>	}
339	>	near_dprod[i] = dprod;
340	>	near_patch[i] = p;
341	>	break;
342	>	}
343	>	}
344	>	/* Get solar radiance */
345	>	double sun_radiance[NSSAMP] = {0};
346	>	get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
347	>	if (ccover > 0) {
348	>	double zenithbr = get_overcast_zenith_brightness(sundir);
349	>	double skybr = get_overcast_brightness(sundir[2], zenithbr);
350	>	int l;
351	>	for (l = 0; l < NSSAMP; ++l) {
352	>	sun_radiance[l] = wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
353	>	}
354	>	}
355	>	/* Normalize */
356	>	double sum = 0.0;
357	>	for (i = 0; i < NSSAMP; ++i) {
358	>	sum += sun_radiance[i];
359	>	}
360	>	double mean = sum / NSSAMP;
361	>
362	>	double intensity = mean * WVLSPAN;
363	>	if (dirnorm > 0) {
364	>	intensity = dirnorm / SOLOMG / WHTEFFICACY;
365	>	}
366	>	double dir_ratio = 1.;
367	>	if (mean > 0)
368	>	dir_ratio = intensity / mean;
369	>	for (i = 0; i < NSSAMP; ++i) {
370	>	sun_radiance[i] *= dir_ratio;
371	>	}
372	>
373	>	/* weight by proximity */
374	>	wtot = 0;
375	>	for (i = nsuns; i--;)
376	>	wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
377	>	/* add to nearest patch radiances */
378	>	for (i = nsuns; i--;) {
379	>	float *pdest = parr + NSSAMP * near_patch[i];
380	>	int k;
381	>	for (k = 0; k < NSSAMP; k++) {
382	>	*pdest++ = sun_radiance[k] * wta[i] / wtot;
383	>	}
384	>	}
385		}
386
329	–	void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m, double ccover,
330	–	float *parr) {
331	–	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};
387
388	<	get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
389	<	for (int k = 0; k < NSSAMP; ++k) {
390	<	sky_radiance[k] *= WVLSPAN;
391	<	}
388	>	void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m,
389	>	DATARRAY *irrad_clear, double ccover, double dif_ratio,
390	>	double overcast_zenithbr, FVECT view_point, float *parr)
391	>	{
392	>	double  mu_sky; /* Sun-sky point azimuthal angle */
393	>	double  sspa;   /* Sun-sky point angle */
394	>	int      i;
395	>	for (i = 1; i < nskypatch; i++) {
396	>	FVECT   rdir_sky;
397	>	vectorize(rh_palt[i], rh_pazi[i], rdir_sky);
398	>	mu_sky = fdot(view_point, rdir_sky) / ER;
399	>	sspa = fdot(rdir_sky, sundir);
400
401	<	if (ccover > 0) {
402	<	double zenithbr = get_zenith_brightness(sundir);
403	<	double grndbr = zenithbr * GNORM;
404	<	double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
405	<	for (int k = 0; k < NSSAMP; ++k) {
406	<	sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
354	<	}
355	<	}
401	>	SCOLOR  sky_radiance = {0};
402	>	get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
403	>	int     k;
404	>	for (k = 0; k < NSSAMP; ++k) {
405	>	sky_radiance[k] *= WVLSPAN;
406	>	}
407
408	<	for (int k = 0; k < NSSAMP; ++k) {
409	<	parr[NSSAMP * i + k] = sky_radiance[k];
410	<	}
411	<	}
408	>	if (ccover > 0) {
409	>	double skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr);
410	>	for (k = 0; k < NSSAMP; ++k) {
411	>	sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
412	>	}
413	>	}
414	>
415	>	/* calibration */
416	>	for (k = 0; k < NSSAMP; ++k) {
417	>	sky_radiance[k] *= dif_ratio;
418	>	}
419	>
420	>	for (k = 0; k < NSSAMP; ++k) {
421	>	parr[NSSAMP * i + k] = sky_radiance[k];
422	>	}
423	>	}
424		}
425
363	–	/* Return maximum of two doubles */
364	–	static inline double dmax(double a, double b) { return (a > b) ? a : b; }
426
427		/* Compute sky patch radiance values (modified by GW) */
428		void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
429	<	DATARRAY *irrad, double ccover, float *parr) {
430	<	int index; /* Category index */
431	<	int i;
432	<	float sun_zenith;
433	<	SCOLOR sky_radiance = {0};
434	<	SCOLOR ground_radiance = {0};
435	<	SCOLR sky_sclr = {0};
436	<	SCOLR ground_sclr = {0};
437	<	FVECT view_point = {0, 0, ER};
438	<	const double radius = VLEN(view_point);
439	<	const double sun_ct = fdot(view_point, sundir) / radius;
440	<	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);
429	>	DATARRAY *irrad, double ccover, double difhor, FVECT view_point, float *parr)
430	>	{
431	>	float sun_zenith;
432	>	SCOLOR sky_radiance = {0};
433	>	SCOLOR ground_radiance = {0};
434	>	SCOLR sky_sclr = {0};
435	>	SCOLR ground_sclr = {0};
436	>	const double radius = VLEN(view_point);
437	>	const double sun_ct = fdot(view_point, sundir) / radius;
438	>	const FVECT rdir_grnd = {0, 0, -1};
439	>	const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
440	>	const double nu_grnd = fdot(rdir_grnd, sundir);
441
442	<	/* Calculate sun zenith angle (don't let it dip below horizon) */
443	<	/* Also limit minimum angle to keep circumsolar off zenith */
444	<	if (altitude <= 0.0)
445	<	sun_zenith = DegToRad(90.0);
446	<	else if (altitude >= DegToRad(87.0))
447	<	sun_zenith = DegToRad(3.0);
448	<	else
449	<	sun_zenith = DegToRad(90.0) - altitude;
442	>	/* Calculate sun zenith angle (don't let it dip below horizon) */
443	>	/* Also limit minimum angle to keep circumsolar off zenith */
444	>	if (altitude <= 0.0)
445	>	sun_zenith = deg_to_rad(90.0);
446	>	else if (altitude >= deg_to_rad(87.0))
447	>	sun_zenith = deg_to_rad(3.0);
448	>	else
449	>	sun_zenith = deg_to_rad(90.0) - altitude;
450
451	<	/* Compute ground radiance (include solar contribution if any) */
452	<	get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
453	<	mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
454	<	for (int j = 0; j < NSSAMP; j++) {
455	<	parr[j] *= WVLSPAN;
456	<	}
457	<	calc_sky_patch_radiance(scat, scat1m, ccover, parr);
451	>	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
452	>
453	>	/* diffuse calibration factor */
454	>	double dif_ratio = 1;
455	>	if (difhor > 0) {
456	>	DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad, radius, sun_ct);
457	>	double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
458	>	double diffuse_irradiance = 0;
459	>	int l;
460	>	for (l = 0; l < NSSAMP; ++l) {
461	>	diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20;  /* 20nm interval */
462	>	}
463	>	/* free(indirect_irradiance_clear); */
464	>	diffuse_irradiance = wmean2(diffuse_irradiance, overcast_ghi, ccover);
465	>	if (diffuse_irradiance > 0) {
466	>	dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;      /* fudge */
467	>	}
468	>	}
469	>
470	>	/* Compute ground radiance (include solar contribution if any) */
471	>	get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
472	>	mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
473	>	int j;
474	>	for (j = 0; j < NSSAMP; j++) {
475	>	parr[j] *= WVLSPAN;
476	>	}
477	>	calc_sky_patch_radiance(scat, scat1m, irrad, ccover, dif_ratio, overcast_zenithbr, view_point, parr);
478		}
479
480	<	int main(int argc, char *argv[]) {
480	>	int main(int argc, char *argv[])
481	>	{
482	>	EPWheader       *epw = NULL;                    /* EPW/WEA input file */
483	>	EPWrecord       erec;                                   /* current EPW/WEA input record */
484	>	int                     doheader = 1;                   /* output header? */
485	>	double          rotation = 0.0;                 /* site rotation (degrees) */
486	>	double          elevation = 0;                  /* site elevation (meters) */
487	>	int                     leap_day = 0;                   /* add leap day? */
488	>	int                     sun_hours_only = 0;             /* only output sun hours? */
489	>	int                     dir_is_horiz;                   /* direct is meas. on horizontal? */
490	>	int                     ntsteps = 0;                    /* number of time steps */
491	>	int                     tstorage = 0;                   /* number of allocated time steps */
492	>	int                     nstored = 0;                    /* number of time steps in matrix */
493	>	int                     last_monthly = 0;               /* month of last report */
494	>	double          dni;                                    /* direct normal illuminance */
495	>	double          dhi;                                    /* diffuse horizontal illuminance */
496
497	<	char buf[256];
498	<	int doheader = 1; /* output header? */
499	<	double rotation = 0.0;
500	<	double elevation = 0;
501	<	int leap_day = 0;       /* add leap day? */
502	<	int sun_hours_only = 0; /* only output sun hours? */
503	<	float *mtx_data = NULL;
504	<	int ntsteps = 0;      /* number of time steps */
505	<	int tstorage = 0;     /* number of allocated time steps */
506	<	int nstored = 0;      /* number of time steps in matrix */
507	<	int last_monthly = 0; /* month of last report */
508	<	int mo, da;
509	<	double hr, aod, cc;
416	<	double dni, dhi;
417	<	int mtx_offset = 0;
418	<	int i, j;
419	<	char lstag[3];
420	<	char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
421	<	char *ddir = ".";
422	<	char mie_name[20] = "mie_ca";
423	<	int num_threads = 1;
424	<	int sorder = 4;
425	<	int solar_only = 0;
426	<	int sky_only = 0;
427	<	FVECT view_point = {0, 0, ER};
497	>	float           *mtx_data = NULL;
498	>	int                     mtx_offset = 0;
499	>	double          timeinterval = 0;
500	>	char            lstag[3];
501	>	char            *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
502	>	char            *ddir = ".";
503	>	char            mie_name[20] = "mie_ca";
504	>	int                     num_threads = 1;
505	>	int                     sorder = 4;
506	>	int                     solar_only = 0;
507	>	int                     sky_only = 0;
508	>	int                     i, j;
509	>	FVECT           view_point      = {0, 0, ER};
510
511	<	progname = argv[0];
511	>	progname = argv[0];
512
513	<	for (i = 1; i < argc && argv[i][0] == '-'; i++) {
514	<	switch (argv[i][1]) {
515	<	case 'd': /* solar (direct) only */
516	<	solar_only = 1;
517	<	break;
518	<	case 's': /* sky only (no direct) */
519	<	sky_only = 1;
520	<	break;
521	<	case 'g':
522	<	grefl = atof(argv[++i]);
523	<	break;
524	<	case 'm':
525	<	rhsubdiv = atoi(argv[++i]);
526	<	break;
527	<	case 'n':
528	<	num_threads = atoi(argv[++i]);
529	<	break;
530	<	case 'r': /* rotate distribution */
531	<	if (argv[i][2] && argv[i][2] != 'z')
532	<	goto userr;
533	<	rotation = atof(argv[++i]);
534	<	break;
535	<	case 'u': /* solar hours only */
536	<	sun_hours_only = 1;
537	<	break;
538	<	case 'p':
539	<	ddir = argv[++i];
540	<	break;
541	<	case 'v': /* verbose progress reports */
542	<	verbose++;
543	<	break;
544	<	case 'h': /* turn off header */
545	<	doheader = 0;
546	<	break;
547	<	case '5': /* 5-phase calculation */
548	<	nsuns = 1;
549	<	fixed_sun_sa = PI / 360. * atof(argv[++i]);
550	<	if (fixed_sun_sa <= 0) {
551	<	fprintf(stderr,
552	<	"%s: missing solar disk size argument for '-5' option\n",
553	<	progname);
554	<	exit(1);
555	<	}
556	<	fixed_sun_sa *= fixed_sun_sa * PI;
557	<	break;
558	<	case 'o': /* output format */
559	<	switch (argv[i][2]) {
560	<	case 'f':
561	<	case 'd':
562	<	case 'a':
563	<	outfmt = argv[i][2];
564	<	break;
565	<	default:
566	<	goto userr;
567	<	}
568	<	break;
569	<	default:
570	<	goto userr;
571	<	}
572	<	}
573	<	if (i < argc - 1)
574	<	goto userr;
575	<	if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
576	<	fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
577	<	exit(1);
578	<	}
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;
513	>	for (i = 1; i < argc && argv[i][0] == '-'; i++) {
514	>	switch (argv[i][1]) {
515	>	case 'd': /* solar (direct) only */
516	>	solar_only = 1;
517	>	break;
518	>	case 's': /* sky only (no direct) */
519	>	sky_only = 1;
520	>	break;
521	>	case 'g':
522	>	grefl = atof(argv[++i]);
523	>	break;
524	>	case 'm':
525	>	rhsubdiv = atoi(argv[++i]);
526	>	break;
527	>	case 'n':
528	>	num_threads = atoi(argv[++i]);
529	>	break;
530	>	case 'r': /* rotate distribution */
531	>	if (argv[i][2] && argv[i][2] != 'z')
532	>	goto userr;
533	>	rotation = atof(argv[++i]);
534	>	break;
535	>	case 'u': /* solar hours only */
536	>	sun_hours_only = 1;
537	>	break;
538	>	case 'p':
539	>	ddir = argv[++i];
540	>	break;
541	>	case 'v': /* verbose progress reports */
542	>	verbose++;
543	>	break;
544	>	case 'h': /* turn off header */
545	>	doheader = 0;
546	>	break;
547	>	case '5': /* 5-phase calculation */
548	>	nsuns = 1;
549	>	fixed_sun_sa = PI / 360. * atof(argv[++i]);
550	>	if (fixed_sun_sa <= 0) {
551	>	fprintf(
552	>	stderr,
553	>	"%s: missing solar disk size argument for '-5' option\n",
554	>	progname);
555	>	exit(1);
556	>	}
557	>	fixed_sun_sa *= fixed_sun_sa * PI;
558	>	break;
559	>	case 'i':
560	>	timeinterval = atof(argv[++i]);
561	>	break;
562	>	case 'o': /* output format */
563	>	switch (argv[i][2]) {
564	>	case 'f':
565	>	case 'd':
566	>	case 'a':
567	>	outfmt = argv[i][2];
568	>	break;
569	>	default:
570	>	goto userr;
571	>	}
572	>	break;
573	>	default:
574	>	goto userr;
575	>	}
576	>	}
577	>	if (i < argc - 1)
578	>	goto userr;
579
580	<	rh_init();
581	<	if (verbose) {
582	<	fprintf(stderr, "%s: location '%s'\n", progname, buf);
583	<	fprintf(stderr, "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
584	<	progname, s_latitude, s_longitude);
585	<	if (rotation != 0)
586	<	fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
587	<	}
580	>	epw = EPWopen(argv[i]);
581	>	if (epw == NULL)
582	>	exit(1);
583	>	if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
584	>	fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
585	>	exit(1);
586	>	}
587	>	if (verbose) {
588	>	if (i == argc - 1)
589	>	fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
590	>	else
591	>	fprintf(stderr, "%s: reading weather tape from <stdin>\n", progname);
592	>	}
593	>	s_latitude = epw->loc.latitude;
594	>	s_longitude = -epw->loc.longitude;
595	>	s_meridian = -15.*epw->loc.timezone;
596	>	elevation = epw->loc.elevation;
597	>	switch (epw->isWEA) {           /* translate units */
598	>	case WEAnot:
599	>	case WEAradnorm:
600	>	input = 1;              /* radiometric quantities */
601	>	dir_is_horiz = 0;       /* direct is perpendicular meas. */
602	>	break;
603	>	case WEAradhoriz:
604	>	input = 1;              /* radiometric quantities */
605	>	dir_is_horiz = 1;       /* solar measured horizontally */
606	>	break;
607	>	case WEAphotnorm:
608	>	input = 2;              /* photometric quantities */
609	>	dir_is_horiz = 0;       /* direct is perpendicular meas. */
610	>	break;
611	>	default:
612	>	goto fmterr;
613	>	}
614
615	<	s_latitude = DegToRad(s_latitude);
527	<	s_longitude = DegToRad(s_longitude);
528	<	s_meridian = DegToRad(s_meridian);
529	<	/* initial allocation */
530	<	mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
615	>	rh_init();
616
617	<	/* Load mie density data */
618	<	DATARRAY *mie_dp = getdata(mie_path);
619	<	if (mie_dp == NULL) {
620	<	fprintf(stderr, "Error reading mie data\n");
621	<	return 0;
622	<	}
617	>	if (verbose) {
618	>	fprintf(stderr, "%s: location '%s %s'\n", progname, epw->loc.city, epw->loc.country);
619	>	fprintf(
620	>	stderr,
621	>	"%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
622	>	progname, s_latitude, s_longitude);
623	>	if (rotation != 0)
624	>	fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
625	>	}
626
627	<	while (scanf("%d %d %lf %lf %lf %lf %lf\n", &mo, &da, &hr, &dni, &dhi, &aod,
628	<	&cc) == 7) {
629	<	double sda, sta;
630	<	int sun_in_sky;
631	<	/* compute solar position */
544	<	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.));
627	>	s_latitude = deg_to_rad(s_latitude);
628	>	s_longitude = deg_to_rad(s_longitude);
629	>	s_meridian = deg_to_rad(s_meridian);
630	>	/* initial allocation */
631	>	mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
632
633	<	azimuth = sazi(sda, hr + sta) + PI - DegToRad(rotation);
633	>	/* Load mie density data */
634	>	DATARRAY *mie_dp = getdata(mie_path);
635	>	if (mie_dp == NULL) {
636	>	fprintf(stderr, "Error reading mie data\n");
637	>	return 0;
638	>	}
639
640	<	vectorize(altitude, azimuth, sundir);
641	<	if (sun_hours_only && sundir[2] <= 0.) {
642	<	continue; /* skipping nighttime points */
643	<	}
644	<	sun_ct = fdot(view_point, sundir) / ER;
640	>	if (epw->isWEA == WEAnot) {
641	>	fprintf(stderr, "EPW input\n");
642	>	} else if (epw->isWEA != WEAphotnorm) {
643	>	fprintf(stderr, "need WEA in photopic unit\n");
644	>	exit(1);
645	>	}
646
647	<	mtx_offset = NSSAMP * nskypatch * nstored;
648	<	nstored += 1;
649	<	/* make space for next row */
650	<	if (nstored > tstorage) {
651	<	tstorage += (tstorage >> 1) + nstored + 7;
652	<	mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
653	<	}
654	<	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);
647	>	while ((j = EPWread(epw, &erec)) > 0) {
648	>	const int          mo = erec.date.month+1;
649	>	const int          da = erec.date.day;
650	>	const double    hr = erec.date.hour;
651	>	double aod = erec.optdepth;
652	>	double cc = erec.skycover;
653	>	double            sda, sta, st;
654	>	int                      sun_in_sky;
655
656	<	clear_atmos.beta_m = mie_dp;
656	>	if (aod == 0.0) {
657	>	aod = AOD0_CA;
658	>	fprintf(stderr, "aod is zero, using default value %.3f\n", AOD0_CA);
659	>	}
660	>	/* compute solar position */
661	>	if ((mo == 2) & (da == 29)) {
662	>	julian_date = 60;
663	>	leap_day = 1;
664	>	} else
665	>	julian_date = jdate(mo, da) + leap_day;
666	>	sda = sdec(julian_date);
667	>	sta = stadj(julian_date);
668	>	st = hr + sta;
669	>	if (timeinterval > 0) {
670	>	if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120)
671	>	st = (st + timeinterval/120 + solar_sunrise(mo, da))/2;
672	>	else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120)
673	>	st = (st - timeinterval/120 + solar_sunset(mo, da))/2;
674	>	}
675	>	altitude = salt(sda, st);
676	>	sun_in_sky = (altitude > -deg_to_rad(SUN_ANG_DEG / 2.));
677
678	<	char gsdir[PATH_MAX];
581	<	size_t siz = strlen(ddir);
582	<	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);
678	>	azimuth = sazi(sda, st) + PI - deg_to_rad(rotation);
679
680	<	if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
681	<	getpath(clear_paths.scat, ".", R_OK) == NULL ||
682	<	getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
683	<	getpath(clear_paths.irrad, ".", R_OK) == NULL) {
684	<	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	>	vectorize(altitude, azimuth, sundir);
681	>	if (sun_hours_only && !sun_in_sky) {
682	>	continue; /* skipping nighttime points */
683	>	}
684	>	sun_ct = fdot(view_point, sundir) / ER;
685
686	<	DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
687	<	DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
604	<	DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
605	<	DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
686	>	dni = erec.dirillum;
687	>	dhi = erec.diffillum;
688
689	<	if (!solar_only)
690	<	compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
609	<	cc, mtx_data + mtx_offset);
610	<	if (!sky_only)
611	<	add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
612	<	cc, mtx_data + mtx_offset);
613	<	/* 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 */
689	>	mtx_offset = NSSAMP * nskypatch * nstored;
690	>	nstored += 1;
691
692	<	freedata(tau_clear_dp);
693	<	freedata(irrad_clear_dp);
694	<	freedata(scat_clear_dp);
695	<	freedata(scat1m_clear_dp);
696	<	}
697	<	freedata(mie_dp);
698	<	if (!ntsteps) {
699	<	fprintf(stderr, "%s: no valid time steps on input\n", progname);
700	<	exit(1);
701	<	}
702	<	/* check for junk at end */
703	<	while ((i = fgetc(stdin)) != EOF)
704	<	if (!isspace(i)) {
705	<	fprintf(stderr, "%s: warning - unexpected data past EOT: ", progname);
706	<	buf[0] = i;
707	<	buf[1] = '\0';
708	<	fgets(buf + 1, sizeof(buf) - 1, stdin);
709	<	fputs(buf, stderr);
710	<	fputc('\n', stderr);
711	<	break;
712	<	}
713	<	/* write out matrix */
714	<	if (outfmt != 'a')
715	<	SET_FILE_BINARY(stdout);
692	>	/* make space for next row */
693	>	if (nstored > tstorage) {
694	>	tstorage += (tstorage >> 1) + nstored + 7;
695	>	mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
696	>	}
697	>	ntsteps++; /* keep count of time steps */
698	>
699	>	/* compute sky patch values */
700	>	Atmosphere clear_atmos = init_atmos(aod, grefl);
701	>	int is_summer = (mo >= SUMMER_START && mo <= SUMMER_END);
702	>	if (s_latitude < 0) {
703	>	is_summer = !is_summer;
704	>	}
705	>	set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
706	>
707	>	clear_atmos.beta_m = mie_dp;
708	>
709	>	char gsdir[PATH_MAX];
710	>	size_t siz = strlen(ddir);
711	>	if (ISDIRSEP(ddir[siz - 1]))
712	>	ddir[siz - 1] = '\0';
713	>	snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
714	>	if (!make_directory(gsdir)) {
715	>	fprintf(stderr, "Failed creating atmos_data directory");
716	>	exit(1);
717	>	}
718	>	DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
719	>
720	>	if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
721	>	getpath(clear_paths.scat, ".", R_OK) == NULL ||
722	>	getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
723	>	getpath(clear_paths.irrad, ".", R_OK) == NULL) {
724	>	fprintf(stderr, "# Pre-computing...\n");
725	>	if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
726	>	fprintf(stderr, "Pre-compute failed\n");
727	>	return 0;
728	>	}
729	>	}
730	>
731	>	DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
732	>	DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
733	>	DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
734	>	DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
735	>
736	>	if (!solar_only)
737	>	compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
738	>	cc, dhi, view_point, mtx_data + mtx_offset);
739	>	if (!sky_only)
740	>	add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
741	>	cc, dni, mtx_data + mtx_offset);
742	>	/* monthly reporting */
743	>	if (verbose && mo != last_monthly)
744	>	fprintf(stderr, "%s: stepping through month %d...\n", progname,
745	>	last_monthly = mo);
746	>	}
747	>	if (j != EOF) {
748	>	fprintf(stderr, "%s: error on input\n", progname);
749	>	exit(1);
750	>	}
751	>	EPWclose(epw); epw = NULL;
752	>	freedata(mie_dp);
753	>	if (!ntsteps) {
754	>	fprintf(stderr, "%s: no valid time steps on input\n", progname);
755	>	exit(1);
756	>	}
757	>	/* write out matrix */
758	>	if (outfmt != 'a')
759	>	SET_FILE_BINARY(stdout);
760		#ifdef getc_unlocked
761	<	flockfile(stdout);
761	>	flockfile(stdout);
762		#endif
763	<	if (verbose)
764	<	fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
765	<	outfmt == 'a' ? "" : "binary ", nstored);
766	<	if (doheader) {
767	<	newheader("RADIANCE", stdout);
768	<	printargs(argc, argv, stdout);
769	<	printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude),
770	<	-RadToDeg(s_longitude));
771	<	printf("NROWS=%d\n", nskypatch);
772	<	printf("NCOLS=%d\n", nstored);
773	<	printf("NCOMP=%d\n", NSSAMP);
774	<	if ((outfmt == 'f') | (outfmt == 'd'))
775	<	fputendian(stdout);
776	<	fputformat((char *)getfmtname(outfmt), stdout);
777	<	putchar('\n');
778	<	}
779	<	/* patches are rows (outer sort) */
780	<	for (i = 0; i < nskypatch; i++) {
781	<	mtx_offset = NSSAMP * i;
782	<	switch (outfmt) {
783	<	case 'a':
784	<	for (j = 0; j < nstored; j++) {
785	<	for (int k = 0; k < NSSAMP; k++) {
786	<	printf("%.3g \n", mtx_data[mtx_offset + k]);
787	<	}
788	<	printf("\n");
789	<	mtx_offset += NSSAMP * nskypatch;
790	<	}
791	<	if (nstored > 1)
792	<	fputc('\n', stdout);
793	<	break;
794	<	case 'f':
795	<	for (j = 0; j < nstored; j++) {
796	<	putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
797	<	mtx_offset += NSSAMP * nskypatch;
798	<	}
799	<	break;
800	<	case 'd':
801	<	for (j = 0; j < nstored; j++) {
802	<	double ment[NSSAMP];
803	<	for (j = 0; j < NSSAMP; j++)
804	<	ment[j] = mtx_data[mtx_offset + j];
805	<	putbinary(ment, sizeof(double), NSSAMP, stdout);
806	<	mtx_offset += NSSAMP * nskypatch;
807	<	}
808	<	break;
809	<	}
810	<	if (ferror(stdout))
811	<	goto writerr;
812	<	}
813	<	alldone:
814	<	if (fflush(NULL) == EOF)
815	<	goto writerr;
816	<	if (verbose)
703	<	fprintf(stderr, "%s: done.\n", progname);
704	<	exit(0);
763	>	if (verbose)
764	>	fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
765	>	outfmt == 'a' ? "" : "binary ", nstored);
766	>	if (doheader) {
767	>	newheader("RADIANCE", stdout);
768	>	printargs(argc, argv, stdout);
769	>	printf("LATLONG= %.8f %.8f\n", rad_to_deg(s_latitude),
770	>	-rad_to_deg(s_longitude));
771	>	printf("NROWS=%d\n", nskypatch);
772	>	printf("NCOLS=%d\n", nstored);
773	>	printf("NCOMP=%d\n", NSSAMP);
774	>	float wvsplit[4] = {380, 480, 588, 780};
775	>	fputwlsplit(wvsplit, stdout);
776	>	if ((outfmt == 'f') | (outfmt == 'd'))
777	>	fputendian(stdout);
778	>	fputformat((char *)getfmtname(outfmt), stdout);
779	>	putchar('\n');
780	>	}
781	>	/* patches are rows (outer sort) */
782	>	for (i = 0; i < nskypatch; i++) {
783	>	mtx_offset = NSSAMP * i;
784	>	switch (outfmt) {
785	>	case 'a':
786	>	for (j = 0; j < nstored; j++) {
787	>	int k;
788	>	for (k = 0; k < NSSAMP; k++) {
789	>	printf("%.3g ", mtx_data[mtx_offset + k]);
790	>	}
791	>	printf("\n");
792	>	mtx_offset += NSSAMP * nskypatch;
793	>	}
794	>	if (nstored > 1)
795	>	fputc('\n', stdout);
796	>	break;
797	>	case 'f':
798	>	for (j = 0; j < nstored; j++) {
799	>	putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
800	>	mtx_offset += NSSAMP * nskypatch;
801	>	}
802	>	break;
803	>	case 'd':
804	>	for (j = 0; j < nstored; j++) {
805	>	double ment[NSSAMP];
806	>	for (j = 0; j < NSSAMP; j++)
807	>	ment[j] = mtx_data[mtx_offset + j];
808	>	putbinary(ment, sizeof(double), NSSAMP, stdout);
809	>	mtx_offset += NSSAMP * nskypatch;
810	>	}
811	>	break;
812	>	}
813	>	if (ferror(stdout))
814	>	goto writerr;
815	>	}
816	>	return 0;
817		userr:
818	<	fprintf(stderr,
819	<	"Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-r "
820	<	"deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
821	<	progname);
822	<	exit(1);
818	>	fprintf(stderr,
819	>	"Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-r "
820	>	"deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
821	>	progname);
822	>	exit(1);
823		fmterr:
824	<	fprintf(stderr, "%s: weather tape format error in header\n", progname);
825	<	exit(1);
824	>	fprintf(stderr, "%s: weather tape format error in header\n", progname);
825	>	exit(1);
826		writerr:
827	<	fprintf(stderr, "%s: write error on output\n", progname);
828	<	exit(1);
827	>	fprintf(stderr, "%s: write error on output\n", progname);
828	>	exit(1);
829		}

Diff Legend

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