ViewVC Help

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

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

Comparing ray/src/gen/gendaymtx.c (file contents):
Revision 2.2 by greg, Fri Jan 18 19:56:03 2013 UTC vs.
Revision 2.15 by greg, Tue Jun 17 18:51:47 2014 UTC

#	Line 86 | Line 86 | static const char RCSid[] = "$Id$";
86		#include <string.h>
87		#include <ctype.h>
88		#include "rtmath.h"
89	+	#include "platform.h"
90		#include "color.h"
91
92		char *progname;                                                         /* Program name */
#	Line 108 | Line 109 | double sky_clearness;                  /* Sky clearness */
109		double solar_rad;                       /* Solar radiance */
110		double sun_zenith;                      /* Sun zenith angle (radians) */
111		int     input = 0;                              /* Input type */
112	+	int     output = 0;                             /* Output type */
113
114		extern double dmax( double, double );
115		extern double CalcAirMass();
#	Line 208 | Line 210 | static const CategoryBounds SkyClearCat[8] =
210		{ 1.950, 2.800 },
211		{ 2.800, 4.500 },
212		{ 4.500, 6.200 },
213	<	{ 6.200, 12.00 }        /* Clear */
213	>	{ 6.200, 12.01 }        /* Clear */
214		};
215
216		/* Luminous efficacy model coefficients */
#	Line 246 | Line 248 | static const ModelCoeff DirectLumEff[8] =
248		{ 101.18,  1.58, -1.10,  -8.29 }
249		};
250
251	+	#ifndef NSUNPATCH
252	+	#define NSUNPATCH       4               /* max. # patches to spread sun into */
253	+	#endif
254	+
255		extern int jdate(int month, int day);
256		extern double stadj(int  jd);
257		extern double sdec(int  jd);
#	Line 256 | Line 262 | extern double  s_latitude;
262		extern double  s_longitude;
263		extern double  s_meridian;
264
265	<	double          grefl = 0.2;            /* diffuse ground reflectance */
265	>	int             nsuns = NSUNPATCH;      /* number of sun patches to use */
266	>	double          fixed_sun_sa = -1;      /* fixed solid angle per sun? */
267
268		int             verbose = 0;            /* progress reports to stderr? */
269
#	Line 264 | Line 271 | int            outfmt = 'a';           /* output format */
271
272		int             rhsubdiv = 1;           /* Reinhart sky subdivisions */
273
274	<	float           skycolor[3] = {.96, 1.004, 1.118};      /* sky coloration */
274	>	COLOR           skycolor = {.96, 1.004, 1.118}; /* sky coloration */
275	>	COLOR           suncolor = {1., 1., 1.};        /* sun color */
276	>	COLOR           grefl = {.2, .2, .2};           /* ground reflectance */
277
269	–	int             do_sun = 1;             /* output direct solar contribution? */
270	–
278		int             nskypatch;              /* number of Reinhart patches */
279		float           *rh_palt;               /* sky patch altitudes (radians) */
280		float           *rh_pazi;               /* sky patch azimuths (radians) */
#	Line 286 | Line 293 | extern int     rh_init(void);
293		extern float *  resize_dmatrix(float *mtx_data, int nsteps, int npatch);
294		extern void     AddDirect(float *parr);
295
296	+
297	+	static const char *
298	+	getfmtname(int fmt)
299	+	{
300	+	switch (fmt) {
301	+	case 'a':
302	+	return("ascii");
303	+	case 'f':
304	+	return("float");
305	+	case 'd':
306	+	return("double");
307	+	}
308	+	return("unknown");
309	+	}
310	+
311	+
312		int
313		main(int argc, char *argv[])
314		{
315		char    buf[256];
316	+	int     doheader = 1;           /* output header? */
317	+	double  rotation = 0;           /* site rotation (degrees) */
318		double  elevation;              /* site elevation (meters) */
319		int     dir_is_horiz;           /* direct is meas. on horizontal? */
320		float   *mtx_data = NULL;       /* our matrix data */
#	Line 305 | Line 330 | main(int argc, char *argv[])
330		/* get options */
331		for (i = 1; i < argc && argv[i][0] == '-'; i++)
332		switch (argv[i][1]) {
333	<	case 'g':
334	<	grefl = atof(argv[++i]);
333	>	case 'g':                       /* ground reflectance */
334	>	grefl[0] = atof(argv[++i]);
335	>	grefl[1] = atof(argv[++i]);
336	>	grefl[2] = atof(argv[++i]);
337		break;
338	<	case 'v':
338	>	case 'v':                       /* verbose progress reports */
339		verbose++;
340		break;
341	<	case 'o':
341	>	case 'h':                       /* turn off header */
342	>	doheader = 0;
343	>	break;
344	>	case 'o':                       /* output format */
345		switch (argv[i][2]) {
346		case 'f':
347		case 'd':
#	Line 322 | Line 352 | main(int argc, char *argv[])
352		goto userr;
353		}
354		break;
355	<	case 'm':
355	>	case 'O':                       /* output type */
356	>	switch (argv[i][2]) {
357	>	case '0':
358	>	output = 0;
359	>	break;
360	>	case '1':
361	>	output = 1;
362	>	break;
363	>	default:
364	>	goto userr;
365	>	}
366	>	if (argv[i][3])
367	>	goto userr;
368	>	break;
369	>	case 'm':                       /* Reinhart subdivisions */
370		rhsubdiv = atoi(argv[++i]);
371		break;
372	<	case 'c':
372	>	case 'c':                       /* sky color */
373		skycolor[0] = atof(argv[++i]);
374		skycolor[1] = atof(argv[++i]);
375		skycolor[2] = atof(argv[++i]);
376		break;
377	<	case 'd':
334	<	do_sun = 1;
377	>	case 'd':                       /* solar (direct) only */
378		skycolor[0] = skycolor[1] = skycolor[2] = 0;
379	+	if (suncolor[1] <= 1e-4)
380	+	suncolor[0] = suncolor[1] = suncolor[2] = 1;
381		break;
382	<	case 's':
383	<	do_sun = 0;
382	>	case 's':                       /* sky only (no direct) */
383	>	suncolor[0] = suncolor[1] = suncolor[2] = 0;
384		if (skycolor[1] <= 1e-4)
385		skycolor[0] = skycolor[1] = skycolor[2] = 1;
386		break;
387	+	case 'r':                       /* rotate distribution */
388	+	if (argv[i][2] && argv[i][2] != 'z')
389	+	goto userr;
390	+	rotation = atof(argv[++i]);
391	+	break;
392	+	case '5':                       /* 5-phase calculation */
393	+	nsuns = 1;
394	+	fixed_sun_sa = 6.797e-05;
395	+	break;
396		default:
397		goto userr;
398		}
#	Line 393 | Line 447 | main(int argc, char *argv[])
447		progname, s_latitude, s_longitude);
448		fprintf(stderr, "%s: %d sky patches per time step\n",
449		progname, nskypatch);
450	+	if (rotation != 0)
451	+	fprintf(stderr, "%s: rotating output %.0f degrees\n",
452	+	progname, rotation);
453		}
454		/* convert quantities to radians */
455		s_latitude = DegToRad(s_latitude);
#	Line 400 | Line 457 | main(int argc, char *argv[])
457		s_meridian = DegToRad(s_meridian);
458		/* process each time step in tape */
459		while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) {
460	+	static int      step_alloc = 0;
461		double          sda, sta;
462		/* make space for next time step */
463		mtx_offset = 3*nskypatch*ntsteps++;
464	<	mtx_data = resize_dmatrix(mtx_data, ntsteps, nskypatch);
464	>	if (ntsteps > step_alloc) {
465	>	step_alloc += step_alloc>>1 + 8;
466	>	mtx_data = resize_dmatrix(mtx_data, step_alloc, nskypatch);
467	>	}
468		if (dif <= 1e-4) {
469		memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch);
470		continue;
#	Line 416 | Line 477 | main(int argc, char *argv[])
477		sda = sdec(julian_date);
478		sta = stadj(julian_date);
479		altitude = salt(sda, hr+sta);
480	<	azimuth = sazi(sda, hr+sta) + PI;
480	>	azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation);
481		/* convert measured values */
482		if (dir_is_horiz && altitude > 0.)
483		dir /= sin(altitude);
#	Line 429 | Line 490 | main(int argc, char *argv[])
490		}
491		/* compute sky patch values */
492		ComputeSky(mtx_data+mtx_offset);
493	<	if (do_sun)
433	<	AddDirect(mtx_data+mtx_offset);
493	>	AddDirect(mtx_data+mtx_offset);
494		}
495		/* check for junk at end */
496		while ((i = fgetc(stdin)) != EOF)
#	Line 443 | Line 503 | main(int argc, char *argv[])
503		break;
504		}
505		/* write out matrix */
506	+	if (outfmt != 'a')
507	+	SET_FILE_BINARY(stdout);
508		#ifdef getc_unlocked
509		flockfile(stdout);
510		#endif
511		if (verbose)
512		fprintf(stderr, "%s: writing %smatrix with %d time steps...\n",
513		progname, outfmt=='a' ? "" : "binary ", ntsteps);
514	+	if (doheader) {
515	+	newheader("RADIANCE", stdout);
516	+	printargs(argc, argv, stdout);
517	+	printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude),
518	+	-RadToDeg(s_longitude));
519	+	printf("NROWS=%d\n", nskypatch);
520	+	printf("NCOLS=%d\n", ntsteps);
521	+	printf("NCOMP=3\n");
522	+	fputformat(getfmtname(outfmt), stdout);
523	+	putchar('\n');
524	+	}
525		/* patches are rows (outer sort) */
526		for (i = 0; i < nskypatch; i++) {
527		mtx_offset = 3*i;
528		switch (outfmt) {
529		case 'a':
530		for (j = 0; j < ntsteps; j++) {
531	<	printf("%.3e %.3e %.3e\n", mtx_data[mtx_offset],
531	>	printf("%.3g %.3g %.3g\n", mtx_data[mtx_offset],
532		mtx_data[mtx_offset+1],
533		mtx_data[mtx_offset+2]);
534		mtx_offset += 3*nskypatch;
#	Line 490 | Line 563 | main(int argc, char *argv[])
563		fprintf(stderr, "%s: done.\n", progname);
564		exit(0);
565		userr:
566	<	fprintf(stderr, "Usage: %s [-v][-d|-s][-m N][-g refl][-c r g b][-o{f|d}] [tape.wea]\n",
566	>	fprintf(stderr, "Usage: %s [-v][-h][-d|-s][-r deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
567		progname);
568		exit(1);
569		fmterr:
#	Line 511 | Line 584 | ComputeSky(float *parr)
584		{
585		int index;                      /* Category index */
586		double norm_diff_illum;         /* Normalized diffuse illuimnance */
514	–	double zlumin;                  /* Zenith luminance */
587		int i;
588
589		/* Calculate atmospheric precipitable water content */
590		apwc = CalcPrecipWater(dew_point);
591
592	<	/* Limit solar altitude to keep circumsolar off zenith */
593	<	if (altitude > DegToRad(87.0))
594	<	altitude = DegToRad(87.0);
592	>	/* Calculate sun zenith angle (don't let it dip below horizon) */
593	>	/* Also limit minimum angle to keep circumsolar off zenith */
594	>	if (altitude <= 0.0)
595	>	sun_zenith = DegToRad(90.0);
596	>	else if (altitude >= DegToRad(87.0))
597	>	sun_zenith = DegToRad(3.0);
598	>	else
599	>	sun_zenith = DegToRad(90.0) - altitude;
600
524	–	/* Calculate sun zenith angle */
525	–	sun_zenith = DegToRad(90.0) - altitude;
526	–
601		/* Compute the inputs for the calculation of the sky distribution */
602
603		if (input == 0)                                 /* XXX never used */
#	Line 542 | Line 616 | ComputeSky(float *parr)
616		sky_brightness = CalcSkyBrightness();
617		sky_clearness =  CalcSkyClearness();
618
619	+	/* Limit sky clearness */
620	+	if (sky_clearness > 11.9)
621	+	sky_clearness = 11.9;
622	+
623	+	/* Limit sky brightness */
624	+	if (sky_brightness < 0.01)
625	+	sky_brightness = 0.01;
626	+
627		/* Calculate illuminance */
628		index = GetCategoryIndex();
629		diff_illum = diff_irrad * CalcDiffuseIllumRatio(index);
#	Line 553 | Line 635 | ComputeSky(float *parr)
635		index = CalcSkyParamFromIllum();
636		}
637
638	+	if (output == 1) {                      /* hack for solar radiance */
639	+	diff_illum = diff_irrad * WHTEFFICACY;
640	+	dir_illum = dir_irrad * WHTEFFICACY;
641	+	}
642	+
643		if (bright(skycolor) <= 1e-4) {                 /* 0 sky component? */
644		memset(parr, 0, sizeof(float)*3*nskypatch);
645		return;
646		}
647		/* Compute ground radiance (include solar contribution if any) */
648	<	parr[0] = diff_illum * (1./PI/WHTEFFICACY);
648	>	parr[0] = diff_illum;
649		if (altitude > 0)
650	<	parr[0] += dir_illum * sin(altitude) * (1./PI/WHTEFFICACY);
651	<	parr[2] = parr[1] = parr[0];
650	>	parr[0] += dir_illum * sin(altitude);
651	>	parr[2] = parr[1] = parr[0] *= (1./PI/WHTEFFICACY);
652	>	multcolor(parr, grefl);
653
654		/* Calculate Perez sky model parameters */
655		CalcPerezParam(sun_zenith, sky_clearness, sky_brightness, index);
#	Line 572 | Line 660 | ComputeSky(float *parr)
660		/* Calculate relative horizontal illuminance */
661		norm_diff_illum = CalcRelHorzIllum(parr);
662
663	+	/* Check for zero sky -- make uniform in that case */
664	+	if (norm_diff_illum <= FTINY) {
665	+	for (i = 1; i < nskypatch; i++)
666	+	setcolor(parr+3*i, 1., 1., 1.);
667	+	norm_diff_illum = PI;
668	+	}
669		/* Normalization coefficient */
670		norm_diff_illum = diff_illum / norm_diff_illum;
671
578	–	/* Calculate relative zenith luminance */
579	–	zlumin = CalcRelLuminance(sun_zenith, 0.0);
580	–
581	–	/* Calculate absolute zenith illuminance */
582	–	zlumin *= norm_diff_illum;
583	–
672		/* Apply to sky patches to get absolute radiance values */
673		for (i = 1; i < nskypatch; i++) {
674	<	scalecolor(parr+3*i, zlumin*(1./WHTEFFICACY));
674	>	scalecolor(parr+3*i, norm_diff_illum*(1./WHTEFFICACY));
675		multcolor(parr+3*i, skycolor);
676		}
677		}
#	Line 593 | Line 681 | void
681		AddDirect(float *parr)
682		{
683		FVECT   svec;
684	<	double  near_dprod[4];
685	<	int     near_patch[4];
686	<	double  wta[4], wtot;
684	>	double  near_dprod[NSUNPATCH];
685	>	int     near_patch[NSUNPATCH];
686	>	double  wta[NSUNPATCH], wtot;
687		int     i, j, p;
688
689	<	if (!do_sun || dir_illum < 1e-4)
689	>	if (dir_illum <= 1e-4 || bright(suncolor) <= 1e-4)
690		return;
691	<	/* identify 4 closest patches */
692	<	for (i = 4; i--; )
691	>	/* identify nsuns closest patches */
692	>	if (nsuns > NSUNPATCH)
693	>	nsuns = NSUNPATCH;
694	>	else if (nsuns <= 0)
695	>	nsuns = 1;
696	>	for (i = nsuns; i--; )
697		near_dprod[i] = -1.;
698		vector(svec, altitude, azimuth);
699		for (p = 1; p < nskypatch; p++) {
#	Line 609 | Line 701 | AddDirect(float *parr)
701		double  dprod;
702		rh_vector(pvec, p);
703		dprod = DOT(pvec, svec);
704	<	for (i = 0; i < 4; i++)
704	>	for (i = 0; i < nsuns; i++)
705		if (dprod > near_dprod[i]) {
706	<	for (j = 4; --j > i; ) {
706	>	for (j = nsuns; --j > i; ) {
707		near_dprod[j] = near_dprod[j-1];
708		near_patch[j] = near_patch[j-1];
709		}
#	Line 621 | Line 713 | AddDirect(float *parr)
713		}
714		}
715		wtot = 0;                       /* weight by proximity */
716	<	for (i = 4; i--; )
716	>	for (i = nsuns; i--; )
717		wtot += wta[i] = 1./(1.002 - near_dprod[i]);
718		/* add to nearest patch radiances */
719	<	for (i = 4; i--; ) {
719	>	for (i = nsuns; i--; ) {
720		float   *pdest = parr + 3*near_patch[i];
721	<	float   val_add = wta[i] * dir_illum /
722	<	(WHTEFFICACY * wtot * rh_dom[near_patch[i]]);
723	<	*pdest++ += val_add;
724	<	*pdest++ += val_add;
725	<	*pdest++ += val_add;
721	>	float   val_add = wta[i] * dir_illum / (WHTEFFICACY * wtot);
722	>
723	>	val_add /= (fixed_sun_sa > 0)   ? fixed_sun_sa
724	>	: rh_dom[near_patch[i]] ;
725	>	*pdest++ += val_add*suncolor[0];
726	>	*pdest++ += val_add*suncolor[1];
727	>	*pdest++ += val_add*suncolor[2];
728		}
729		}
730
#	Line 665 | Line 759 | rh_init(void)
759		for (i = 0; i < NROW*rhsubdiv; i++) {
760		const float     ralt = alpha*(i + .5);
761		const int       ninrow = tnaz[i/rhsubdiv]*rhsubdiv;
762	<	const float     dom = (sin(alpha*(i+1)) - sin(alpha*i))/ninrow;
762	>	const float     dom = 2.*PI*(sin(alpha*(i+1)) - sin(alpha*i)) /
763	>	(double)ninrow;
764		for (j = 0; j < ninrow; j++) {
765		rh_palt[p] = ralt;
766		rh_pazi[p] = 2.*PI * j / (double)ninrow;
#	Line 770 | Line 865 | double CalcSkyClearness()
865		double sz_cubed;        /* Sun zenith angle cubed */
866
867		/* Calculate sun zenith angle cubed */
868	<	sz_cubed = pow(sun_zenith, 3.0);
868	>	sz_cubed = sun_zenith*sun_zenith*sun_zenith;
869
870		return ((diff_irrad + dir_irrad) / diff_irrad + 1.041 *
871		sz_cubed) / (1.0 + 1.041 * sz_cubed);
#	Line 801 | Line 896 | double CalcDiffuseIrradiance()
896		double CalcDirectIrradiance()
897		{
898		return CalcDiffuseIrradiance() * ((sky_clearness - 1.0) * (1 + 1.041
899	<	* pow(sun_zenith, 3.0)));
899	>	* sun_zenith*sun_zenith*sun_zenith));
900		}
901
902		/* Calculate sky brightness and clearness from illuminance values */
#	Line 827 | Line 922 | int CalcSkyParamFromIllum()
922		sky_clearness = 12.0;
923
924		/* Limit sky brightness */
925	<	if (sky_brightness < 0.05)
925	>	if (sky_brightness < 0.01)
926		sky_brightness = 0.01;
927
928		while (((fabs(diff_irrad - test1) > 10.0) ||
#	Line 851 | Line 946 | int CalcSkyParamFromIllum()
946		sky_clearness = 12.0;
947
948		/* Limit sky brightness */
949	<	if (sky_brightness < 0.05)
949	>	if (sky_brightness < 0.01)
950		sky_brightness = 0.01;
951		}
952
#	Line 937 | Line 1032 | double CalcRelHorzIllum( float *parr )
1032		double rh_illum = 0.0;  /* Relative horizontal illuminance */
1033
1034		for (i = 1; i < nskypatch; i++)
1035	<	rh_illum += parr[3*i+1] * rh_cos(i);
1035	>	rh_illum += parr[3*i+1] * rh_cos(i) * rh_dom[i];
1036
1037	<	return rh_illum * (2.0 * PI / (nskypatch-1));
1037	>	return rh_illum;
1038		}
1039
1040		/* Calculate earth orbit eccentricity correction factor */

Diff Legend

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