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.5 by greg, Thu Apr 10 23:30:58 2025 UTC vs.
Revision 1.6 by greg, Tue Apr 15 20:15:50 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
32	–	#define rh_vector(v, i) vector(v, rh_palt[i], rh_pazi[i])
33	–
34	–	#define rh_cos(i) tsin(rh_palt[i])
35	–
36	–	#define solar_minute(jd, hr) ((24 * 60) * ((jd) - 1) + (int)((hr) * 60. + .5))
37	–
38	–
39	–	char *progname;
40	–
41	–	double altitude;   /* Solar altitude (radians) */
42	–	double azimuth;    /* Solar azimuth (radians) */
43	–	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 */
47	–	FVECT sundir;
48	–
49	–	const double ARCTIC_LAT = 67.;
50	–	const double TROPIC_LAT = 23.;
51	–	const int SUMMER_START = 4;
52	–	const int SUMMER_END = 9;
53	–	const double GNORM = 0.777778;
54	–
55	–	const double D65EFF = 203.; /* standard illuminant D65 */
56	–
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 */
63	<	#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
74	–	int nsuns = NSUNPATCH;    /* number of sun patches to use */
75	–	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	<
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);
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
127	–	static inline double wmean2(const double a, const double b, const double x) {
128	–	return a * (1 - x) + b * x;
129	–	}
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
137
137	–	static double get_overcast_zenith_brightness(const double sundir[3]) {
138	–	double zenithbr;
139	–	if (sundir[2] < 0) {
140	–	zenithbr = 0;
141	–	} else {
142	–	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
143	–	}
144	–	return zenithbr;
145	–	}
146	–
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	<	double
155	<	solar_sunset(int month, int day)
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 + (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15));
151	>	return 12 + (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
152		}
153
154
155	<	double
164	<	solar_sunrise(int month, int day)
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 - (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15));
159	>	return 12 - (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
160		}
161
162	<	int rh_init(void) {
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
277	<	static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
278	<	const int is_summer,
279	<	const double s_latitude) {
280	<	/* Set rayleigh density profile */
281	<	if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
282	<	tag[0] = 's';
283	<	if (is_summer) {
284	<	tag[1] = 's';
285	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
286	<	atmos->beta_r0 = BR0_SS;
287	<	} else {
288	<	tag[1] = 'w';
289	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
290	<	atmos->beta_r0 = BR0_SW;
291	<	}
292	<	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
293	<	tag[0] = 'm';
294	<	if (is_summer) {
295	<	tag[1] = 's';
296	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
297	<	atmos->beta_r0 = BR0_MS;
298	<	} else {
299	<	tag[1] = 'w';
300	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
301	<	atmos->beta_r0 = BR0_MW;
302	<	}
303	<	} else {
304	<	tag[0] = 't';
305	<	tag[1] = 'r';
306	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
307	<	atmos->beta_r0 = BR0_T;
308	<	}
309	<	tag[2] = '\0';
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, double dirnorm, 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_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;
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	<	}
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	<	}
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
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};
387
388	<	get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
389	<	for (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 skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr);
403	<	int k;
404	<	for (k = 0; k < NSSAMP; ++k) {
405	<	sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
406	<	}
403	<	}
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	<	/* calibration */
409	<	for (k = 0; k < NSSAMP; ++k) {
410	<	sky_radiance[k] *= dif_ratio;
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	<	for (k = 0; k < NSSAMP; ++k) {
416	<	parr[NSSAMP * i + k] = sky_radiance[k];
417	<	}
418	<	}
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
416	–	/* Return maximum of two doubles */
417	–	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, double difhor, float *parr) {
430	<	float sun_zenith;
431	<	SCOLOR sky_radiance = {0};
432	<	SCOLOR ground_radiance = {0};
433	<	SCOLR sky_sclr = {0};
434	<	SCOLR ground_sclr = {0};
435	<	FVECT view_point = {0, 0, ER};
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);
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	<	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
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	<	}
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, parr);
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	<	EPWheader   *epw = NULL;    /* EPW/WEA input file */
498	<	EPWrecord   erec;           /* current EPW/WEA input record */
499	<	int         doheader = 1; /* output header? */
500	<	double      rotation = 0.0;
501	<	double      elevation = 0;
502	<	int         leap_day = 0;       /* add leap day? */
503	<	int         sun_hours_only = 0; /* only output sun hours? */
504	<	int         dir_is_horiz;           /* direct is meas. on horizontal? */
505	<	float       *mtx_data = NULL;
506	<	int         ntsteps = 0;      /* number of time steps */
507	<	int         tstorage = 0;     /* number of allocated time steps */
508	<	int         nstored = 0;      /* number of time steps in matrix */
509	<	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;
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(
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;
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	<	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	<	}
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	<	rh_init();
615	>	rh_init();
616
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	<	}
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	<	s_latitude = DegToRad(s_latitude);
628	<	s_longitude = DegToRad(s_longitude);
629	<	s_meridian = DegToRad(s_meridian);
630	<	/* initial allocation */
631	<	mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);
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	<	/* 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	<	}
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	<	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	<	}
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	<	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;
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	<	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;
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 > -DegToRad(SUN_ANG_DEG / 2.));
675	>	altitude = salt(sda, st);
676	>	sun_in_sky = (altitude > -deg_to_rad(SUN_ANG_DEG / 2.));
677
678	<	azimuth = sazi(sda, st) + PI - DegToRad(rotation);
678	>	azimuth = sazi(sda, st) + PI - deg_to_rad(rotation);
679
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;
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	<	dni = erec.dirillum;
687	<	dhi = erec.diffillum;
678	<	printf("%d %d %f %f %f %f %f\n", mo, da, hr, dni, dhi, aod, cc);
686	>	dni = erec.dirillum;
687	>	dhi = erec.diffillum;
688
689	<	mtx_offset = NSSAMP * nskypatch * nstored;
690	<	nstored += 1;
689	>	mtx_offset = NSSAMP * nskypatch * nstored;
690	>	nstored += 1;
691
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 */
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);
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;
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);
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	<	}
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);
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, 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);
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	<	int k;
786	<	for (k = 0; k < NSSAMP; k++) {
787	<	printf("%.3g ", mtx_data[mtx_offset + k]);
788	<	}
789	<	printf("\n");
790	<	mtx_offset += NSSAMP * nskypatch;
791	<	}
792	<	if (nstored > 1)
793	<	fputc('\n', stdout);
794	<	break;
795	<	case 'f':
796	<	for (j = 0; j < nstored; j++) {
797	<	putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
798	<	mtx_offset += NSSAMP * nskypatch;
799	<	}
800	<	break;
801	<	case 'd':
802	<	for (j = 0; j < nstored; j++) {
803	<	double ment[NSSAMP];
804	<	for (j = 0; j < NSSAMP; j++)
805	<	ment[j] = mtx_data[mtx_offset + j];
806	<	putbinary(ment, sizeof(double), NSSAMP, stdout);
807	<	mtx_offset += NSSAMP * nskypatch;
808	<	}
809	<	break;
810	<	}
811	<	if (ferror(stdout))
812	<	goto writerr;
813	<	}
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