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.2 by greg, Fri Aug 2 19:01:13 2024 UTC vs.
Revision 1.9 by greg, Wed Jul 2 01:56:35 2025 UTC

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

Diff Legend

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