--- ray/src/gen/gendaymtx.c 2019/06/25 00:09:45 2.27 +++ ray/src/gen/gendaymtx.c 2020/01/07 18:26:55 2.35 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: gendaymtx.c,v 2.27 2019/06/25 00:09:45 greg Exp $"; +static const char RCSid[] = "$Id: gendaymtx.c,v 2.35 2020/01/07 18:26:55 greg Exp $"; #endif /* * gendaymtx.c @@ -81,16 +81,13 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.27 2 /* Include files */ #define _USE_MATH_DEFINES -#include #include -#include #include +#include "platform.h" #include "rtmath.h" #include "rtio.h" -#include "resolu.h" -#include "platform.h" #include "color.h" -#include "resolu.h" +#include "sun.h" char *progname; /* Program name */ char errmsg[128]; /* Error message buffer */ @@ -255,15 +252,7 @@ static const ModelCoeff DirectLumEff[8] = #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; +#define SUN_ANG_DEG 0.533 /* sun full-angle in degrees */ int nsuns = NSUNPATCH; /* number of sun patches to use */ double fixed_sun_sa = -1; /* fixed solid angle per sun? */ @@ -294,6 +283,7 @@ float *rh_dom; /* sky patch solid angle (sr) */ extern int rh_init(void); extern float * resize_dmatrix(float *mtx_data, int nsteps, int npatch); +extern void OutputSun(int id, FILE *fp); extern void AddDirect(float *parr); @@ -320,13 +310,13 @@ main(int argc, char *argv[]) double rotation = 0; /* site rotation (degrees) */ double elevation; /* site elevation (meters) */ int dir_is_horiz; /* direct is meas. on horizontal? */ + FILE *sunsfp = NULL; /* output file for individual suns */ float *mtx_data = NULL; /* our matrix data */ - int ntsteps = 0; /* number of time steps */ int avgSky = 0; /* compute average sky r.t. matrix? */ - int nrows = 0; /* number of rows in matrix */ - int rows_alloc = 0; /* number of allocated rows */ + 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 */ @@ -377,24 +367,28 @@ main(int argc, char *argv[]) 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 'n': /* no matrix output */ + avgSky = -1; + rhsubdiv = 1; + /* fall through */ case 'd': /* solar (direct) only */ skycolor[0] = skycolor[1] = skycolor[2] = 0; - if (suncolor[1] <= 1e-4) { - inconsistent = 1; - suncolor[0] = suncolor[1] = suncolor[2] = 1; + grefl[0] = grefl[1] = grefl[2] = 0; + break; + case 'D': /* output suns to file */ + sunsfp = fopen(argv[++i], "w"); + if (!sunsfp) { + fprintf(stderr, "%s: cannot open '%s' for output\n", + progname, argv[i]); + exit(1); } break; case 's': /* sky only (no direct) */ suncolor[0] = suncolor[1] = suncolor[2] = 0; - 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') @@ -406,7 +400,7 @@ main(int argc, char *argv[]) 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]); + progname); exit(1); } fixed_sun_sa *= fixed_sun_sa*PI; @@ -419,9 +413,6 @@ main(int argc, char *argv[]) } 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]); @@ -469,8 +460,12 @@ main(int argc, char *argv[]) fprintf(stderr, "%s: location '%s'\n", progname, buf); fprintf(stderr, "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n", progname, s_latitude, s_longitude); - fprintf(stderr, "%s: %d sky patches per time step\n", - progname, nskypatch); + if (avgSky >= 0) + fprintf(stderr, "%s: %d sky patches\n", + progname, nskypatch); + if (sunsfp) + fprintf(stderr, "%s: outputting suns to file\n", + progname); if (rotation != 0) fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation); @@ -480,27 +475,24 @@ main(int argc, char *argv[]) s_longitude = DegToRad(s_longitude); s_meridian = DegToRad(s_meridian); /* initial allocation */ - mtx_data = resize_dmatrix(mtx_data, rows_alloc=2, nskypatch); + 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; - mtx_offset = 3*nskypatch*nrows; - nrows += !avgSky | !nrows; + mtx_offset = 3*nskypatch*nstored; + nstored += !avgSky | !nstored; /* make space for next row */ - if (nrows > rows_alloc) { - rows_alloc += (rows_alloc>>1) + nrows + 7; - mtx_data = resize_dmatrix(mtx_data, rows_alloc, nskypatch); + 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) { - if (!avgSky | !nrows) + if (!avgSky | !mtx_offset) memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch); continue; } - if (verbose && mo != last_monthly) - fprintf(stderr, "%s: stepping through month %d...\n", - progname, last_monthly=mo); /* compute solar position */ julian_date = jdate(mo, da); sda = sdec(julian_date); @@ -519,11 +511,26 @@ main(int argc, char *argv[]) } /* compute sky patch values */ ComputeSky(mtx_data+mtx_offset); + + if (sunsfp) /* output sun if indicated */ + OutputSun(ntsteps, sunsfp); + + if (avgSky < 0) /* no matrix? */ + continue; + AddDirect(mtx_data+mtx_offset); /* update cumulative sky? */ for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; ) mtx_data[i] += mtx_data[mtx_offset+i]; + /* monthly reporting */ + if (verbose && mo != last_monthly) + fprintf(stderr, "%s: stepping through month %d...\n", + progname, last_monthly=mo); } + if (!ntsteps) { + fprintf(stderr, "%s: no valid time steps on input\n", progname); + exit(1); + } /* check for junk at end */ while ((i = fgetc(stdin)) != EOF) if (!isspace(i)) { @@ -534,10 +541,10 @@ 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); - } + + if (avgSky < 0) /* no matrix output? */ + goto alldone; + dif = 1./(double)ntsteps; /* average sky? */ for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; ) mtx_data[i] *= dif; @@ -549,15 +556,17 @@ main(int argc, char *argv[]) #endif if (verbose) fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", - progname, outfmt=='a' ? "" : "binary ", nrows); + 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", nrows); + printf("NCOLS=%d\n", nstored); printf("NCOMP=3\n"); + if ((outfmt == 'f') | (outfmt == 'd')) + fputendian(stdout); fputformat((char *)getfmtname(outfmt), stdout); putchar('\n'); } @@ -566,24 +575,24 @@ main(int argc, char *argv[]) mtx_offset = 3*i; switch (outfmt) { case 'a': - for (j = 0; j < nrows; 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 (nrows > 1) + if (nstored > 1) fputc('\n', stdout); break; case 'f': - for (j = 0; j < nrows; j++) { + 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 < nrows; j++) { + for (j = 0; j < nstored; j++) { double ment[3]; ment[0] = mtx_data[mtx_offset]; ment[1] = mtx_data[mtx_offset+1]; @@ -596,17 +605,18 @@ main(int argc, char *argv[]) if (ferror(stdout)) goto writerr; } - if (fflush(stdout) == EOF) +alldone: + if (fflush(NULL) == EOF) goto writerr; if (verbose) fprintf(stderr, "%s: done.\n", progname); exit(0); userr: - 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", + fprintf(stderr, "Usage: %s [-v][-h][-A][-d|-s|-n][-D file][-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: - fprintf(stderr, "%s: input weather tape format error\n", progname); + fprintf(stderr, "%s: weather tape format error in header\n", progname); exit(1); writerr: fprintf(stderr, "%s: write error on output\n", progname); @@ -678,11 +688,6 @@ ComputeSky(float *parr) 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; - } /* Compute ground radiance (include solar contribution if any) */ parr[0] = diff_illum; if (altitude > 0) @@ -690,6 +695,10 @@ ComputeSky(float *parr) parr[2] = parr[1] = parr[0] *= (1./PI/WHTEFFICACY); multcolor(parr, grefl); + if (bright(skycolor) <= 1e-4) { /* 0 sky component? */ + memset(parr+3, 0, sizeof(float)*3*(nskypatch-1)); + return; + } /* Calculate Perez sky model parameters */ CalcPerezParam(sun_zenith, sky_clearness, sky_brightness, index); @@ -765,6 +774,26 @@ AddDirect(float *parr) *pdest++ += val_add*suncolor[1]; *pdest++ += val_add*suncolor[2]; } +} + +/* Output a sun to indicated file if appropriate for this time step */ +void +OutputSun(int id, FILE *fp) +{ + double srad; + FVECT sv; + + if ((altitude <= 0) | (dir_illum <= 1e-4)) + return; + + srad = DegToRad(SUN_ANG_DEG/2.); + srad = dir_illum/(WHTEFFICACY * PI*srad*srad); + vector(sv, altitude, azimuth); + fprintf(fp, "\nvoid light solar%d\n0\n0\n", id); + fprintf(fp, "3 %.3e %.3e %.3e\n", srad*suncolor[0], + srad*suncolor[1], srad*suncolor[2]); + fprintf(fp, "\nsolar%d source sun%d\n0\n0\n", id, id); + fprintf(fp, "4 %.6f %.6f %.6f %.4f\n", sv[0], sv[1], sv[2], SUN_ANG_DEG); } /* Initialize Reinhart sky patch positions (GW) */