ViewVC Help

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

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

Comparing ray/src/gen/genssky.c (file contents):
Revision 2.8 by greg, Sat Jun 7 05:09:45 2025 UTC vs.
Revision 2.9 by greg, Fri Jul 11 18:12:25 2025 UTC

#	Line 1 | Line 1
1	–	#include "color.h"
1		#ifndef lint
2	<	static const char RCSid[] =
4	<	"$Id$";
2	>	static const char RCSid[] = "$Id$";
3		#endif
4		/* Main function for generating spectral sky */
5		/* Cloudy sky computed as weight average of clear and cie overcast sky */
6
7	+	#include "atmos.h"
8		#include "copyright.h"
9	+	#include "color.h"
10		#include "paths.h"
11	–	#include "atmos.h"
11		#include "resolu.h"
12		#include "rtio.h"
13		#include <ctype.h>
#	Line 30 | Line 29 | const double D65EFF = 203.; /* standard illuminant D65
29
30		/* Mean normalized relative daylight spectra where CCT = 6415K for overcast; */
31		const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
32	<	1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
33	<	1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
34	<	0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
32	>	1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
33	>	1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
34	>	0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
35
36		/* European and North American zones */
37	<	struct {
38	<	char zname[8]; /* time zone name (all caps) */
39	<	float zmer;    /* standard meridian */
40	<	} tzone[] = {{"YST", 135},   {"YDT", 120},   {"PST", 120},  {"PDT", 105},
41	<	{"MST", 105},   {"MDT", 90},    {"CST", 90},   {"CDT", 75},
42	<	{"EST", 75},    {"EDT", 60},    {"AST", 60},   {"ADT", 45},
43	<	{"NST", 52.5},  {"NDT", 37.5},  {"GMT", 0},    {"BST", -15},
44	<	{"CET", -15},   {"CEST", -30},  {"EET", -30},  {"EEST", -45},
45	<	{"AST", -45},   {"ADT", -60},   {"GST", -60},  {"GDT", -75},
46	<	{"IST", -82.5}, {"IDT", -97.5}, {"JST", -135}, {"NDT", -150},
47	<	{"NZST", -180}, {"NZDT", -195}, {"", 0}};
37	>	struct
38	>	{
39	>	char zname[8];                 /* time zone name (all caps) */
40	>	float zmer;                 /* standard meridian */
41	>	} tzone[] = {{"YST", 135},       {"YDT", 120},   {"PST", 120},  {"PDT", 105},
42	>	{"MST", 105},   {"MDT", 90},    {"CST", 90},   {"CDT", 75},
43	>	{"EST", 75},    {"EDT", 60},    {"AST", 60},   {"ADT", 45},
44	>	{"NST", 52.5},  {"NDT", 37.5},  {"GMT", 0},    {"BST", -15},
45	>	{"CET", -15},   {"CEST", -30},  {"EET", -30},  {"EEST", -45},
46	>	{"AST", -45},   {"ADT", -60},   {"GST", -60},  {"GDT", -75},
47	>	{"IST", -82.5}, {"IDT", -97.5}, {"JST", -135}, {"NDT", -150},
48	>	{"NZST", -180}, {"NZDT", -195}, {"", 0}};
49
50	<	static int make_directory(const char *path) {
50	>	static int
51	>	make_directory
52	>	(
53	>	const char *path
54	>	)
55	>	{
56		#ifdef _WIN32
57	<	if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
58	<	return 1;
59	<	}
60	<	return 0;
57	>	if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
58	>	return 1;
59	>	}
60	>	return 0;
61	>
62		#else
63	<	if (mkdir(path, 0777) == 0 || errno == EEXIST) {
64	<	return 1;
65	<	}
66	<	return 0;
63	>	if (mkdir(path, 0777) == 0 || errno == EEXIST) {
64	>	return 1;
65	>	}
66	>	return 0;
67	>
68		#endif
69		}
70
71	<	inline static float deg2rad(float deg) { return deg * (PI / 180.); }
71	>	inline static float
72	>	deg2rad
73	>	(
74	>	float deg
75	>	)
76	>	{
77	>	return deg * (PI / 180.);
78	>	}
79
80	<	static int cvthour(char *hs, int *tsolar, double *hour) {
81	<	char *cp = hs;
82	<	int i, j;
80	>	static int
81	>	cvthour
82	>	(
83	>	char *hs,
84	>	int *tsolar,
85	>	double *hour
86	>	)
87	>	{
88	>	char *cp = hs;
89	>	int i, j;
90
91	<	if ((*tsolar = *cp == '+'))
92	<	cp++; /* solar time? */
93	<	while (isdigit(*cp))
94	<	cp++;
95	<	if (*cp == ':')
96	<	*hour = atoi(hs) + atoi(++cp) / 60.0;
97	<	else {
98	<	*hour = atof(hs);
99	<	if (*cp == '.')
100	<	cp++;
101	<	}
102	<	while (isdigit(*cp))
103	<	cp++;
104	<	if (!*cp)
105	<	return (0);
106	<	if (*tsolar || !isalpha(*cp)) {
107	<	fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
108	<	exit(1);
109	<	}
110	<	i = 0;
111	<	do {
112	<	for (j = 0; cp[j]; j++)
113	<	if (toupper(cp[j]) != tzone[i].zname[j])
114	<	break;
115	<	if (!cp[j] && !tzone[i].zname[j]) {
116	<	s_meridian = tzone[i].zmer * (PI / 180);
117	<	return (1);
118	<	}
119	<	} while (tzone[i++].zname[0]);
91	>	if ((*tsolar = *cp == '+')) {
92	>	cp++;                                 /* solar time? */
93	>	}
94	>	while (isdigit(*cp)) {
95	>	cp++;
96	>	}
97	>	if (*cp == ':') {
98	>	*hour = atoi(hs) + atoi(++cp) / 60.0;
99	>	}else{
100	>	*hour = atof(hs);
101	>	if (*cp == '.') {
102	>	cp++;
103	>	}
104	>	}
105	>	while (isdigit(*cp)) {
106	>	cp++;
107	>	}
108	>	if (!*cp) {
109	>	return (0);
110	>	}
111	>	if (*tsolar || !isalpha(*cp)) {
112	>	fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
113	>	exit(1);
114	>	}
115	>	i = 0;
116	>	do {
117	>	for (j = 0; cp[j]; j++) {
118	>	if (toupper(cp[j]) != tzone[i].zname[j]) {
119	>	break;
120	>	}
121	>	}
122	>	if (!cp[j] && !tzone[i].zname[j]) {
123	>	s_meridian = tzone[i].zmer * (PI / 180);
124	>	return (1);
125	>	}
126	>	} while (tzone[i++].zname[0]);
127
128	<	fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
129	<	fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
130	<	for (i = 1; tzone[i].zname[0]; i++)
131	<	fprintf(stderr, " %s", tzone[i].zname);
132	<	putc('\n', stderr);
133	<	exit(1);
128	>	fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
129	>	fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
130	>	for (i = 1; tzone[i].zname[0]; i++) {
131	>	fprintf(stderr, " %s", tzone[i].zname);
132	>	}
133	>	putc('\n', stderr);
134	>	exit(1);
135		}
136
137	<	static void basename(const char *path, char *output, size_t outsize) {
138	<	const char *last_slash = strrchr(path, '/');
139	<	const char *last_backslash = strrchr(path, '\\');
140	<	const char *filename = path;
141	<	const char *last_dot;
137	>	static void
138	>	basename
139	>	(
140	>	const char *path,
141	>	char *output,
142	>	size_t outsize
143	>	)
144	>	{
145	>	const char *last_slash = strrchr(path, '/');
146	>	const char *last_backslash = strrchr(path, '\\');
147	>	const char *filename = path;
148	>	const char *last_dot;
149
150	<	if (last_slash && last_backslash) {
151	<	filename =
152	<	(last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
153	<	} else if (last_slash) {
154	<	filename = last_slash + 1;
155	<	} else if (last_backslash) {
156	<	filename = last_backslash + 1;
157	<	}
150	>	if (last_slash && last_backslash) {
151	>	filename =
152	>	(last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
153	>	} else if (last_slash) {
154	>	filename = last_slash + 1;
155	>	} else if (last_backslash) {
156	>	filename = last_backslash + 1;
157	>	}
158
159	<	last_dot = strrchr(filename, '.');
160	<	if (last_dot) {
161	<	size_t length = last_dot - filename;
162	<	if (length < outsize) {
163	<	strncpy(output, filename, length);
164	<	output[length] = '\0';
165	<	} else {
166	<	strncpy(output, filename, outsize - 1);
167	<	output[outsize - 1] = '\0';
168	<	}
169	<	}
159	>	last_dot = strrchr(filename, '.');
160	>	if (last_dot) {
161	>	size_t length = last_dot - filename;
162	>	if (length < outsize) {
163	>	strncpy(output, filename, length);
164	>	output[length] = '\0';
165	>	} else {
166	>	strncpy(output, filename, outsize - 1);
167	>	output[outsize - 1] = '\0';
168	>	}
169	>	}
170		}
171
172	<	static char *join_paths(const char *path1, const char *path2) {
173	<	size_t len1 = strlen(path1);
174	<	size_t len2 = strlen(path2);
175	<	int need_separator = (path1[len1 - 1] != DIRSEP);
172	>	static char *
173	>	join_paths
174	>	(
175	>	const char *path1,
176	>	const char *path2
177	>	)
178	>	{
179	>	size_t len1 = strlen(path1);
180	>	size_t len2 = strlen(path2);
181	>	int need_separator = (path1[len1 - 1] != DIRSEP);
182
183	<	char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
184	<	if (!result)
185	<	return NULL;
183	>	char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
184	>	if (!result) {
185	>	return NULL;
186	>	}
187
188	<	strcpy(result, path1);
189	<	if (need_separator) {
190	<	result[len1] = DIRSEP;
191	<	len1++;
192	<	}
193	<	strcpy(result + len1, path2);
188	>	strcpy(result, path1);
189	>	if (need_separator) {
190	>	result[len1] = DIRSEP;
191	>	len1++;
192	>	}
193	>	strcpy(result + len1, path2);
194
195	<	return result;
195	>	return result;
196		}
197
198	<	static inline double wmean2(const double a, const double b, const double x) {
199	<	return a * (1 - x) + b * x;
198	>	static inline double
199	>	wmean2
200	>	(
201	>	const double a,
202	>	const double b,
203	>	const double x
204	>	)
205	>	{
206	>	return a * (1 - x) + b * x;
207		}
208
209	<	static inline double wmean(const double a, const double x, const double b,
210	<	const double y) {
211	<	return (a * x + b * y) / (a + b);
209	>	static inline double
210	>	wmean
211	>	(
212	>	const double a,
213	>	const double x,
214	>	const double b,
215	>	const double y
216	>	)
217	>	{
218	>	return (a * x + b * y) / (a + b);
219		}
220
221	<	static double get_overcast_zenith_brightness(const double sundir[3]) {
222	<	double zenithbr;
223	<	if (sundir[2] < 0) {
224	<	zenithbr = 0;
225	<	} else {
226	<	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
227	<	}
228	<	return zenithbr;
221	>	static double
222	>	get_overcast_zenith_brightness
223	>	(
224	>	const double sundir[3]
225	>	)
226	>	{
227	>	double zenithbr;
228	>	if (sundir[2] < 0) {
229	>	zenithbr = 0;
230	>	} else {
231	>	zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
232	>	}
233	>	return zenithbr;
234		}
235
236		/* from gensky.c */
237	<	static double get_overcast_brightness(const double dz, const double zenithbr) {
238	<	double groundbr = zenithbr * GNORM;
239	<	return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
240	<	pow(dz + 1.01, -10), groundbr);
237	>	static double
238	>	get_overcast_brightness
239	>	(
240	>	const double dz,
241	>	const double zenithbr
242	>	)
243	>	{
244	>	double groundbr = zenithbr * GNORM;
245	>	return wmean(
246	>	pow(dz + 1.01, 10),
247	>	zenithbr * (1 + 2 * dz) / 3,
248	>	pow(dz + 1.01, -10),
249	>	groundbr);
250		}
251
252	<	static void write_header(const int argc, char **argv, const double cloud_cover,
253	<	const double grefl, const int res) {
254	<	int i;
255	<	printf("# ");
256	<	for (i = 0; i < argc; i++) {
257	<	printf("%s ", argv[i]);
258	<	}
259	<	printf("\n");
260	<	printf(
261	<	"#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
262	<	cloud_cover, grefl, res);
252	>	static void
253	>	write_header
254	>	(
255	>	const int argc,
256	>	char **argv,
257	>	const double cloud_cover,
258	>	const double grefl,
259	>	const int res
260	>	)
261	>	{
262	>	int i;
263	>	printf("# ");
264	>	for (i = 0; i < argc; i++) {
265	>	printf("%s ", argv[i]);
266	>	}
267	>	printf("\n");
268	>	printf(
269	>	"#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: "
270	>	"%d\n\n",
271	>	cloud_cover,
272	>	grefl,
273	>	res);
274		}
275
276	<	static void write_rad(const double *sun_radiance, const double intensity,
277	<	const FVECT sundir, const char *ddir,
278	<	const char *skyfile) {
279	<	if (sundir[2] > 0) {
280	<	printf("void spectrum sunrad\n0\n0\n22 380 780 ");
281	<	int i;
282	<	for (i = 0; i < NSSAMP; ++i) {
283	<	printf("%.3f ", sun_radiance[i]);
284	<	}
285	<	printf("\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n", intensity,
286	<	intensity, intensity);
287	<	printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
288	<	sundir[2]);
289	<	}
290	<	printf("void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
291	<	"'1-Acos(Dz)/PI'\n0\n0\n\n",
292	<	skyfile);
276	>	static void
277	>	write_rad
278	>	(
279	>	const double *sun_radiance,
280	>	const double intensity,
281	>	const FVECT sundir,
282	>	const char *ddir,
283	>	const char *skyfile
284	>	)
285	>	{
286	>	if (sundir[2] > 0) {
287	>	printf("void spectrum sunrad\n0\n0\n22 380 780 ");
288	>	int i;
289	>	for (i = 0; i < NSSAMP; ++i) {
290	>	printf("%.3f ", sun_radiance[i]);
291	>	}
292	>	printf(
293	>	"\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n",
294	>	intensity,
295	>	intensity,
296	>	intensity);
297	>	printf(
298	>	"solar source sun\n0\n0\n4 %f %f %f 0.533\n\n",
299	>	sundir[0],
300	>	sundir[1],
301	>	sundir[2]);
302	>	}
303	>	printf(
304	>	"void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
305	>	"'1-Acos(Dz)/PI'\n0\n0\n\n",
306	>	skyfile);
307		}
308
309	<	static void write_hsr_header(FILE *fp, RESOLU *res) {
310	<	float wvsplit[4] = {380, 480, 588, 780};
311	<	newheader("RADIANCE", fp);
312	<	fputncomp(NSSAMP, fp);
313	<	fputwlsplit(wvsplit, fp);
314	<	fputformat(SPECFMT, fp);
315	<	fputc('\n', fp);
316	<	fputsresolu(res, fp);
309	>	static void
310	>	write_hsr_header
311	>	(
312	>	FILE *fp,
313	>	RESOLU *res
314	>	)
315	>	{
316	>	newheader("RADIANCE", fp);
317	>	fputncomp(NSSAMP, fp);
318	>	fputwlsplit(WLPART, fp);
319	>	fputformat(SPECFMT, fp);
320	>	fputc('\n', fp);
321	>	fputsresolu(res, fp);
322		}
323
324	<	static inline float frac(float x) { return x - floor(x); }
324	>	static void
325	>	reverse_array_float
326	>	(
327	>	float arr[],
328	>	int size
329	>	)
330	>	{
331	>	int start = 0;
332	>	int end = size - 1;
333
334	<	int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
335	<	DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
336	<	const double cloud_cover, const FVECT sundir,
337	<	const double grefl, const int res, const char *outname,
338	<	const char *ddir, const double dirnorm, const double difhor) {
339	<	char skyfile[PATH_MAX];
340	<	if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP,
341	<	outname)) {
233	<	fprintf(stderr, "Error setting sky file name\n");
234	<	return 0;
235	<	};
236	<	int xres = res;
237	<	int yres = xres / 2;
238	<	RESOLU rs = {PIXSTANDARD, xres, yres};
239	<	FILE *skyfp = fopen(skyfile, "w");
240	<	write_hsr_header(skyfp, &rs);
334	>	while (start < end) {
335	>	float temp = arr[start];
336	>	arr[start] = arr[end];
337	>	arr[end] = temp;
338	>	start++;
339	>	end--;
340	>	}
341	>	}
342
343	<	CNDX[3] = NSSAMP;
343	>	int
344	>	gen_spect_sky
345	>	(
346	>	DATARRAY *tau_clear,
347	>	DATARRAY *scat_clear,
348	>	DATARRAY *scat1m_clear,
349	>	DATARRAY *irrad_clear,
350	>	const double cloud_cover,
351	>	const FVECT sundir,
352	>	const double grefl,
353	>	const int res,
354	>	const char *outname,
355	>	const char *ddir,
356	>	const double dirnorm,
357	>	const double difhor
358	>	)
359	>	{
360	>	char skyfile[PATH_MAX];
361	>	if (!snprintf(
362	>	skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP, outname)) {
363	>	fprintf(stderr, "Error setting sky file name\n");
364	>	return 0;
365	>	}
366	>	;
367	>	int xres = res;
368	>	int yres = xres / 2;
369	>	RESOLU rs = {PIXSTANDARD, xres, yres};
370	>	FILE *skyfp = fopen(skyfile, "w");
371	>	write_hsr_header(skyfp, &rs);
372
373	<	FVECT view_point = {0, 0, ER + 10};
245	<	const double radius = VLEN(view_point);
246	<	const double sun_ct = fdot(view_point, sundir) / radius;
373	>	CNDX[3] = NSSAMP;
374
375	<	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
376	<	double overcast_grndbr = overcast_zenithbr * GNORM;
375	>	FVECT view_point = {0, 0, ER + 10};
376	>	const double radius = VLEN(view_point);
377	>	const double sun_ct = fdot(view_point, sundir) / radius;
378
379	<	double dif_ratio = 1;
380	<	if (difhor > 0) {
253	<	DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad_clear, radius, sun_ct);
254	<	double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
255	<	double diffuse_irradiance = 0;
256	<	int l;
257	<	for (l = 0; l < NSSAMP; ++l) {
258	<	diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20;  /* 20nm interval */
259	<	}
260	<	free(indirect_irradiance_clear);
261	<	diffuse_irradiance = wmean2(diffuse_irradiance, overcast_ghi, cloud_cover);
262	<	if (diffuse_irradiance > 0) {
263	<	dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;       /* fudge */
264	<	}
265	<	}
266	<	int i, j, k;
267	<	for (j = 0; j < yres; ++j) {
268	<	for (i = 0; i < xres; ++i) {
269	<	SCOLOR radiance = {0};
270	<	SCOLR sky_sclr = {0};
379	>	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
380	>	double overcast_grndbr = overcast_zenithbr * GNORM;
381
382	<	float px = i / (xres - 1.0);
383	<	float py = j / (yres - 1.0);
384	<	float lambda = ((1 - py) * PI) - (PI / 2.0);
385	<	float phi = (px * 2.0 * PI) - PI;
382	>	double dif_ratio = 1;
383	>	if (difhor > 0) {
384	>	DATARRAY *indirect_irradiance_clear =
385	>	get_indirect_irradiance(irrad_clear, radius, sun_ct);
386	>	double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
387	>	double diffuse_irradiance = 0;
388	>	int l;
389	>	for (l = 0; l < NSSAMP; ++l) {
390	>	diffuse_irradiance +=
391	>	indirect_irradiance_clear->arr.d[l] * 20;                                                                 /* 20nm interval */
392	>	}
393	>	free(indirect_irradiance_clear);
394	>	diffuse_irradiance =
395	>	wmean2(diffuse_irradiance, overcast_ghi, cloud_cover);
396	>	if (diffuse_irradiance > 0) {
397	>	dif_ratio =
398	>	difhor / WHTEFFICACY / diffuse_irradiance / 1.15;                                                                 /* fudge */
399	>	}
400	>	}
401	>	int i, j, k;
402	>	for (j = 0; j < yres; ++j) {
403	>	for (i = 0; i < xres; ++i) {
404	>	SCOLOR radiance = {0};
405	>	SCOLR sky_sclr = {0};
406
407	<	FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
408	<	sin(lambda)};
407	>	float px = i / (xres - 1.0);
408	>	float py = j / (yres - 1.0);
409	>	float lambda = ((1 - py) * PI) - (PI / 2.0);
410	>	float phi = (px * 2.0 * PI) - PI;
411
412	<	const double mu = fdot(view_point, rdir) / radius;
413	<	const double nu = fdot(rdir, sundir);
412	>	FVECT rdir = {
413	>	cos(lambda) * cos(phi), cos(lambda) * sin(phi), sin(lambda)
414	>	};
415
416	<	/* hit ground */
417	<	if (rdir[2] < 0) {
285	<	get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
286	<	view_point, rdir, radius, mu, sun_ct, nu, grefl,
287	<	sundir, radiance);
288	<	} else {
289	<	get_sky_radiance(scat_clear, scat1m_clear, radius, mu, sun_ct, nu,
290	<	radiance);
291	<	}
416	>	const double mu = fdot(view_point, rdir) / radius;
417	>	const double nu = fdot(rdir, sundir);
418
419	<	for (k = 0; k < NSSAMP; ++k) {
420	<	radiance[k] *= WVLSPAN;
421	<	}
419	>	/* hit ground */
420	>	if (rdir[2] < 0) {
421	>	get_ground_radiance(
422	>	tau_clear,
423	>	scat_clear,
424	>	scat1m_clear,
425	>	irrad_clear,
426	>	view_point,
427	>	rdir,
428	>	radius,
429	>	mu,
430	>	sun_ct,
431	>	nu,
432	>	grefl,
433	>	sundir,
434	>	radiance);
435	>	} else {
436	>	get_sky_radiance(
437	>	scat_clear, scat1m_clear, radius, mu, sun_ct, nu, radiance);
438	>	}
439
440	<	if (cloud_cover > 0) {
441	<	double skybr = get_overcast_brightness(rdir[2], overcast_zenithbr);
442	<	if (rdir[2] < 0) {
300	<	for (k = 0; k < NSSAMP; ++k) {
301	<	radiance[k] = wmean2(radiance[k], overcast_grndbr * D6415[k], cloud_cover);
302	<	}
303	<	} else {
304	<	for (k = 0; k < NSSAMP; ++k) {
305	<	radiance[k] = wmean2(radiance[k], skybr * D6415[k], cloud_cover);
306	<	}
307	<	}
308	<	}
440	>	for (k = 0; k < NSSAMP; ++k) {
441	>	radiance[k] *= WVLSPAN;
442	>	}
443
444	<	for (k = 0; k < NSSAMP; ++k) {
445	<	radiance[k] *= dif_ratio;
446	<	}
444	>	if (cloud_cover > 0) {
445	>	double skybr =
446	>	get_overcast_brightness(rdir[2], overcast_zenithbr);
447	>	if (rdir[2] < 0) {
448	>	for (k = 0; k < NSSAMP; ++k) {
449	>	radiance[k] = wmean2(
450	>	radiance[k],
451	>	overcast_grndbr * D6415[k],
452	>	cloud_cover);
453	>	}
454	>	} else {
455	>	for (k = 0; k < NSSAMP; ++k) {
456	>	radiance[k] =
457	>	wmean2(radiance[k], skybr * D6415[k], cloud_cover);
458	>	}
459	>	}
460	>	}
461
462	<	scolor2scolr(sky_sclr, radiance, NSSAMP);
463	<	putbinary(sky_sclr, LSCOLR, 1, skyfp);
464	<	}
317	<	}
318	<	fclose(skyfp);
462	>	for (k = 0; k < NSSAMP; ++k) {
463	>	radiance[k] *= dif_ratio;
464	>	}
465
466	<	/* Get solar radiance */
321	<	double sun_radiance[NSSAMP] = {0};
322	<	get_solar_radiance(tau_clear, scat_clear, scat1m_clear, sundir, radius,
323	<	sun_ct, sun_radiance);
324	<	if (cloud_cover > 0) {
325	<	double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
326	<	int i;
327	<	for (i = 0; i < NSSAMP; ++i) {
328	<	sun_radiance[i] =
329	<	wmean2(sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
330	<	}
331	<	}
466	>	reverse_array_float(radiance, NSSAMP);
467
468	<	/* Normalize */
469	<	double sum = 0.0;
470	<	for (i = 0; i < NSSAMP; ++i) {
471	<	sum += sun_radiance[i];
472	<	}
338	<	double mean = sum / NSSAMP;
339	<	for (i = 0; i < NSSAMP; ++i) {
340	<	sun_radiance[i] /= mean;
341	<	}
342	<	double intensity = mean * WVLSPAN;
343	<	if (dirnorm > 0) {
344	<	intensity = dirnorm / SOLOMG / WHTEFFICACY;
345	<	}
468	>	scolor2scolr(sky_sclr, radiance, NSSAMP);
469	>	putbinary(sky_sclr, LSCOLR, 1, skyfp);
470	>	}
471	>	}
472	>	fclose(skyfp);
473
474	<	write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
475	<	return 1;
474	>	/* Get solar radiance */
475	>	double sun_radiance[NSSAMP] = {0};
476	>	get_solar_radiance(
477	>	tau_clear,
478	>	scat_clear,
479	>	scat1m_clear,
480	>	sundir,
481	>	radius,
482	>	sun_ct,
483	>	sun_radiance);
484	>	if (cloud_cover > 0) {
485	>	double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
486	>	int i;
487	>	for (i = 0; i < NSSAMP; ++i) {
488	>	sun_radiance[i] = wmean2(
489	>	sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
490	>	}
491	>	}
492	>
493	>	/* Normalize */
494	>	double sum = 0.0;
495	>	for (i = 0; i < NSSAMP; ++i) {
496	>	sum += sun_radiance[i];
497	>	}
498	>	double mean = sum / NSSAMP;
499	>	for (i = 0; i < NSSAMP; ++i) {
500	>	sun_radiance[i] /= mean;
501	>	}
502	>	double intensity = mean * WVLSPAN;
503	>	if (dirnorm > 0) {
504	>	intensity = dirnorm / SOLOMG / WHTEFFICACY;
505	>	}
506	>
507	>	write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
508	>	return 1;
509		}
510
511	<	static DpPaths get_dppaths(const char *dir, const double aod, const char *mname,
512	<	const char *tag) {
513	<	DpPaths paths;
511	>	static DpPaths
512	>	get_dppaths
513	>	(
514	>	const char *dir,
515	>	const double aod,
516	>	const char *mname,
517	>	const char *tag
518	>	)
519	>	{
520	>	DpPaths paths;
521
522	<	snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag,
523	<	mname, aod);
524	<	snprintf(paths.scat, PATH_MAX, "%s%cscat_%s_%s_%.2f.dat", dir, DIRSEP, tag,
525	<	mname, aod);
526	<	snprintf(paths.scat1m, PATH_MAX, "%s%cscat1m_%s_%s_%.2f.dat", dir, DIRSEP,
527	<	tag, mname, aod);
528	<	snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag,
529	<	mname, aod);
522	>	snprintf(
523	>	paths.tau,
524	>	PATH_MAX,
525	>	"%s%ctau_%s_%s_%.2f.dat",
526	>	dir,
527	>	DIRSEP,
528	>	tag,
529	>	mname,
530	>	aod);
531	>	snprintf(
532	>	paths.scat,
533	>	PATH_MAX,
534	>	"%s%cscat_%s_%s_%.2f.dat",
535	>	dir,
536	>	DIRSEP,
537	>	tag,
538	>	mname,
539	>	aod);
540	>	snprintf(
541	>	paths.scat1m,
542	>	PATH_MAX,
543	>	"%s%cscat1m_%s_%s_%.2f.dat",
544	>	dir,
545	>	DIRSEP,
546	>	tag,
547	>	mname,
548	>	aod);
549	>	snprintf(
550	>	paths.irrad,
551	>	PATH_MAX,
552	>	"%s%cirrad_%s_%s_%.2f.dat",
553	>	dir,
554	>	DIRSEP,
555	>	tag,
556	>	mname,
557	>	aod);
558
559	<	return paths;
559	>	return paths;
560		}
561
562	<	static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
563	<	const int is_summer,
564	<	const double s_latitude) {
565	<	if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
566	<	tag[0] = 's';
567	<	if (is_summer) {
568	<	tag[1] = 's';
569	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
570	<	atmos->beta_r0 = BR0_SS;
571	<	} else {
572	<	tag[1] = 'w';
573	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
574	<	atmos->beta_r0 = BR0_SW;
575	<	}
576	<	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
577	<	tag[0] = 'm';
578	<	if (is_summer) {
579	<	tag[1] = 's';
580	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
581	<	atmos->beta_r0 = BR0_MS;
582	<	} else {
583	<	tag[1] = 'w';
584	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
585	<	atmos->beta_r0 = BR0_MW;
586	<	}
587	<	} else {
588	<	tag[0] = 't';
589	<	tag[1] = 'r';
590	<	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
591	<	atmos->beta_r0 = BR0_T;
592	<	}
593	<	tag[2] = '\0';
562	>	static void
563	>	set_rayleigh_density_profile
564	>	(
565	>	Atmosphere *atmos,
566	>	char *tag,
567	>	const int is_summer,
568	>	const double s_latitude
569	>	)
570	>	{
571	>	if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
572	>	tag[0] = 's';
573	>	if (is_summer) {
574	>	tag[1] = 's';
575	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
576	>	atmos->beta_r0 = BR0_SS;
577	>	} else {
578	>	tag[1] = 'w';
579	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
580	>	atmos->beta_r0 = BR0_SW;
581	>	}
582	>	} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
583	>	tag[0] = 'm';
584	>	if (is_summer) {
585	>	tag[1] = 's';
586	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
587	>	atmos->beta_r0 = BR0_MS;
588	>	} else {
589	>	tag[1] = 'w';
590	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
591	>	atmos->beta_r0 = BR0_MW;
592	>	}
593	>	} else {
594	>	tag[0] = 't';
595	>	tag[1] = 'r';
596	>	atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
597	>	atmos->beta_r0 = BR0_T;
598	>	}
599	>	tag[2] = '\0';
600		}
601
602	<	static Atmosphere init_atmos(const double aod, const double grefl) {
603	<	Atmosphere atmos = {.ozone_density = {.layers =
604	<	{
605	<	{.width = 25000.0,
606	<	.exp_term = 0.0,
607	<	.exp_scale = 0.0,
608	<	.linear_term = 1.0 / 15000.0,
609	<	.constant_term = -2.0 / 3.0},
610	<	{.width = AH,
611	<	.exp_term = 0.0,
612	<	.exp_scale = 0.0,
613	<	.linear_term = -1.0 / 15000.0,
614	<	.constant_term = 8.0 / 3.0},
615	<	}},
616	<	.rayleigh_density = {.layers =
617	<	{
618	<	{.width = AH,
619	<	.exp_term = 1.0,
620	<	.exp_scale = -1.0 / HR_MS,
621	<	.linear_term = 0.0,
622	<	.constant_term = 0.0},
623	<	}},
624	<	.beta_r0 = BR0_MS,
625	<	.beta_scale = aod / AOD0_CA,
626	<	.beta_m = NULL,
627	<	.grefl = grefl};
628	<	return atmos;
602	>	static Atmosphere
603	>	init_atmos
604	>	(
605	>	const double aod,
606	>	const double grefl
607	>	)
608	>	{
609	>	Atmosphere atmos = {
610	>	.ozone_density =
611	>	{.layers =
612	>	{
613	>	{.width = 25000.0,
614	>	.exp_term = 0.0,
615	>	.exp_scale = 0.0,
616	>	.linear_term = 1.0 / 15000.0,
617	>	.constant_term = -2.0 / 3.0},
618	>	{.width = AH,
619	>	.exp_term = 0.0,
620	>	.exp_scale = 0.0,
621	>	.linear_term = -1.0 / 15000.0,
622	>	.constant_term = 8.0 / 3.0},
623	>	}},
624	>	.rayleigh_density =
625	>	{.layers =
626	>	{
627	>	{.width = AH,
628	>	.exp_term = 1.0,
629	>	.exp_scale = -1.0 / HR_MS,
630	>	.linear_term = 0.0,
631	>	.constant_term = 0.0},
632	>	}},
633	>	.beta_r0 = BR0_MS,
634	>	.beta_scale = aod / AOD0_CA,
635	>	.beta_m = NULL,
636	>	.grefl = grefl
637	>	};
638	>	return atmos;
639		}
640
641	<	int main(int argc, char *argv[]) {
642	<	int month, day;
643	<	double hour;
644	<	FVECT sundir;
645	<	int num_threads = 1;
646	<	int sorder = 4;
647	<	int year = 0;
648	<	int tsolar = 0;
649	<	int got_meridian = 0;
650	<	double grefl = 0.2;
651	<	double ccover = 0.0;
652	<	int res = 64;
653	<	double aod = AOD0_CA;
654	<	char *outname = "out";
655	<	char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
656	<	char mie_name[20] = "mie_ca";
657	<	char lstag[3];
658	<	char *ddir = ".";
659	<	int i;
660	<	double dirnorm = 0; /* direct normal illuminance */
661	<	double difhor = 0;  /* diffuse horizontal illuminance */
641	>	int
642	>	main
643	>	(
644	>	int argc,
645	>	char *argv[]
646	>	)
647	>	{
648	>	int month, day;
649	>	double hour;
650	>	FVECT sundir;
651	>	int num_threads = 1;
652	>	int sorder = 4;
653	>	int year = 0;
654	>	int tsolar = 0;
655	>	int got_meridian = 0;
656	>	double grefl = 0.2;
657	>	double ccover = 0.0;
658	>	int res = 64;
659	>	double aod = AOD0_CA;
660	>	char *outname = "out";
661	>	char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
662	>	char mie_name[20] = "mie_ca";
663	>	char lstag[3];
664	>	char *ddir = ".";
665	>	int i;
666	>	double dirnorm = 0;                 /* direct normal illuminance */
667	>	double difhor = 0;                 /* diffuse horizontal illuminance */
668
669	<	fixargv0(argv[0]);
670	<	if (argc == 2 && !strcmp(argv[1], "-defaults")) {
671	<	printf("-i %d\t\t\t\t#scattering order\n", sorder);
672	<	printf("-g %f\t\t\t#ground reflectance\n", grefl);
673	<	printf("-c %f\t\t\t#cloud cover\n", ccover);
674	<	printf("-r %d\t\t\t\t#image resolution\n", res);
675	<	printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
676	<	printf("-f %s\t\t\t\t#output name (-f)\n", outname);
677	<	printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
678	<	exit(0);
679	<	}
669	>	fixargv0(argv[0]);
670	>	if (argc == 2 && !strcmp(argv[1], "-defaults")) {
671	>	printf("-i %d\t\t\t\t#scattering order\n", sorder);
672	>	printf("-g %f\t\t\t#ground reflectance\n", grefl);
673	>	printf("-c %f\t\t\t#cloud cover\n", ccover);
674	>	printf("-r %d\t\t\t\t#image resolution\n", res);
675	>	printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
676	>	printf("-f %s\t\t\t\t#output name (-f)\n", outname);
677	>	printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
678	>	exit(0);
679	>	}
680
681	<	if (argc < 4) {
682	<	fprintf(stderr,
683	<	"Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r "
684	<	"res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
685	<	"-g grefl -f outpath\n",
686	<	argv[0]);
687	<	return 0;
688	<	}
681	>	if (argc < 4) {
682	>	fprintf(
683	>	stderr,
684	>	"Usage: %s month day hour -y year -a lat -o lon -m tz -d aod "
685	>	"-r "
686	>	"res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
687	>	"-g grefl -f outpath\n",
688	>	argv[0]);
689	>	return 0;
690	>	}
691
692	<	month = atoi(argv[1]);
693	<	if (month < 1 || month > 12) {
694	<	fprintf(stderr, "bad month");
695	<	exit(1);
696	<	}
697	<	day = atoi(argv[2]);
698	<	if (day < 1 || day > 31) {
699	<	fprintf(stderr, "bad month");
700	<	exit(1);
701	<	}
702	<	got_meridian = cvthour(argv[3], &tsolar, &hour);
692	>	month = atoi(argv[1]);
693	>	if (month < 1 || month > 12) {
694	>	fprintf(stderr, "bad month");
695	>	exit(1);
696	>	}
697	>	day = atoi(argv[2]);
698	>	if (day < 1 || day > 31) {
699	>	fprintf(stderr, "bad month");
700	>	exit(1);
701	>	}
702	>	got_meridian = cvthour(argv[3], &tsolar, &hour);
703
704	<	if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
705	<	fprintf(stderr, "Cannot compute solar angle\n");
706	<	exit(1);
707	<	}
704	>	if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
705	>	fprintf(stderr, "Cannot compute solar angle\n");
706	>	exit(1);
707	>	}
708
709	<	for (i = 4; i < argc; i++) {
710	<	if (argv[i][0] == '-') {
711	<	switch (argv[i][1]) {
712	<	case 'a':
713	<	s_latitude = atof(argv[++i]) * (PI / 180.0);
714	<	break;
715	<	case 'c':
716	<	ccover = atof(argv[++i]);
717	<	break;
718	<	case 'd':
719	<	aod = atof(argv[++i]);
720	<	break;
721	<	case 'f':
722	<	outname = argv[++i];
723	<	break;
724	<	case 'g':
725	<	grefl = atof(argv[++i]);
726	<	break;
727	<	case 'i':
728	<	sorder = atoi(argv[++i]);
729	<	break;
730	<	case 'l':
731	<	mie_path = argv[++i];
732	<	basename(mie_path, mie_name, sizeof(mie_name));
733	<	break;
734	<	case 'm':
735	<	if (got_meridian) {
736	<	++i;
737	<	break;
738	<	}
739	<	s_meridian = atof(argv[++i]) * (PI / 180.0);
740	<	break;
741	<	case 'n':
742	<	num_threads = atoi(argv[++i]);
743	<	break;
744	<	case 'o':
745	<	s_longitude = atof(argv[++i]) * (PI / 180.0);
746	<	break;
747	<	case 'L':
748	<	dirnorm = atof(argv[++i]);
749	<	difhor = atof(argv[++i]);
750	<	break;
751	<	case 'p':
752	<	ddir = argv[++i];
753	<	break;
754	<	case 'r':
755	<	res = atoi(argv[++i]);
756	<	break;
757	<	case 'y':
758	<	year = atoi(argv[++i]);
759	<	break;
760	<	default:
761	<	fprintf(stderr, "Unknown option %s\n", argv[i]);
762	<	exit(1);
763	<	}
764	<	}
765	<	}
766	<	if (year && (year < 1950) | (year > 2050))
767	<	fprintf(stderr, "%s: warning - year should be in range 1950-2050\n",
768	<	progname);
769	<	if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180)
770	<	fprintf(stderr,
771	<	"%s: warning - %.1f hours btwn. standard meridian and longitude\n",
772	<	progname, (s_longitude - s_meridian) * 12 / PI);
709	>	for (i = 4; i < argc; i++) {
710	>	if (argv[i][0] == '-') {
711	>	switch (argv[i][1]) {
712	>	case 'a':
713	>	s_latitude = atof(argv[++i]) * (PI / 180.0);
714	>	break;
715	>	case 'c':
716	>	ccover = atof(argv[++i]);
717	>	break;
718	>	case 'd':
719	>	aod = atof(argv[++i]);
720	>	break;
721	>	case 'f':
722	>	outname = argv[++i];
723	>	break;
724	>	case 'g':
725	>	grefl = atof(argv[++i]);
726	>	break;
727	>	case 'i':
728	>	sorder = atoi(argv[++i]);
729	>	break;
730	>	case 'l':
731	>	mie_path = argv[++i];
732	>	basename(mie_path, mie_name, sizeof(mie_name));
733	>	break;
734	>	case 'm':
735	>	if (got_meridian) {
736	>	++i;
737	>	break;
738	>	}
739	>	s_meridian = atof(argv[++i]) * (PI / 180.0);
740	>	break;
741	>	case 'n':
742	>	num_threads = atoi(argv[++i]);
743	>	break;
744	>	case 'o':
745	>	s_longitude = atof(argv[++i]) * (PI / 180.0);
746	>	break;
747	>	case 'L':
748	>	dirnorm = atof(argv[++i]);
749	>	difhor = atof(argv[++i]);
750	>	break;
751	>	case 'p':
752	>	ddir = argv[++i];
753	>	break;
754	>	case 'r':
755	>	res = atoi(argv[++i]);
756	>	break;
757	>	case 'y':
758	>	year = atoi(argv[++i]);
759	>	break;
760	>	default:
761	>	fprintf(stderr, "Unknown option %s\n", argv[i]);
762	>	exit(1);
763	>	}
764	>	}
765	>	}
766	>	if (year && (year < 1950) | (year > 2050)) {
767	>	fprintf(
768	>	stderr,
769	>	"%s: warning - year should be in range 1950-2050\n",
770	>	progname);
771	>	}
772	>	if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180) {
773	>	fprintf(
774	>	stderr,
775	>	"%s: warning - %.1f hours btwn. standard meridian and "
776	>	"longitude\n",
777	>	progname,
778	>	(s_longitude - s_meridian) * 12 / PI);
779	>	}
780
781	<	Atmosphere clear_atmos = init_atmos(aod, grefl);
781	>	Atmosphere clear_atmos = init_atmos(aod, grefl);
782
783	<	int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
784	<	if (s_latitude < 0) {
785	<	is_summer = !is_summer;
786	<	}
787	<	set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
783	>	int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
784	>	if (s_latitude < 0) {
785	>	is_summer = !is_summer;
786	>	}
787	>	set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
788
789	<	/* Load mie density data */
790	<	DATARRAY *mie_dp = getdata(mie_path);
791	<	if (mie_dp == NULL) {
792	<	fprintf(stderr, "Error reading mie data\n");
793	<	return 0;
794	<	}
795	<	clear_atmos.beta_m = mie_dp;
789	>	/* Load mie density data */
790	>	DATARRAY *mie_dp = getdata(mie_path);
791	>	if (mie_dp == NULL) {
792	>	fprintf(stderr, "Error reading mie data\n");
793	>	return 0;
794	>	}
795	>	clear_atmos.beta_m = mie_dp;
796
797	<	char gsdir[PATH_MAX];
798	<	size_t siz = strlen(ddir);
799	<	if (ISDIRSEP(ddir[siz - 1]))
800	<	ddir[siz - 1] = '\0';
801	<	snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
802	<	if (!make_directory(gsdir)) {
803	<	fprintf(stderr, "Failed creating atmos_data directory");
804	<	exit(1);
805	<	}
806	<	DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
797	>	char gsdir[PATH_MAX];
798	>	size_t siz = strlen(ddir);
799	>	if (ISDIRSEP(ddir[siz - 1])) {
800	>	ddir[siz - 1] = '\0';
801	>	}
802	>	snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
803	>	if (!make_directory(gsdir)) {
804	>	fprintf(stderr, "Failed creating atmos_data directory");
805	>	exit(1);
806	>	}
807	>	DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
808
809	<	if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
810	<	getpath(clear_paths.scat, ".", R_OK) == NULL ||
811	<	getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
812	<	getpath(clear_paths.irrad, ".", R_OK) == NULL) {
813	<	printf("# Pre-computing...\n");
814	<	if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
815	<	fprintf(stderr, "Pre-compute failed\n");
816	<	return 0;
817	<	}
818	<	}
809	>	if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
810	>	getpath(clear_paths.scat, ".", R_OK) == NULL ||
811	>	getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
812	>	getpath(clear_paths.irrad, ".", R_OK) == NULL) {
813	>	printf("# Pre-computing...\n");
814	>	if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
815	>	fprintf(stderr, "Pre-compute failed\n");
816	>	return 0;
817	>	}
818	>	}
819
820	<	DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
821	<	DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
822	<	DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
823	<	DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
820	>	DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
821	>	DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
822	>	DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
823	>	DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
824
825	<	write_header(argc, argv, ccover, grefl, res);
825	>	write_header(argc, argv, ccover, grefl, res);
826
827	<	if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
828	<	irrad_clear_dp, ccover, sundir, grefl, res, outname, ddir,
829	<	dirnorm, difhor)) {
830	<	fprintf(stderr, "gen_spect_sky failed\n");
831	<	exit(1);
832	<	}
827	>	if (!gen_spect_sky(
828	>	tau_clear_dp,
829	>	scat_clear_dp,
830	>	scat1m_clear_dp,
831	>	irrad_clear_dp,
832	>	ccover,
833	>	sundir,
834	>	grefl,
835	>	res,
836	>	outname,
837	>	ddir,
838	>	dirnorm,
839	>	difhor)) {
840	>	fprintf(stderr, "gen_spect_sky failed\n");
841	>	exit(1);
842	>	}
843
844	<	freedata(mie_dp);
845	<	freedata(tau_clear_dp);
846	<	freedata(scat_clear_dp);
847	<	freedata(irrad_clear_dp);
848	<	freedata(scat1m_clear_dp);
844	>	freedata(mie_dp);
845	>	freedata(tau_clear_dp);
846	>	freedata(scat_clear_dp);
847	>	freedata(irrad_clear_dp);
848	>	freedata(scat1m_clear_dp);
849
850	<	return 1;
850	>	return 1;
851		}

Diff Legend

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