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.2 by greg, Fri Jul 19 23:38:28 2024 UTC vs.
Revision 2.8 by greg, Sat Jun 7 05:09:45 2025 UTC

#	Line 1 | Line 1
1	+	#include "color.h"
2		#ifndef lint
3	<	static const char RCSid[] = "$Id$";
3	>	static const char RCSid[] =
4	>	"$Id$";
5		#endif
6		/* Main function for generating spectral sky */
7		/* Cloudy sky computed as weight average of clear and cie overcast sky */
8
7	–	#include "atmos.h"
9		#include "copyright.h"
10	+	#include "paths.h"
11	+	#include "atmos.h"
12		#include "resolu.h"
13		#include "rtio.h"
11	–	#include "view.h"
14		#include <ctype.h>
15		#ifdef _WIN32
16		#include <windows.h>
#	Line 18 | Line 20 | static const char RCSid[] = "$Id$";
20		#include <sys/types.h>
21		#endif
22
21	–	char *progname;
22	–
23		const double ARCTIC_LAT = 67.;
24		const double TROPIC_LAT = 23.;
25		const int SUMMER_START = 4;
#	Line 61 | Line 61 | static int make_directory(const char *path) {
61		#endif
62		}
63
64	+	inline static float deg2rad(float deg) { return deg * (PI / 180.); }
65	+
66		static int cvthour(char *hs, int *tsolar, double *hour) {
67		char *cp = hs;
68		int i, j;
#	Line 131 | Line 133 | static void basename(const char *path, char *output, s
133		}
134		}
135
136	<	char *join_paths(const char *path1, const char *path2) {
136	>	static char *join_paths(const char *path1, const char *path2) {
137		size_t len1 = strlen(path1);
138		size_t len2 = strlen(path2);
139		int need_separator = (path1[len1 - 1] != DIRSEP);
#	Line 159 | Line 161 | static inline double wmean(const double a, const doubl
161		return (a * x + b * y) / (a + b);
162		}
163
164	<	static double get_zenith_brightness(const double sundir[3]) {
164	>	static double get_overcast_zenith_brightness(const double sundir[3]) {
165		double zenithbr;
166		if (sundir[2] < 0) {
167		zenithbr = 0;
#	Line 178 | Line 180 | static double get_overcast_brightness(const double dz,
180
181		static void write_header(const int argc, char **argv, const double cloud_cover,
182		const double grefl, const int res) {
183	+	int i;
184		printf("# ");
185	<	for (int i = 0; i < argc; i++) {
185	>	for (i = 0; i < argc; i++) {
186		printf("%s ", argv[i]);
187		}
188		printf("\n");
189	<	printf("#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
190	<	cloud_cover, grefl, res);
189	>	printf(
190	>	"#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
191	>	cloud_cover, grefl, res);
192		}
193
194	<	static void write_rad(const double *sun_radiance, const FVECT sundir,
195	<	const char skyfile[PATH_MAX],
196	<	const char grndfile[PATH_MAX]) {
194	>	static void write_rad(const double *sun_radiance, const double intensity,
195	>	const FVECT sundir, const char *ddir,
196	>	const char *skyfile) {
197		if (sundir[2] > 0) {
198		printf("void spectrum sunrad\n0\n0\n22 380 780 ");
199	<	/* Normalize to one */
200	<	double sum = 0.0;
201	<	for (int i = 0; i < NSSAMP; ++i) {
198	<	sum += sun_radiance[i];
199	>	int i;
200	>	for (i = 0; i < NSSAMP; ++i) {
201	>	printf("%.3f ", sun_radiance[i]);
202		}
200	–	double mean = sum / NSSAMP;
201	–	for (int i = 0; i < NSSAMP; ++i) {
202	–	printf("%.3f ", sun_radiance[i] / mean);
203	–	}
204	–	double intensity = mean * WVLSPAN;
203		printf("\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n", intensity,
204		intensity, intensity);
205		printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
206		sundir[2]);
207		}
208	<	printf("void specpict skymap\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
209	<	"-mx\n0\n0\n\n",
208	>	printf("void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
209	>	"'1-Acos(Dz)/PI'\n0\n0\n\n",
210		skyfile);
213	–
214	–	printf("void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
215	–	"-my\n0\n0\n\n",
216	–	grndfile);
217	–	printf("void mixfunc skyfunc\n4 skymap grndmap if(Dz,1,0) .\n0\n0\n");
211		}
212
213		static void write_hsr_header(FILE *fp, RESOLU *res) {
214	<	float wvsplit[4] = {380, 480, 588,
222	<	780}; /* RGB wavelength limits+partitions (nm) */
214	>	float wvsplit[4] = {380, 480, 588, 780};
215		newheader("RADIANCE", fp);
216		fputncomp(NSSAMP, fp);
217		fputwlsplit(wvsplit, fp);
#	Line 228 | Line 220 | static void write_hsr_header(FILE *fp, RESOLU *res) {
220		fputsresolu(res, fp);
221		}
222
223	+	static inline float frac(float x) { return x - floor(x); }
224	+
225		int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
226		DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
227		const double cloud_cover, const FVECT sundir,
228	<	const double grefl, const int res, const char *outname) {
228	>	const double grefl, const int res, const char *outname,
229	>	const char *ddir, const double dirnorm, const double difhor) {
230		char skyfile[PATH_MAX];
231	<	char grndfile[PATH_MAX];
232	<	if (!snprintf(skyfile, sizeof(skyfile), "%s_sky.hsr", outname)) {
231	>	if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP,
232	>	outname)) {
233		fprintf(stderr, "Error setting sky file name\n");
234		return 0;
235		};
236	<	if (!snprintf(grndfile, sizeof(grndfile), "%s_ground.hsr", outname)) {
237	<	fprintf(stderr, "Error setting ground file name\n");
238	<	return 0;
244	<	}
245	<	RESOLU rs = {PIXSTANDARD, res, res};
236	>	int xres = res;
237	>	int yres = xres / 2;
238	>	RESOLU rs = {PIXSTANDARD, xres, yres};
239		FILE *skyfp = fopen(skyfile, "w");
247	–	FILE *grndfp = fopen(grndfile, "w");
248	–	write_hsr_header(grndfp, &rs);
240		write_hsr_header(skyfp, &rs);
250	–	VIEW skyview = {VT_ANG, {0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, 1.,
251	–	180.,   180.,         0.,           0.,           0.,
252	–	0.,     {0., 0., 0.}, {0., 0., 0.}, 0.,           0.};
253	–	VIEW grndview = {
254	–	VT_ANG, {0., 0., 0.}, {0., 0., -1.}, {0., 1., 0.}, 1., 180., 180., 0., 0.,
255	–	0.,     0.,           {0., 0., 0.},  {0., 0., 0.}, 0., 0.};
256	–	setview(&skyview);
257	–	setview(&grndview);
241
242		CNDX[3] = NSSAMP;
243
244	<	FVECT view_point = {0, 0, ER};
244	>	FVECT view_point = {0, 0, ER + 10};
245		const double radius = VLEN(view_point);
246		const double sun_ct = fdot(view_point, sundir) / radius;
247	<	for (unsigned int j = 0; j < res; ++j) {
248	<	for (unsigned int i = 0; i < res; ++i) {
249	<	RREAL loc[2];
250	<	FVECT rorg = {0};
251	<	FVECT rdir_sky = {0};
252	<	FVECT rdir_grnd = {0};
253	<	SCOLOR sky_radiance = {0};
254	<	SCOLOR ground_radiance = {0};
247	>
248	>	double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
249	>	double overcast_grndbr = overcast_zenithbr * GNORM;
250	>
251	>	double dif_ratio = 1;
252	>	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};
273	–	SCOLR ground_sclr = {0};
271
272	<	pix2loc(loc, &rs, i, j);
273	<	viewray(rorg, rdir_sky, &skyview, loc[0], loc[1]);
274	<	viewray(rorg, rdir_grnd, &grndview, loc[0], loc[1]);
272	>	float px = i / (xres - 1.0);
273	>	float py = j / (yres - 1.0);
274	>	float lambda = ((1 - py) * PI) - (PI / 2.0);
275	>	float phi = (px * 2.0 * PI) - PI;
276
277	<	const double mu_sky = fdot(view_point, rdir_sky) / radius;
278	<	const double nu_sky = fdot(rdir_sky, sundir);
277	>	FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
278	>	sin(lambda)};
279
280	<	const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
281	<	const double nu_grnd = fdot(rdir_grnd, sundir);
280	>	const double mu = fdot(view_point, rdir) / radius;
281	>	const double nu = fdot(rdir, sundir);
282
283	<	get_sky_radiance(scat_clear, scat1m_clear, radius, mu_sky, sun_ct, nu_sky,
284	<	sky_radiance);
285	<	get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
286	<	view_point, rdir_grnd, radius, mu_grnd, sun_ct,
287	<	nu_grnd, grefl, sundir, ground_radiance);
283	>	/* hit ground */
284	>	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	>	}
292
293	<	for (int k = 0; k < NSSAMP; ++k) {
294	<	sky_radiance[k] *= WVLSPAN;
293	<	ground_radiance[k] *= WVLSPAN;
293	>	for (k = 0; k < NSSAMP; ++k) {
294	>	radiance[k] *= WVLSPAN;
295		}
296
297		if (cloud_cover > 0) {
298	<	double zenithbr = get_zenith_brightness(sundir);
299	<	double grndbr = zenithbr * GNORM;
300	<	double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
301	<	for (int k = 0; k < NSSAMP; ++k) {
302	<	sky_radiance[k] =
303	<	wmean2(sky_radiance[k], skybr * D6415[k], cloud_cover);
304	<	ground_radiance[k] =
305	<	wmean2(ground_radiance[k], grndbr * D6415[k], cloud_cover);
298	>	double skybr = get_overcast_brightness(rdir[2], overcast_zenithbr);
299	>	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		}
309
310	<	scolor2scolr(sky_sclr, sky_radiance, 20);
311	<	putbinary(sky_sclr, LSCOLR, 1, skyfp);
310	>	for (k = 0; k < NSSAMP; ++k) {
311	>	radiance[k] *= dif_ratio;
312	>	}
313
314	<	scolor2scolr(ground_sclr, ground_radiance, 20);
315	<	putbinary(ground_sclr, LSCOLR, 1, grndfp);
314	>	scolor2scolr(sky_sclr, radiance, NSSAMP);
315	>	putbinary(sky_sclr, LSCOLR, 1, skyfp);
316		}
317		}
318		fclose(skyfp);
316	–	fclose(grndfp);
319
320		/* 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 zenithbr = get_zenith_brightness(sundir);
326	<	double skybr = get_overcast_brightness(sundir[2], zenithbr);
327	<	for (int i = 0; i < NSSAMP; ++i) {
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		}
332
333	<	write_rad(sun_radiance, sundir, skyfile, grndfile);
333	>	/* Normalize */
334	>	double sum = 0.0;
335	>	for (i = 0; i < NSSAMP; ++i) {
336	>	sum += sun_radiance[i];
337	>	}
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	>	}
346	>
347	>	write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
348		return 1;
349		}
350
#	Line 351 | Line 367 | static DpPaths get_dppaths(const char *dir, const doub
367		static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
368		const int is_summer,
369		const double s_latitude) {
354	–	/* Set rayleigh density profile */
370		if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
371		tag[0] = 's';
372		if (is_summer) {
#	Line 413 | Line 428 | static Atmosphere init_atmos(const double aod, const d
428		}
429
430		int main(int argc, char *argv[]) {
416	–	progname = argv[0];
431		int month, day;
432		double hour;
433		FVECT sundir;
#	Line 424 | Line 438 | int main(int argc, char *argv[]) {
438		int got_meridian = 0;
439		double grefl = 0.2;
440		double ccover = 0.0;
441	<	int res = 128;
441	>	int res = 64;
442		double aod = AOD0_CA;
443		char *outname = "out";
444		char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
445		char mie_name[20] = "mie_ca";
446		char lstag[3];
447		char *ddir = ".";
448	+	int i;
449	+	double dirnorm = 0; /* direct normal illuminance */
450	+	double difhor = 0;  /* diffuse horizontal illuminance */
451
452	<	if (!strcmp(argv[1], "-defaults")) {
452	>	fixargv0(argv[0]);
453	>	if (argc == 2 && !strcmp(argv[1], "-defaults")) {
454		printf("-i %d\t\t\t\t#scattering order\n", sorder);
455		printf("-g %f\t\t\t#ground reflectance\n", grefl);
456		printf("-c %f\t\t\t#cloud cover\n", ccover);
#	Line 440 | Line 458 | int main(int argc, char *argv[]) {
458		printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
459		printf("-f %s\t\t\t\t#output name (-f)\n", outname);
460		printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
461	<	exit(1);
461	>	exit(0);
462		}
463
464		if (argc < 4) {
465		fprintf(stderr,
466		"Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r "
467	<	"res -n nproc -c ccover -l mie -g grefl -f outpath\n",
467	>	"res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
468	>	"-g grefl -f outpath\n",
469		argv[0]);
470		return 0;
471		}
#	Line 468 | Line 487 | int main(int argc, char *argv[]) {
487		exit(1);
488		}
489
490	<	for (int i = 4; i < argc; i++) {
490	>	for (i = 4; i < argc; i++) {
491		if (argv[i][0] == '-') {
492		switch (argv[i][1]) {
493		case 'a':
#	Line 506 | Line 525 | int main(int argc, char *argv[]) {
525		case 'o':
526		s_longitude = atof(argv[++i]) * (PI / 180.0);
527		break;
528	+	case 'L':
529	+	dirnorm = atof(argv[++i]);
530	+	difhor = atof(argv[++i]);
531	+	break;
532		case 'p':
533		ddir = argv[++i];
534		break;
#	Line 547 | Line 570 | int main(int argc, char *argv[]) {
570
571		char gsdir[PATH_MAX];
572		size_t siz = strlen(ddir);
573	<	if (ISDIRSEP(ddir[siz-1]))
574	<	ddir[siz-1] = '\0';
573	>	if (ISDIRSEP(ddir[siz - 1]))
574	>	ddir[siz - 1] = '\0';
575		snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
553	–	printf("gsdir: %s\n", gsdir);
576		if (!make_directory(gsdir)) {
577		fprintf(stderr, "Failed creating atmos_data directory");
578		exit(1);
#	Line 576 | Line 598 | int main(int argc, char *argv[]) {
598		write_header(argc, argv, ccover, grefl, res);
599
600		if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
601	<	irrad_clear_dp, ccover, sundir, grefl, res, outname)) {
601	>	irrad_clear_dp, ccover, sundir, grefl, res, outname, ddir,
602	>	dirnorm, difhor)) {
603		fprintf(stderr, "gen_spect_sky failed\n");
604		exit(1);
605		}

Diff Legend

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