--- ray/src/gen/gendaymtx.c 2013/02/05 06:00:19 2.8 +++ ray/src/gen/gendaymtx.c 2019/11/07 23:15:07 2.30 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: gendaymtx.c,v 2.8 2013/02/05 06:00:19 greg Exp $"; +static const char RCSid[] = "$Id: gendaymtx.c,v 2.30 2019/11/07 23:15:07 greg Exp $"; #endif /* * gendaymtx.c @@ -85,8 +85,11 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.8 20 #include #include #include +#include "platform.h" #include "rtmath.h" +#include "rtio.h" #include "color.h" +#include "sun.h" char *progname; /* Program name */ char errmsg[128]; /* Error message buffer */ @@ -108,6 +111,7 @@ double sky_clearness; /* Sky clearness */ double solar_rad; /* Solar radiance */ double sun_zenith; /* Sun zenith angle (radians) */ int input = 0; /* Input type */ +int output = 0; /* Output type */ extern double dmax( double, double ); extern double CalcAirMass(); @@ -208,7 +212,7 @@ static const CategoryBounds SkyClearCat[8] = { 1.950, 2.800 }, { 2.800, 4.500 }, { 4.500, 6.200 }, - { 6.200, 12.00 } /* Clear */ + { 6.200, 12.01 } /* Clear */ }; /* Luminous efficacy model coefficients */ @@ -247,18 +251,11 @@ static const ModelCoeff DirectLumEff[8] = }; #ifndef NSUNPATCH -#define NSUNPATCH 4 /* # patches to spread sun into */ +#define NSUNPATCH 4 /* max. # patches to spread sun into */ #endif -extern int jdate(int month, int day); -extern double stadj(int jd); -extern double sdec(int jd); -extern double salt(double sd, double st); -extern double sazi(double sd, double st); - /* sun calculation constants */ -extern double s_latitude; -extern double s_longitude; -extern double s_meridian; +int nsuns = NSUNPATCH; /* number of sun patches to use */ +double fixed_sun_sa = -1; /* fixed solid angle per sun? */ int verbose = 0; /* progress reports to stderr? */ @@ -288,16 +285,37 @@ extern int rh_init(void); extern float * resize_dmatrix(float *mtx_data, int nsteps, int npatch); extern void AddDirect(float *parr); + +static const char * +getfmtname(int fmt) +{ + switch (fmt) { + case 'a': + return("ascii"); + case 'f': + return("float"); + case 'd': + return("double"); + } + return("unknown"); +} + + int main(int argc, char *argv[]) { char buf[256]; + int doheader = 1; /* output header? */ double rotation = 0; /* site rotation (degrees) */ double elevation; /* site elevation (meters) */ int dir_is_horiz; /* direct is meas. on horizontal? */ float *mtx_data = NULL; /* our matrix data */ - int ntsteps = 0; /* number of rows in matrix */ + int avgSky = 0; /* compute average sky r.t. matrix? */ + int ntsteps = 0; /* number of time steps */ + int tstorage = 0; /* number of allocated time steps */ + int nstored = 0; /* number of time steps in matrix */ int last_monthly = 0; /* month of last report */ + int inconsistent = 0; /* inconsistent options set? */ int mo, da; /* month (1-12) and day (1-31) */ double hr; /* hour (local standard time) */ double dir, dif; /* direct and diffuse values */ @@ -316,6 +334,9 @@ main(int argc, char *argv[]) case 'v': /* verbose progress reports */ verbose++; break; + case 'h': /* turn off header */ + doheader = 0; + break; case 'o': /* output format */ switch (argv[i][2]) { case 'f': @@ -327,34 +348,69 @@ main(int argc, char *argv[]) goto userr; } break; + case 'O': /* output type */ + switch (argv[i][2]) { + case '0': + output = 0; + break; + case '1': + output = 1; + break; + default: + goto userr; + } + if (argv[i][3]) + goto userr; + break; case 'm': /* Reinhart subdivisions */ rhsubdiv = atoi(argv[++i]); break; case 'c': /* sky color */ + inconsistent |= (skycolor[1] <= 1e-4); skycolor[0] = atof(argv[++i]); skycolor[1] = atof(argv[++i]); skycolor[2] = atof(argv[++i]); break; case 'd': /* solar (direct) only */ skycolor[0] = skycolor[1] = skycolor[2] = 0; - if (suncolor[1] <= 1e-4) + if (suncolor[1] <= 1e-4) { + inconsistent = 1; suncolor[0] = suncolor[1] = suncolor[2] = 1; + } break; case 's': /* sky only (no direct) */ suncolor[0] = suncolor[1] = suncolor[2] = 0; - if (skycolor[1] <= 1e-4) + if (skycolor[1] <= 1e-4) { + inconsistent = 1; skycolor[0] = skycolor[1] = skycolor[2] = 1; + } break; case 'r': /* rotate distribution */ if (argv[i][2] && argv[i][2] != 'z') goto userr; rotation = atof(argv[++i]); break; + case '5': /* 5-phase calculation */ + nsuns = 1; + fixed_sun_sa = PI/360.*atof(argv[++i]); + if (fixed_sun_sa <= 0) { + fprintf(stderr, "%s: missing solar disk size argument for '-5' option\n", + argv[0]); + exit(1); + } + fixed_sun_sa *= fixed_sun_sa*PI; + break; + case 'A': /* compute average sky */ + avgSky = 1; + break; default: goto userr; } if (i < argc-1) goto userr; + if (inconsistent) + fprintf(stderr, "%s: WARNING: inconsistent -s, -d, -c options!\n", + progname); if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) { fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]); @@ -412,14 +468,23 @@ main(int argc, char *argv[]) s_latitude = DegToRad(s_latitude); s_longitude = DegToRad(s_longitude); s_meridian = DegToRad(s_meridian); + /* initial allocation */ + mtx_data = resize_dmatrix(mtx_data, tstorage=2, nskypatch); /* process each time step in tape */ while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) { double sda, sta; - /* make space for next time step */ - mtx_offset = 3*nskypatch*ntsteps++; - mtx_data = resize_dmatrix(mtx_data, ntsteps, nskypatch); + + mtx_offset = 3*nskypatch*nstored; + nstored += !avgSky | !nstored; + /* make space for next row */ + if (nstored > tstorage) { + tstorage += (tstorage>>1) + nstored + 7; + mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch); + } + ntsteps++; /* keep count of time steps */ if (dif <= 1e-4) { - memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch); + if (!avgSky | !mtx_offset) + memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch); continue; } if (verbose && mo != last_monthly) @@ -444,6 +509,9 @@ main(int argc, char *argv[]) /* compute sky patch values */ ComputeSky(mtx_data+mtx_offset); AddDirect(mtx_data+mtx_offset); + /* update cumulative sky? */ + for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; ) + mtx_data[i] += mtx_data[mtx_offset+i]; } /* check for junk at end */ while ((i = fgetc(stdin)) != EOF) @@ -455,41 +523,63 @@ main(int argc, char *argv[]) fputs(buf, stderr); fputc('\n', stderr); break; } + if (!ntsteps) { + fprintf(stderr, "%s: no valid time steps on input\n", progname); + exit(1); + } + dif = 1./(double)ntsteps; /* average sky? */ + for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; ) + mtx_data[i] *= dif; /* write out matrix */ + if (outfmt != 'a') + SET_FILE_BINARY(stdout); #ifdef getc_unlocked flockfile(stdout); #endif if (verbose) fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", - progname, outfmt=='a' ? "" : "binary ", ntsteps); + progname, outfmt=='a' ? "" : "binary ", nstored); + if (doheader) { + newheader("RADIANCE", stdout); + printargs(argc, argv, stdout); + printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude), + -RadToDeg(s_longitude)); + printf("NROWS=%d\n", nskypatch); + printf("NCOLS=%d\n", nstored); + printf("NCOMP=3\n"); + if ((outfmt == 'f') | (outfmt == 'd')) + fputendian(stdout); + fputformat((char *)getfmtname(outfmt), stdout); + putchar('\n'); + } /* patches are rows (outer sort) */ for (i = 0; i < nskypatch; i++) { mtx_offset = 3*i; switch (outfmt) { case 'a': - for (j = 0; j < ntsteps; j++) { + for (j = 0; j < nstored; j++) { printf("%.3g %.3g %.3g\n", mtx_data[mtx_offset], mtx_data[mtx_offset+1], mtx_data[mtx_offset+2]); mtx_offset += 3*nskypatch; } - if (ntsteps > 1) + if (nstored > 1) fputc('\n', stdout); break; case 'f': - for (j = 0; j < ntsteps; j++) { - fwrite(mtx_data+mtx_offset, sizeof(float), 3, + for (j = 0; j < nstored; j++) { + putbinary(mtx_data+mtx_offset, sizeof(float), 3, stdout); mtx_offset += 3*nskypatch; } break; case 'd': - for (j = 0; j < ntsteps; j++) { + for (j = 0; j < nstored; j++) { double ment[3]; ment[0] = mtx_data[mtx_offset]; ment[1] = mtx_data[mtx_offset+1]; ment[2] = mtx_data[mtx_offset+2]; - fwrite(ment, sizeof(double), 3, stdout); + putbinary(ment, sizeof(double), 3, stdout); mtx_offset += 3*nskypatch; } break; @@ -503,7 +593,7 @@ main(int argc, char *argv[]) fprintf(stderr, "%s: done.\n", progname); exit(0); userr: - fprintf(stderr, "Usage: %s [-v][-d|-s][-r deg][-m N][-g r g b][-c r g b][-o{f|d}] [tape.wea]\n", + fprintf(stderr, "Usage: %s [-v][-h][-A][-d|-s][-r deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n", progname); exit(1); fmterr: @@ -556,6 +646,14 @@ ComputeSky(float *parr) sky_brightness = CalcSkyBrightness(); sky_clearness = CalcSkyClearness(); + /* Limit sky clearness */ + if (sky_clearness > 11.9) + sky_clearness = 11.9; + + /* Limit sky brightness */ + if (sky_brightness < 0.01) + sky_brightness = 0.01; + /* Calculate illuminance */ index = GetCategoryIndex(); diff_illum = diff_irrad * CalcDiffuseIllumRatio(index); @@ -567,6 +665,11 @@ ComputeSky(float *parr) index = CalcSkyParamFromIllum(); } + if (output == 1) { /* hack for solar radiance */ + diff_illum = diff_irrad * WHTEFFICACY; + dir_illum = dir_irrad * WHTEFFICACY; + } + if (bright(skycolor) <= 1e-4) { /* 0 sky component? */ memset(parr, 0, sizeof(float)*3*nskypatch); return; @@ -587,6 +690,12 @@ ComputeSky(float *parr) /* Calculate relative horizontal illuminance */ norm_diff_illum = CalcRelHorzIllum(parr); + /* Check for zero sky -- make uniform in that case */ + if (norm_diff_illum <= FTINY) { + for (i = 1; i < nskypatch; i++) + setcolor(parr+3*i, 1., 1., 1.); + norm_diff_illum = PI; + } /* Normalization coefficient */ norm_diff_illum = diff_illum / norm_diff_illum; @@ -609,8 +718,12 @@ AddDirect(float *parr) if (dir_illum <= 1e-4 || bright(suncolor) <= 1e-4) return; - /* identify NSUNPATCH closest patches */ - for (i = NSUNPATCH; i--; ) + /* identify nsuns closest patches */ + if (nsuns > NSUNPATCH) + nsuns = NSUNPATCH; + else if (nsuns <= 0) + nsuns = 1; + for (i = nsuns; i--; ) near_dprod[i] = -1.; vector(svec, altitude, azimuth); for (p = 1; p < nskypatch; p++) { @@ -618,9 +731,9 @@ AddDirect(float *parr) double dprod; rh_vector(pvec, p); dprod = DOT(pvec, svec); - for (i = 0; i < NSUNPATCH; i++) + for (i = 0; i < nsuns; i++) if (dprod > near_dprod[i]) { - for (j = NSUNPATCH; --j > i; ) { + for (j = nsuns; --j > i; ) { near_dprod[j] = near_dprod[j-1]; near_patch[j] = near_patch[j-1]; } @@ -630,13 +743,15 @@ AddDirect(float *parr) } } wtot = 0; /* weight by proximity */ - for (i = NSUNPATCH; i--; ) + for (i = nsuns; i--; ) wtot += wta[i] = 1./(1.002 - near_dprod[i]); /* add to nearest patch radiances */ - for (i = NSUNPATCH; i--; ) { + for (i = nsuns; i--; ) { float *pdest = parr + 3*near_patch[i]; - float val_add = wta[i] * dir_illum / - (WHTEFFICACY * wtot * rh_dom[near_patch[i]]); + float val_add = wta[i] * dir_illum / (WHTEFFICACY * wtot); + + val_add /= (fixed_sun_sa > 0) ? fixed_sun_sa + : rh_dom[near_patch[i]] ; *pdest++ += val_add*suncolor[0]; *pdest++ += val_add*suncolor[1]; *pdest++ += val_add*suncolor[2]; @@ -780,7 +895,7 @@ double CalcSkyClearness() double sz_cubed; /* Sun zenith angle cubed */ /* Calculate sun zenith angle cubed */ - sz_cubed = pow(sun_zenith, 3.0); + sz_cubed = sun_zenith*sun_zenith*sun_zenith; return ((diff_irrad + dir_irrad) / diff_irrad + 1.041 * sz_cubed) / (1.0 + 1.041 * sz_cubed); @@ -811,7 +926,7 @@ double CalcDiffuseIrradiance() double CalcDirectIrradiance() { return CalcDiffuseIrradiance() * ((sky_clearness - 1.0) * (1 + 1.041 - * pow(sun_zenith, 3.0))); + * sun_zenith*sun_zenith*sun_zenith)); } /* Calculate sky brightness and clearness from illuminance values */ @@ -837,7 +952,7 @@ int CalcSkyParamFromIllum() sky_clearness = 12.0; /* Limit sky brightness */ - if (sky_brightness < 0.05) + if (sky_brightness < 0.01) sky_brightness = 0.01; while (((fabs(diff_irrad - test1) > 10.0) || @@ -850,7 +965,9 @@ int CalcSkyParamFromIllum() /* Convert illuminance to irradiance */ index = GetCategoryIndex(); diff_irrad = diff_illum / CalcDiffuseIllumRatio(index); - dir_irrad = dir_illum / CalcDirectIllumRatio(index); + dir_irrad = CalcDirectIllumRatio(index); + if (dir_irrad > 0.1) + dir_irrad = dir_illum / dir_irrad; /* Calculate sky brightness and clearness */ sky_brightness = CalcSkyBrightness(); @@ -861,7 +978,7 @@ int CalcSkyParamFromIllum() sky_clearness = 12.0; /* Limit sky brightness */ - if (sky_brightness < 0.05) + if (sky_brightness < 0.01) sky_brightness = 0.01; }