ViewVC Help

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

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

Comparing ray/src/gen/gensdaymtx.c (file contents):
Revision 1.5 by greg, Thu Apr 10 23:30:58 2025 UTC vs.
Revision 1.13 by greg, Fri Aug 22 01:04:30 2025 UTC

#	Line 21 | Line 21 | static const char RCSid[] = "$Id$";
21		#include "sun.h"
22		#include "loadEPW.h"
23
24	–	#ifndef M_PI
25	–	#define M_PI 3.14159265358579
26	–	#endif
24
25	<	#define vector(v, alt, azi)                                                    \
26	<	((v)[1] = cos(alt), (v)[0] = (v)[1] * sin(azi), (v)[1] *= cos(azi),          \
27	<	(v)[2] = sin(alt))
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	<	#define rh_vector(v, i) vector(v, rh_palt[i], rh_pazi[i])
33	>	/* Mean normalized relative daylight spectra where CCT = 6415K for overcast */
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
39	>	};
40
41	<	#define rh_cos(i) tsin(rh_palt[i])
41	>	enum {
42	>	NSUNPATCH = 4          /* max. # patches to spread sun into */
43	>	};
44
45	<	#define solar_minute(jd, hr) ((24 * 60) * ((jd) - 1) + (int)((hr) * 60. + .5))
46	<
47	<
48	<	char *progname;
49	<
50	<	double altitude;   /* Solar altitude (radians) */
51	<	double azimuth;    /* Solar azimuth (radians) */
52	<	int julian_date;   /* Julian date */
44	<	double sun_zenith; /* Sun zenith angle (radians) */
45	<	int input = 0;     /* Input type */
46	<	int output = 0;    /* Output type */
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	<	const double ARCTIC_LAT = 67.;
57	<	const double TROPIC_LAT = 23.;
58	<	const int SUMMER_START = 4;
59	<	const int SUMMER_END = 9;
60	<	const double GNORM = 0.777778;
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
55	–	const double D65EFF = 203.; /* standard illuminant D65 */
67
68	<	/* Mean normalized relative daylight spectra where CCT = 6415K for overcast */
69	<	const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
70	<	1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
71	<	1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
72	<	0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
73	<	/* Degrees into radians */
74	<	#define DegToRad(deg) ((deg) * (PI / 180.))
68	>	static inline double
69	>	deg_to_rad
70	>	(
71	>	double deg
72	>	)
73	>	{
74	>	return deg * (PI / 180.);
75	>	}
76
77	<	/* Radiuans into degrees */
78	<	#define RadToDeg(rad) ((rad) * (180. / PI))
77	>	static inline double
78	>	rad_to_deg
79	>	(
80	>	double rad
81	>	)
82	>	{
83	>	return rad * (180. / PI);
84	>	}
85
86	<	#ifndef NSUNPATCH
87	<	#define NSUNPATCH 4 /* max. # patches to spread sun into */
88	<	#endif
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);
96	>	v[1] *= cos(azimuth);
97	>	v[2] = sin(altitude);
98	>	}
99
100	<	#define SUN_ANG_DEG 0.533 /* sun full-angle in degrees */
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	<	int nsuns = NSUNPATCH;    /* number of sun patches to use */
112	<	double fixed_sun_sa = -1; /* fixed solid angle per sun? */
113	<
114	<	int verbose = 0; /* progress reports to stderr? */
115	<
116	<	int outfmt = 'a'; /* output format */
117	<
118	<	int rhsubdiv = 1; /* Reinhart sky subdivisions */
119	<
120	<	COLOR skycolor = {.96, 1.004, 1.118}; /* sky coloration */
84	<	COLOR suncolor = {1., 1., 1.};        /* sun color */
85	<	double grefl = .2;                    /* ground reflectance */
86	<
87	<	int nskypatch;  /* number of Reinhart patches */
88	<	float *rh_palt; /* sky patch altitudes (radians) */
89	<	float *rh_pazi; /* sky patch azimuths (radians) */
90	<	float *rh_dom;  /* sky patch solid angle (sr) */
91	<
92	<	double sun_ct;
93	<
94	<	inline void vectorize(double altitude, double azimuth, FVECT v) {
95	<	v[1] = cos(altitude);
96	<	v[0] = (v)[1] * sin(azimuth);
97	<	v[1] *= cos(azimuth);
98	<	v[2] = sin(altitude);
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;
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;
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) {
145	<	switch (fmt) {
146	<	case 'a':
147	<	return ("ascii");
148	<	case 'f':
149	<	return ("float");
150	<	case 'd':
151	<	return ("double");
152	<	}
124	<	return ("unknown");
125	<	}
144	>	static const char *
145	>	getfmtname
146	>	(
147	>	int fmt
148	>	)
149	>	{
150	>	switch (fmt) {
151	>	case 'a':
152	>	return ("ascii");
153
154	<	static inline double wmean2(const double a, const double b, const double x) {
155	<	return a * (1 - x) + b * x;
129	<	}
154	>	case 'f':
155	>	return ("float");
156
157	<	static inline double wmean(const double a, const double x, const double b,
158	<	const double y) {
159	<	return (a * x + b * y) / (a + b);
157	>	case 'd':
158	>	return ("double");
159	>	}
160	>	return ("unknown");
161		}
162
163
164	<	static double get_overcast_zenith_brightness(const double sundir[3]) {
165	<	double zenithbr;
166	<	if (sundir[2] < 0) {
167	<	zenithbr = 0;
168	<	} else {
169	<	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
170	<	}
171	<	return zenithbr;
164	>	static double
165	>	get_overcast_zenith_brightness
166	>	(
167	>	const double sundir[3]
168	>	)
169	>	{
170	>	double zenithbr;
171	>	if (sundir[2] < 0) {
172	>	zenithbr = 0;
173	>	} else {
174	>	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFFICACY;
175	>	}
176	>	return zenithbr;
177		}
178
179	+
180		/* from gensky.c */
181	<	static double get_overcast_brightness(const double dz, const double zenithbr) {
182	<	double groundbr = zenithbr * GNORM;
183	<	return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
184	<	pow(dz + 1.01, -10), groundbr);
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);
192		}
193
194		double
195	<	solar_sunset(int month, int day)
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))));
203	<	return(12 + (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15));
203	>	return 12 + (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
204		}
205
206
207		double
208	<	solar_sunrise(int month, int day)
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 - (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15));
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};
227	<	const double alpha = (PI / 2.) / (NROW * rhsubdiv + .5);
228	<	int p, i, j;
229	<	/* allocate patch angle arrays */
230	<	nskypatch = 0;
231	<	for (p = 0; p < NROW; p++)
232	<	nskypatch += tnaz[p];
233	<	nskypatch *= rhsubdiv * rhsubdiv;
234	<	nskypatch += 2;
235	<	rh_palt = (float *)malloc(sizeof(float) * nskypatch);
236	<	rh_pazi = (float *)malloc(sizeof(float) * nskypatch);
237	<	rh_dom = (float *)malloc(sizeof(float) * nskypatch);
238	<	if ((rh_palt == NULL) | (rh_pazi == NULL) | (rh_dom == NULL)) {
239	<	fprintf(stderr, "%s: out of memory in rh_init()\n", progname);
240	<	exit(1);
241	<	}
242	<	rh_palt[0] = -PI / 2.; /* ground & zenith patches */
243	<	rh_pazi[0] = 0.;
244	<	rh_dom[0] = 2. * PI;
245	<	rh_palt[nskypatch - 1] = PI / 2.;
246	<	rh_pazi[nskypatch - 1] = 0.;
247	<	rh_dom[nskypatch - 1] = 2. * PI * (1. - cos(alpha * .5));
248	<	p = 1; /* "normal" patches */
249	<	for (i = 0; i < NROW * rhsubdiv; i++) {
250	<	const float ralt = alpha * (i + .5);
251	<	const int ninrow = tnaz[i / rhsubdiv] * rhsubdiv;
252	<	const float dom =
253	<	2. * PI * (sin(alpha * (i + 1)) - sin(alpha * i)) / (double)ninrow;
254	<	for (j = 0; j < ninrow; j++) {
255	<	rh_palt[p] = ralt;
256	<	rh_pazi[p] = 2. * PI * j / (double)ninrow;
257	<	rh_dom[p++] = dom;
258	<	}
259	<	}
260	<	return nskypatch;
226	>	static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6};
227	>	const double alpha = (PI / 2.) / (NROW * rhsubdiv + .5);
228	>	int p, i, j;
229	>	/* allocate patch angle arrays */
230	>	nskypatch = 0;
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);
237	>	rh_pazi = (float *)malloc(sizeof(float) * nskypatch);
238	>	rh_dom = (float *)malloc(sizeof(float) * nskypatch);
239	>	if ((rh_palt == NULL) | (rh_pazi == NULL) | (rh_dom == NULL)) {
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 */
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 */
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;
255	>	for (j = 0; j < ninrow; j++) {
256	>	rh_palt[p] = ralt;
257	>	rh_pazi[p] = 2. * PI * j / (double)ninrow;
258	>	rh_dom[p++] = dom;
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) {
268	<	if (mtx_data == NULL)
269	<	mtx_data = (float *)malloc(sizeof(float) * NSSAMP * nsteps * npatch);
270	<	else
271	<	mtx_data =
272	<	(float *)realloc(mtx_data, sizeof(float) * NSSAMP * nsteps * npatch);
273	<	if (mtx_data == NULL) {
274	<	fprintf(stderr, "%s: out of memory in resize_dmatrix(%d,%d)\n", progname,
275	<	nsteps, npatch);
276	<	exit(1);
277	<	}
278	<	return (mtx_data);
267	>	float *
268	>	resize_dmatrix
269	>	(
270	>	float *mtx_data,
271	>	int nsteps,
272	>	int npatch
273	>	)
274	>	{
275	>	if (mtx_data == NULL) {
276	>	mtx_data = (float * ) malloc(sizeof(float) * NSSAMP * nsteps * npatch);
277	>	}else{
278	>	mtx_data = (float * ) realloc(mtx_data,
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;
288		}
289
290	<	static Atmosphere init_atmos(const double aod, const double grefl) {
291	<	Atmosphere atmos = {.ozone_density = {.layers =
292	<	{
293	<	{.width = 25000.0,
294	<	.exp_term = 0.0,
295	<	.exp_scale = 0.0,
296	<	.linear_term = 1.0 / 15000.0,
297	<	.constant_term = -2.0 / 3.0},
298	<	{.width = AH,
299	<	.exp_term = 0.0,
300	<	.exp_scale = 0.0,
301	<	.linear_term = -1.0 / 15000.0,
302	<	.constant_term = 8.0 / 3.0},
303	<	}},
304	<	.rayleigh_density = {.layers =
305	<	{
306	<	{.width = AH,
307	<	.exp_term = 1.0,
308	<	.exp_scale = -1.0 / HR_MS,
309	<	.linear_term = 0.0,
310	<	.constant_term = 0.0},
311	<	}},
312	<	.beta_r0 = BR0_MS,
313	<	.beta_scale = aod / AOD0_CA,
314	<	.beta_m = NULL,
315	<	.grefl = grefl};
316	<	return atmos;
290	>	static Atmosphere
291	>	init_atmos
292	>	(
293	>	const double aod,
294	>	const double grefl
295	>	)
296	>	{
297	>	Atmosphere atmos = {
298	>	.ozone_density = {
299	>	.layers = {
300	>	{
301	>	.width = 25000.0,
302	>	.exp_term = 0.0,
303	>	.exp_scale = 0.0,
304	>	.linear_term = 1.0 / 15000.0,
305	>	.constant_term = -2.0 / 3.0
306	>	},
307	>	{
308	>	.width = AH,
309	>	.exp_term = 0.0,
310	>	.exp_scale = 0.0,
311	>	.linear_term = -1.0 / 15000.0,
312	>	.constant_term = 8.0 / 3.0
313	>	},
314	>	}
315	>	},
316	>	.rayleigh_density = {
317	>	.layers = {
318	>	{
319	>	.width = AH,
320	>	.exp_term = 1.0,
321	>	.exp_scale = -1.0 / HR_MS,
322	>	.linear_term = 0.0,
323	>	.constant_term = 0.0
324	>	},
325	>	}
326	>	},
327	>	.beta_r0 = BR0_MS,
328	>	.beta_scale = aod / AOD0_CA,
329	>	.beta_m = NULL,
330	>	.grefl = grefl
331	>	};
332	>	return atmos;
333		}
334
335	<	static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
336	<	const char *tag) {
337	<	DpPaths paths;
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", dir, DIRSEP, tag,
347	<	mname, aod);
348	<	snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
349	<	mname, aod);
350	<	snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
351	<	tag, mname, aod);
352	<	snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
353	<	mname, aod);
346	>	snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat",
347	>	dir, DIRSEP, tag, mname, aod);
348	>	snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat",
349	>	dir, DIRSEP, tag, mname, aod);
350	>	snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat",
351	>	dir, DIRSEP, tag, mname, aod);
352	>	snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat",
353	>	dir, DIRSEP, tag, mname, aod);
354
355	<	return paths;
355	>	return paths;
356		}
357
358	<	static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
359	<	const int is_summer,
360	<	const double s_latitude) {
361	<	/* Set rayleigh density profile */
362	<	if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
363	<	tag[0] = 's';
364	<	if (is_summer) {
365	<	tag[1] = 's';
366	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
367	<	atmos->beta_r0 = BR0_SS;
368	<	} else {
369	<	tag[1] = 'w';
370	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
371	<	atmos->beta_r0 = BR0_SW;
372	<	}
373	<	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
374	<	tag[0] = 'm';
375	<	if (is_summer) {
376	<	tag[1] = 's';
377	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
378	<	atmos->beta_r0 = BR0_MS;
379	<	} else {
380	<	tag[1] = 'w';
381	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
382	<	atmos->beta_r0 = BR0_MW;
383	<	}
384	<	} else {
385	<	tag[0] = 't';
386	<	tag[1] = 'r';
387	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
388	<	atmos->beta_r0 = BR0_T;
389	<	}
390	<	tag[2] = '\0';
358	>
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) {
370	>	tag[0] = 's';
371	>	if (is_summer) {
372	>	tag[1] = 's';
373	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
374	>	atmos->beta_r0 = BR0_SS;
375	>	} else {
376	>	tag[1] = 'w';
377	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
378	>	atmos->beta_r0 = BR0_SW;
379	>	}
380	>	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
381	>	tag[0] = 'm';
382	>	if (is_summer) {
383	>	tag[1] = 's';
384	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
385	>	atmos->beta_r0 = BR0_MS;
386	>	} else {
387	>	tag[1] = 'w';
388	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
389	>	atmos->beta_r0 = BR0_MW;
390	>	}
391	>	} else {
392	>	tag[0] = 't';
393	>	tag[1] = 'r';
394	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
395	>	atmos->beta_r0 = BR0_T;
396	>	}
397	>	tag[2] = '\0';
398		}
399
400	+
401		/* Add in solar direct to nearest sky patches (GW) */
402	<	void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
403	<	DATARRAY *irrad, double ccover, double dirnorm, float *parr) {
404	<	FVECT svec;
405	<	double near_dprod[NSUNPATCH];
406	<	int near_patch[NSUNPATCH];
407	<	double wta[NSUNPATCH], wtot;
408	<	int i, j, p;
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];
416	>	int near_patch[NSUNPATCH];
417	>	double wta[NSUNPATCH], wtot;
418	>	int i, j, p;
419
420	<	/* identify nsuns closest patches */
421	<	for (i = nsuns; i--;)
422	<	near_dprod[i] = -1.;
423	<	vectorize(altitude, azimuth, svec);
424	<	for (p = 1; p < nskypatch; p++) {
425	<	FVECT pvec;
426	<	double dprod;
427	<	vectorize(rh_palt[p], rh_pazi[p], pvec);
428	<	dprod = DOT(pvec, svec);
429	<	for (i = 0; i < nsuns; i++)
430	<	if (dprod > near_dprod[i]) {
431	<	for (j = nsuns; --j > i;) {
432	<	near_dprod[j] = near_dprod[j - 1];
433	<	near_patch[j] = near_patch[j - 1];
434	<	}
435	<	near_dprod[i] = dprod;
436	<	near_patch[i] = p;
437	<	break;
438	<	}
439	<	}
440	<	/* Get solar radiance */
441	<	double sun_radiance[NSSAMP] = {0};
442	<	get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
443	<	if (ccover > 0) {
444	<	double zenithbr = get_overcast_zenith_brightness(sundir);
445	<	double skybr = get_overcast_brightness(sundir[2], zenithbr);
446	<	int l;
447	<	for (l = 0; l < NSSAMP; ++l) {
448	<	sun_radiance[l] = wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
449	<	}
450	<	}
451	<	/* Normalize */
452	<	double sum = 0.0;
453	<	for (i = 0; i < NSSAMP; ++i) {
454	<	sum += sun_radiance[i];
455	<	}
456	<	double mean = sum / NSSAMP;
420	>	/* identify nsuns closest patches */
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++) {
431	>	if (dprod > near_dprod[i]) {
432	>	for (j = nsuns; --j > i;) {
433	>	near_dprod[j] = near_dprod[j - 1];
434	>	near_patch[j] = near_patch[j - 1];
435	>	}
436	>	near_dprod[i] = dprod;
437	>	near_patch[i] = p;
438	>	break;
439	>	}
440	>	}
441	>	}
442	>	/* Get solar radiance */
443	>	double sun_radiance[NSSAMP] = {0};
444	>	get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
445	>	if (ccover > 0) {
446	>	double zenithbr = get_overcast_zenith_brightness(sundir);
447	>	double skybr = get_overcast_brightness(sundir[2], zenithbr);
448	>	int l;
449	>	for (l = 0; l < NSSAMP; ++l) {
450	>	sun_radiance[l] = wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
451	>	}
452	>	}
453	>	/* Normalize */
454	>	double sum = 0.0;
455	>	for (i = 0; i < NSSAMP; ++i) {
456	>	sum += sun_radiance[i];
457	>	}
458	>	double mean = sum / NSSAMP;
459
460	<	double intensity = mean * WVLSPAN;
461	<	if (dirnorm > 0) {
462	<	intensity = dirnorm / SOLOMG / WHTEFFICACY;
463	<	}
464	<	double dir_ratio = 1.;
465	<	if (mean > 0)
466	<	dir_ratio = intensity / mean;
467	<	for (i = 0; i < NSSAMP; ++i) {
468	<	sun_radiance[i] *= dir_ratio;
469	<	}
460	>	double intensity = mean * WVLSPAN;
461	>	if (dirnorm >= 0) {
462	>	intensity = dirnorm / SOLOMG / WHTEFFICACY;
463	>	}
464	>	double dir_ratio = 1.;
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--;)
475	<	wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
476	<	/* add to nearest patch radiances */
477	<	for (i = nsuns; i--;) {
478	<	float *pdest = parr + NSSAMP * near_patch[i];
479	<	int k;
480	<	for (k = 0; k < NSSAMP; k++) {
481	<	*pdest++ = sun_radiance[k] * wta[i] / wtot;
482	<	}
483	<	}
472	>	/* weight by proximity */
473	>	wtot = 0;
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];
480	>	int k;
481	>	for (k = 0; k < NSSAMP; k++) {
482	>	*pdest++ += sun_radiance[k] * wta[i] / wtot;
483	>	}
484	>	}
485		}
486
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};
487
488	<	get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
489	<	for (k = 0; k < NSSAMP; ++k) {
490	<	sky_radiance[k] *= WVLSPAN;
491	<	}
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;
504	>	for (i = 1; i < nskypatch; i++) {
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	<	if (ccover > 0) {
511	<	double skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr);
512	<	int k;
513	<	for (k = 0; k < NSSAMP; ++k) {
514	<	sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
515	<	}
403	<	}
510	>	SCOLOR sky_radiance = {0};
511	>	get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
512	>	int k;
513	>	for (k = 0; k < NSSAMP; ++k) {
514	>	sky_radiance[k] *= WVLSPAN;
515	>	}
516
517	<	/* calibration */
518	<	for (k = 0; k < NSSAMP; ++k) {
519	<	sky_radiance[k] *= dif_ratio;
520	<	}
517	>	if (ccover > 0) {
518	>	double skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr);
519	>	for (k = 0; k < NSSAMP; ++k) {
520	>	sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
521	>	}
522	>	}
523
524	<	for (k = 0; k < NSSAMP; ++k) {
525	<	parr[NSSAMP * i + k] = sky_radiance[k];
526	<	}
527	<	}
524	>	/* calibration */
525	>	for (k = 0; k < NSSAMP; ++k) {
526	>	sky_radiance[k] *= dif_ratio;
527	>	}
528	>
529	>	for (k = 0; k < NSSAMP; ++k) {
530	>	parr[NSSAMP * i + k] = sky_radiance[k];
531	>	}
532	>	}
533		}
534
416	–	/* Return maximum of two doubles */
417	–	static inline double dmax(double a, double b) { return (a > b) ? a : b; }
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, float *parr) {
539	<	float sun_zenith;
540	<	SCOLOR sky_radiance = {0};
541	<	SCOLOR ground_radiance = {0};
542	<	SCOLR sky_sclr = {0};
543	<	SCOLR ground_sclr = {0};
544	<	FVECT view_point = {0, 0, ER};
545	<	const double radius = VLEN(view_point);
546	<	const double sun_ct = fdot(view_point, sundir) / radius;
547	<	const FVECT rdir_grnd = {0, 0, -1};
548	<	const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
549	<	const double nu_grnd = fdot(rdir_grnd, sundir);
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};
552	>	SCOLOR ground_radiance = {0};
553	>	SCOLR sky_sclr = {0};
554	>	SCOLR ground_sclr = {0};
555	>	const double radius = VLEN(view_point);
556	>	const double sun_ct = fdot(view_point, sundir) / radius;
557	>	const FVECT rdir_grnd = {0, 0, -1};
558	>	const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
559	>	const double nu_grnd = fdot(rdir_grnd, sundir);
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)
564	<	sun_zenith = DegToRad(90.0);
565	<	else if (altitude >= DegToRad(87.0))
566	<	sun_zenith = DegToRad(3.0);
567	<	else
568	<	sun_zenith = DegToRad(90.0) - altitude;
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) {
564	>	sun_zenith = deg_to_rad(90.0);
565	>	}else if (altitude >= deg_to_rad(87.0)) {
566	>	sun_zenith = deg_to_rad(3.0);
567	>	}else{
568	>	sun_zenith = deg_to_rad(90.0) - altitude;
569	>	}
570
571	<	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
571	>	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
572
573	<	/* diffuse calibration factor */
574	<	double dif_ratio = 1;
575	<	if (difhor > 0) {
576	<	DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad, radius, sun_ct);
577	<	double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
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 */
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 */
587	<	}
588	<	}
573	>	/* diffuse calibration factor */
574	>	double dif_ratio = 1;
575	>	if (difhor >= 0) {
576	>	DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad, radius, sun_ct);
577	>	double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
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 */
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 */
587	>	}
588	>	}
589
590	<	/* Compute ground radiance (include solar contribution if any) */
591	<	get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
592	<	mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
593	<	int j;
594	<	for (j = 0; j < NSSAMP; j++) {
595	<	parr[j] *= WVLSPAN;
596	<	}
597	<	calc_sky_patch_radiance(scat, scat1m, irrad, ccover, dif_ratio, overcast_zenithbr, parr);
590	>	/* Compute ground radiance (include solar contribution if any) */
591	>	get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
592	>	mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
593	>	int j;
594	>	for (j = 0; j < NSSAMP; j++) {
595	>	parr[j] *= WVLSPAN;
596	>	}
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 = -1.0;                                                 /* direct normal illuminance */
620	>	double dhi = -1.0;                                                 /* diffuse horizontal illuminance */
621
622	<	EPWheader   *epw = NULL;    /* EPW/WEA input file */
623	<	EPWrecord   erec;           /* current EPW/WEA input record */
624	<	int         doheader = 1; /* output header? */
625	<	double      rotation = 0.0;
626	<	double      elevation = 0;
627	<	int         leap_day = 0;       /* add leap day? */
628	<	int         sun_hours_only = 0; /* only output sun hours? */
629	<	int         dir_is_horiz;           /* direct is meas. on horizontal? */
630	<	float       *mtx_data = NULL;
631	<	int         ntsteps = 0;      /* number of time steps */
632	<	int         tstorage = 0;     /* number of allocated time steps */
633	<	int         nstored = 0;      /* number of time steps in matrix */
634	<	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;
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	<	progname = argv[0];
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 */
641	<	solar_only = 1;
642	<	break;
643	<	case 's': /* sky only (no direct) */
644	<	sky_only = 1;
645	<	break;
646	<	case 'g':
647	<	grefl = atof(argv[++i]);
648	<	break;
649	<	case 'm':
650	<	rhsubdiv = atoi(argv[++i]);
651	<	break;
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')
657	<	goto userr;
658	<	rotation = atof(argv[++i]);
659	<	break;
660	<	case 'u': /* solar hours only */
661	<	sun_hours_only = 1;
662	<	break;
663	<	case 'p':
664	<	ddir = argv[++i];
665	<	break;
666	<	case 'v': /* verbose progress reports */
667	<	verbose++;
668	<	break;
669	<	case 'h': /* turn off header */
670	<	doheader = 0;
671	<	break;
672	<	case '5': /* 5-phase calculation */
673	<	nsuns = 1;
674	<	fixed_sun_sa = PI / 360. * atof(argv[++i]);
675	<	if (fixed_sun_sa <= 0) {
676	<	fprintf(
677	<	stderr,
678	<	"%s: missing solar disk size argument for '-5' option\n",
679	<	progname);
680	<	exit(1);
681	<	}
682	<	fixed_sun_sa *= fixed_sun_sa * PI;
683	<	break;
684	<	case 'i':
685	<	timeinterval = atof(argv[++i]);
686	<	break;
687	<	case 'o': /* output format */
688	<	switch (argv[i][2]) {
689	<	case 'f':
690	<	case 'd':
691	<	case 'a':
692	<	outfmt = argv[i][2];
693	<	break;
694	<	default:
695	<	goto userr;
696	<	}
697	<	break;
698	<	default:
699	<	goto userr;
700	<	}
701	<	}
702	<	if (i < argc - 1)
703	<	goto userr;
638	>	for (i = 1; i < argc && argv[i][0] == '-'; i++) {
639	>	switch (argv[i][1]) {
640	>	case 'd':         /* solar (direct) only */
641	>	solar_only = 1;
642	>	break;
643	>	case 's':         /* sky only (no direct) */
644	>	sky_only = 1;
645	>	break;
646	>	case 'g':
647	>	grefl = atof(argv[++i]);
648	>	break;
649	>	case 'm':
650	>	rhsubdiv = atoi(argv[++i]);
651	>	break;
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') {
657	>	goto userr;
658	>	}
659	>	rotation = atof(argv[++i]);
660	>	break;
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 */
668	>	verbose++;
669	>	break;
670	>	case 'h':         /* turn off header */
671	>	doheader = 0;
672	>	break;
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,
679	>	"%s: missing solar disk size argument for '-5' option\n",
680	>	progname);
681	>	exit(1);
682	>	}
683	>	fixed_sun_sa *= fixed_sun_sa * PI;
684	>	break;
685	>	case 'i':
686	>	timeinterval = atof(argv[++i]);
687	>	break;
688	>	case 'o':         /* output format */
689	>	switch (argv[i][2]) {
690	>	case 'f':
691	>	case 'd':
692	>	case 'a':
693	>	outfmt = argv[i][2];
694	>	break;
695	>	default:
696	>	goto userr;
697	>	}
698	>	break;
699	>	default:
700	>	goto userr;
701	>	}
702	>	}
703	>	if (i < argc - 1) {
704	>	goto userr;
705	>	}
706
707	<	epw = EPWopen(argv[i]);
708	<	if (epw == NULL)
709	<	exit(1);
710	<	if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
711	<	fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
712	<	exit(1);
713	<	}
714	<	if (verbose) {
715	<	if (i == argc - 1)
716	<	fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
717	<	else
718	<	fprintf(stderr, "%s: reading weather tape from <stdin>\n", progname);
719	<	}
720	<	s_latitude = epw->loc.latitude;
721	<	s_longitude = -epw->loc.longitude;
722	<	s_meridian = -15.*epw->loc.timezone;
723	<	elevation = epw->loc.elevation;
724	<	switch (epw->isWEA) {               /* translate units */
725	<	case WEAnot:
726	<	case WEAradnorm:
727	<	input = 1;              /* radiometric quantities */
728	<	dir_is_horiz = 0;    /* direct is perpendicular meas. */
729	<	break;
730	<	case WEAradhoriz:
731	<	input = 1;        /* radiometric quantities */
732	<	dir_is_horiz = 1;    /* solar measured horizontally */
733	<	break;
734	<	case WEAphotnorm:
735	<	input = 2;        /* photometric quantities */
736	<	dir_is_horiz = 0;    /* direct is perpendicular meas. */
737	<	break;
738	<	default:
739	<	goto fmterr;
740	<	}
707	>	epw = EPWopen(argv[i]);
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) {
717	>	fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
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 */
727	>	case WEAnot:
728	>	case WEAradnorm:
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 */
735	>	break;
736	>	case WEAphotnorm:
737	>	input = 2;                      /* photometric quantities */
738	>	dir_is_horiz = 0;               /* direct is perpendicular meas. */
739	>	break;
740	>	default:
741	>	goto fmterr;
742	>	}
743
744	<	rh_init();
744	>	rh_init();
745
746	<	if (verbose) {
747	<	fprintf(stderr, "%s: location '%s %s'\n", progname, epw->loc.city, epw->loc.country);
748	<	fprintf(
749	<	stderr,
750	<	"%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
751	<	progname, s_latitude, s_longitude);
752	<	if (rotation != 0)
753	<	fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
754	<	}
746	>	if (verbose) {
747	>	fprintf(stderr, "%s: location '%s %s'\n", progname, epw->loc.city, epw->loc.country);
748	>	fprintf(
749	>	stderr,
750	>	"%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
751	>	progname, s_latitude, s_longitude);
752	>	if (rotation != 0) {
753	>	fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
754	>	}
755	>	}
756
757	<	s_latitude = DegToRad(s_latitude);
758	<	s_longitude = DegToRad(s_longitude);
759	<	s_meridian = DegToRad(s_meridian);
760	<	/* initial allocation */
761	<	mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
757	>	s_latitude = deg_to_rad(s_latitude);
758	>	s_longitude = deg_to_rad(s_longitude);
759	>	s_meridian = deg_to_rad(s_meridian);
760	>	/* initial allocation */
761	>	mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
762
763	<	/* Load mie density data */
764	<	DATARRAY *mie_dp = getdata(mie_path);
765	<	if (mie_dp == NULL) {
766	<	fprintf(stderr, "Error reading mie data\n");
767	<	return 0;
768	<	}
763	>	/* Load mie density data */
764	>	DATARRAY *mie_dp = getdata(mie_path);
765	>	if (mie_dp == NULL) {
766	>	fprintf(stderr, "Error reading mie data\n");
767	>	return 1;
768	>	}
769
770	<	if (epw->isWEA == WEAnot) {
771	<	fprintf(stderr, "EPW input\n");
772	<	} else if (epw->isWEA != WEAphotnorm) {
773	<	fprintf(stderr, "need WEA in photopic unit\n");
774	<	exit(1);
775	<	}
770	>	if (epw->isWEA == WEAnot) {
771	>	fprintf(stderr, "EPW input\n");
772	>	} else if (epw->isWEA != WEAphotnorm) {
773	>	fprintf(stderr, "need WEA in photopic unit\n");
774	>	exit(1);
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;
782	<	double cc = erec.skycover;
783	<	double          sda, sta, st;
784	<	int             sun_in_sky;
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 * M_PER_KM;
782	>	if (aod >= 999.0) {
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) {
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;
793
794	<	if (aod == 0.0) {
795	<	aod = AOD0_CA;
796	<	fprintf(stderr, "aod is zero, using default value %.3f\n", AOD0_CA);
797	<	}
798	<	/* compute solar position */
799	<	if ((mo == 2) & (da == 29)) {
800	<	julian_date = 60;
801	<	leap_day = 1;
802	<	} else
803	<	julian_date = jdate(mo, da) + leap_day;
656	<	sda = sdec(julian_date);
657	<	sta = stadj(julian_date);
658	<	st = hr + sta;
794	>	/* compute solar position */
795	>	if ((mo == 2) & (da == 29)) {
796	>	julian_date = 60;
797	>	leap_day = 1;
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 > -DegToRad(SUN_ANG_DEG / 2.));
811	>	altitude = salt(sda, st);
812	>	sun_in_sky = (altitude > -deg_to_rad(SUN_ANG_DEG / 2.));
813
814	<	azimuth = sazi(sda, st) + PI - DegToRad(rotation);
814	>	azimuth = sazi(sda, st) + PI - deg_to_rad(rotation);
815
816	<	vectorize(altitude, azimuth, sundir);
817	<	if (sun_hours_only && !sun_in_sky) {
818	<	continue; /* skipping nighttime points */
819	<	}
820	<	sun_ct = fdot(view_point, sundir) / ER;
816	>	vectorize(altitude, azimuth, sundir);
817	>	if (sun_hours_only && !sun_in_sky) {
818	>	continue;             /* skipping nighttime points */
819	>	}
820	>	sun_ct = fdot(view_point, sundir) / ER;
821
822	<	dni = erec.dirillum;
823	<	dhi = erec.diffillum;
678	<	printf("%d %d %f %f %f %f %f\n", mo, da, hr, dni, dhi, aod, cc);
822	>	dni = erec.dirillum;
823	>	dhi = erec.diffillum;
824
825	<	mtx_offset = NSSAMP * nskypatch * nstored;
826	<	nstored += 1;
825	>	mtx_offset = NSSAMP * nskypatch * nstored;
826	>	nstored += 1;
827
828	<	/* make space for next row */
829	<	if (nstored > tstorage) {
830	<	tstorage += (tstorage >> 1) + nstored + 7;
831	<	mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
832	<	}
833	<	ntsteps++; /* keep count of time steps */
828	>	/* make space for next row */
829	>	if (nstored > tstorage) {
830	>	tstorage += (tstorage >> 1) + nstored + 7;
831	>	mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
832	>	}
833	>	ntsteps++;         /* keep count of time steps */
834
835	<	/* compute sky patch values */
836	<	Atmosphere clear_atmos = init_atmos(aod, grefl);
837	<	int is_summer = (mo >= SUMMER_START && mo <= SUMMER_END);
838	<	if (s_latitude < 0) {
839	<	is_summer = !is_summer;
840	<	}
841	<	set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
835	>	/* compute sky patch values */
836	>	Atmosphere clear_atmos = init_atmos(aod, grefl);
837	>	int is_summer = (mo >= SUMMER_START && mo <= SUMMER_END);
838	>	if (s_latitude < 0) {
839	>	is_summer = !is_summer;
840	>	}
841	>	set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
842
843	<	clear_atmos.beta_m = mie_dp;
843	>	clear_atmos.beta_m = mie_dp;
844
845	<	char gsdir[PATH_MAX];
846	<	size_t siz = strlen(ddir);
847	<	if (ISDIRSEP(ddir[siz - 1]))
848	<	ddir[siz - 1] = '\0';
849	<	snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
850	<	if (!make_directory(gsdir)) {
851	<	fprintf(stderr, "Failed creating atmos_data directory");
852	<	exit(1);
853	<	}
854	<	DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
845	>	char gsdir[PATH_MAX];
846	>	size_t siz = strlen(ddir);
847	>	if (ISDIRSEP(ddir[siz - 1])) {
848	>	ddir[siz - 1] = '\0';
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);
854	>	}
855	>	DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
856
857	<	if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
858	<	getpath(clear_paths.scat, ".", R_OK) == NULL ||
859	<	getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
860	<	getpath(clear_paths.irrad, ".", R_OK) == NULL) {
861	<	fprintf(stderr, "# Pre-computing...\n");
862	<	if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
863	<	fprintf(stderr, "Pre-compute failed\n");
864	<	return 0;
865	<	}
866	<	}
857	>	if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
858	>	getpath(clear_paths.scat, ".", R_OK) == NULL ||
859	>	getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
860	>	getpath(clear_paths.irrad, ".", R_OK) == NULL) {
861	>	fprintf(stderr, "# Pre-computing...\n");
862	>	if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
863	>	fprintf(stderr, "Pre-compute failed\n");
864	>	return 1;
865	>	}
866	>	}
867
868	<	DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
869	<	DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
870	<	DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
871	<	DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
868	>	DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
869	>	DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
870	>	DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
871	>	DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
872
873	<	if (!solar_only)
874	<	compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
875	<	cc, dhi, mtx_data + mtx_offset);
876	<	if (!sky_only)
877	<	add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
878	<	cc, dni, mtx_data + mtx_offset);
879	<	/* monthly reporting */
880	<	if (verbose && mo != last_monthly)
881	<	fprintf(stderr, "%s: stepping through month %d...\n", progname,
882	<	last_monthly = mo);
883	<	}
884	<	if (j != EOF) {
885	<	fprintf(stderr, "%s: error on input\n", progname);
886	<	exit(1);
887	<	}
888	<	EPWclose(epw); epw = NULL;
889	<	freedata(mie_dp);
890	<	if (!ntsteps) {
891	<	fprintf(stderr, "%s: no valid time steps on input\n", progname);
892	<	exit(1);
893	<	}
894	<	/* write out matrix */
895	<	if (outfmt != 'a')
896	<	SET_FILE_BINARY(stdout);
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	>	}
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);
880	>	}
881	>	/* monthly reporting */
882	>	if (verbose && mo != last_monthly) {
883	>	fprintf(stderr, "%s: stepping through month %d...\n", progname,
884	>	last_monthly = mo);
885	>	}
886	>	}
887	>	if (j != EOF) {
888	>	fprintf(stderr, "%s: error on input\n", progname);
889	>	exit(1);
890	>	}
891	>	EPWclose(epw); epw = NULL;
892	>	freedata(mie_dp);
893	>	if (!ntsteps) {
894	>	fprintf(stderr, "%s: no valid time steps on input\n", progname);
895	>	exit(1);
896	>	}
897	>	/* write out matrix */
898	>	if (outfmt != 'a') {
899	>	SET_FILE_BINARY(stdout);
900	>	}
901		#ifdef getc_unlocked
902	<	flockfile(stdout);
902	>	flockfile(stdout);
903		#endif
904	<	if (verbose)
905	<	fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
906	<	outfmt == 'a' ? "" : "binary ", nstored);
907	<	if (doheader) {
908	<	newheader("RADIANCE", stdout);
909	<	printargs(argc, argv, stdout);
910	<	printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude),
911	<	-RadToDeg(s_longitude));
912	<	printf("NROWS=%d\n", nskypatch);
913	<	printf("NCOLS=%d\n", nstored);
914	<	printf("NCOMP=%d\n", NSSAMP);
915	<	if ((outfmt == 'f') | (outfmt == 'd'))
916	<	fputendian(stdout);
917	<	fputformat((char *)getfmtname(outfmt), stdout);
918	<	putchar('\n');
919	<	}
920	<	/* patches are rows (outer sort) */
921	<	for (i = 0; i < nskypatch; i++) {
922	<	mtx_offset = NSSAMP * i;
923	<	switch (outfmt) {
924	<	case 'a':
925	<	for (j = 0; j < nstored; j++) {
926	<	int k;
927	<	for (k = 0; k < NSSAMP; k++) {
928	<	printf("%.3g ", mtx_data[mtx_offset + k]);
929	<	}
930	<	printf("\n");
931	<	mtx_offset += NSSAMP * nskypatch;
932	<	}
933	<	if (nstored > 1)
934	<	fputc('\n', stdout);
935	<	break;
936	<	case 'f':
937	<	for (j = 0; j < nstored; j++) {
938	<	putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
939	<	mtx_offset += NSSAMP * nskypatch;
940	<	}
941	<	break;
942	<	case 'd':
943	<	for (j = 0; j < nstored; j++) {
944	<	double ment[NSSAMP];
945	<	for (j = 0; j < NSSAMP; j++)
946	<	ment[j] = mtx_data[mtx_offset + j];
947	<	putbinary(ment, sizeof(double), NSSAMP, stdout);
948	<	mtx_offset += NSSAMP * nskypatch;
949	<	}
950	<	break;
951	<	}
952	<	if (ferror(stdout))
953	<	goto writerr;
954	<	}
904	>	if (verbose) {
905	>	fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
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));
913	>	printf("NROWS=%d\n", nskypatch);
914	>	printf("NCOLS=%d\n", nstored);
915	>	printf("NCOMP=%d\n", NSSAMP);
916	>	fputwlsplit(WLPART, stdout);
917	>	if ((outfmt == 'f') | (outfmt == 'd')) {
918	>	fputendian(stdout);
919	>	}
920	>	fputformat((char *)getfmtname(outfmt), stdout);
921	>	putchar('\n');
922	>	}
923	>	/* patches are rows (outer sort) */
924	>	for (i = 0; i < nskypatch; i++) {
925	>	mtx_offset = NSSAMP * i;
926	>	switch (outfmt) {
927	>	case 'a':
928	>	for (j = 0; j < nstored; j++) {
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) {
936	>	fputc('\n', stdout);
937	>	}
938	>	break;
939	>	case 'f':
940	>	for (j = 0; j < nstored; j++) {
941	>	float ment[NSSAMP];
942	>	for (k = 0; k < NSSAMP; k++) {
943	>	ment[NSSAMP-1 - 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 (k = 0; k < NSSAMP; k++) {
953	>	ment[NSSAMP-1 - k] = 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)) {
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);
971	<	exit(1);
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);
971	>	exit(1);
972		fmterr:
973	<	fprintf(stderr, "%s: weather tape format error in header\n", progname);
974	<	exit(1);
973	>	fprintf(stderr, "%s: weather tape format error in header\n", progname);
974	>	exit(1);
975		writerr:
976	<	fprintf(stderr, "%s: write error on output\n", progname);
977	<	exit(1);
976	>	fprintf(stderr, "%s: write error on output\n", progname);
977	>	exit(1);
978		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)