--- ray/src/gen/gendaymtx.c 2020/01/07 18:26:55 2.35 +++ ray/src/gen/gendaymtx.c 2022/03/11 18:04:39 2.39 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: gendaymtx.c,v 2.35 2020/01/07 18:26:55 greg Exp $"; +static const char RCSid[] = "$Id: gendaymtx.c,v 2.39 2022/03/11 18:04:39 greg Exp $"; #endif /* * gendaymtx.c @@ -89,8 +89,7 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.35 2 #include "color.h" #include "sun.h" -char *progname; /* Program name */ -char errmsg[128]; /* Error message buffer */ +char *progname; /* Program name */ const double DC_SolarConstantE = 1367.0; /* Solar constant W/m^2 */ const double DC_SolarConstantL = 127.5; /* Solar constant klux */ @@ -135,7 +134,6 @@ extern void ComputeSky( float *parr ); /* Radiuans into degrees */ #define RadToDeg(rad) ((rad)*(180./PI)) - /* Perez sky model coefficients */ /* Reference: Perez, R., R. Seals, and J. Michalsky, 1993. "All- */ @@ -272,18 +270,20 @@ float *rh_palt; /* sky patch altitudes (radians) */ float *rh_pazi; /* sky patch azimuths (radians) */ float *rh_dom; /* sky patch solid angle (sr) */ -#define vector(v,alt,azi) ( (v)[1] = tcos(alt), \ - (v)[0] = (v)[1]*tsin(azi), \ - (v)[1] *= tcos(azi), \ - (v)[2] = tsin(alt) ) +#define vector(v,alt,azi) ( (v)[1] = cos(alt), \ + (v)[0] = (v)[1]*sin(azi), \ + (v)[1] *= cos(azi), \ + (v)[2] = sin(alt) ) #define rh_vector(v,i) vector(v,rh_palt[i],rh_pazi[i]) #define rh_cos(i) tsin(rh_palt[i]) +#define solar_minute(jd,hr) ((24*60)*((jd)-1)+(int)((hr)*60.+.5)) + 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 OutputSun(int id, int goodsun, FILE *fp, FILE *mfp); extern void AddDirect(float *parr); @@ -309,8 +309,11 @@ main(int argc, char *argv[]) int doheader = 1; /* output header? */ double rotation = 0; /* site rotation (degrees) */ double elevation; /* site elevation (meters) */ + int leap_day = 0; /* add leap day? */ + int sun_hours_only = 0; /* only output sun hours? */ int dir_is_horiz; /* direct is meas. on horizontal? */ FILE *sunsfp = NULL; /* output file for individual suns */ + FILE *modsfp = NULL; /* modifier output file */ float *mtx_data = NULL; /* our matrix data */ int avgSky = 0; /* compute average sky r.t. matrix? */ int ntsteps = 0; /* number of time steps */ @@ -371,6 +374,19 @@ main(int argc, char *argv[]) skycolor[1] = atof(argv[++i]); skycolor[2] = atof(argv[++i]); break; + case 'D': /* output suns to file */ + if (strcmp(argv[++i], "-")) { + sunsfp = fopen(argv[i], "w"); + if (sunsfp == NULL) { + fprintf(stderr, + "%s: cannot open '%s' for output\n", + progname, argv[i]); + exit(1); + } + break; /* still may output matrix */ + } + sunsfp = stdout; /* sending to stdout, so... */ + /* fall through */ case 'n': /* no matrix output */ avgSky = -1; rhsubdiv = 1; @@ -379,9 +395,8 @@ main(int argc, char *argv[]) skycolor[0] = skycolor[1] = skycolor[2] = 0; grefl[0] = grefl[1] = grefl[2] = 0; break; - case 'D': /* output suns to file */ - sunsfp = fopen(argv[++i], "w"); - if (!sunsfp) { + case 'M': /* send sun modifiers to file */ + if ((modsfp = fopen(argv[++i], "w")) == NULL) { fprintf(stderr, "%s: cannot open '%s' for output\n", progname, argv[i]); exit(1); @@ -390,6 +405,9 @@ main(int argc, char *argv[]) case 's': /* sky only (no direct) */ suncolor[0] = suncolor[1] = suncolor[2] = 0; break; + case 'u': /* solar hours only */ + sun_hours_only = 1; + break; case 'r': /* rotate distribution */ if (argv[i][2] && argv[i][2] != 'z') goto userr; @@ -418,6 +436,9 @@ main(int argc, char *argv[]) progname, argv[i]); exit(1); } + if ((modsfp != NULL) & (sunsfp == NULL)) + fprintf(stderr, "%s: warning -M output will be empty without -D\n", + progname); if (verbose) { if (i == argc-1) fprintf(stderr, "%s: reading weather tape '%s'\n", @@ -464,8 +485,8 @@ main(int argc, char *argv[]) fprintf(stderr, "%s: %d sky patches\n", progname, nskypatch); if (sunsfp) - fprintf(stderr, "%s: outputting suns to file\n", - progname); + fprintf(stderr, "%s: outputting suns to %s\n", + progname, sunsfp==stdout ? "stdout" : "file"); if (rotation != 0) fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation); @@ -479,6 +500,20 @@ main(int argc, char *argv[]) /* process each time step in tape */ while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) { double sda, sta; + int sun_in_sky; + /* compute solar position */ + if ((mo == 2) & (da == 29)) { + julian_date = 60; + leap_day = 1; + } else + julian_date = jdate(mo, da) + leap_day; + sda = sdec(julian_date); + sta = stadj(julian_date); + altitude = salt(sda, hr+sta); + sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG/2.)); + if (sun_hours_only && !sun_in_sky) + continue; /* skipping nighttime points */ + azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); mtx_offset = 3*nskypatch*nstored; nstored += !avgSky | !nstored; @@ -488,19 +523,23 @@ main(int argc, char *argv[]) mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch); } ntsteps++; /* keep count of time steps */ - if (dif <= 1e-4) { + + if (dir+dif <= 1e-4) { /* effectively nighttime? */ if (!avgSky | !mtx_offset) - memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch); + memset(mtx_data+mtx_offset, 0, + sizeof(float)*3*nskypatch); + /* output black sun? */ + if (sunsfp && sun_in_sky) + OutputSun(solar_minute(julian_date,hr), 0, + sunsfp, modsfp); continue; } - /* compute solar position */ - julian_date = jdate(mo, da); - sda = sdec(julian_date); - sta = stadj(julian_date); - altitude = salt(sda, hr+sta); - azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); + if (!sun_in_sky && dir > (input==1 ? 20. : 20.*WHTEFFICACY)) + fprintf(stderr, + "%s: warning - unusually bright at %.1f on %d-%d\n", + progname, hr, mo, da); /* convert measured values */ - if (dir_is_horiz && altitude > 0.) + if (dir_is_horiz && altitude > FTINY) dir /= sin(altitude); if (input == 1) { dir_irrad = dir; @@ -511,10 +550,11 @@ main(int argc, char *argv[]) } /* compute sky patch values */ ComputeSky(mtx_data+mtx_offset); + /* output sun if requested */ + if (sunsfp && sun_in_sky) + OutputSun(solar_minute(julian_date,hr), 1, + sunsfp, modsfp); - if (sunsfp) /* output sun if indicated */ - OutputSun(ntsteps, sunsfp); - if (avgSky < 0) /* no matrix? */ continue; @@ -526,6 +566,7 @@ main(int argc, char *argv[]) if (verbose && mo != last_monthly) fprintf(stderr, "%s: stepping through month %d...\n", progname, last_monthly=mo); + /* note whether leap-day was given */ } if (!ntsteps) { fprintf(stderr, "%s: no valid time steps on input\n", progname); @@ -612,7 +653,7 @@ alldone: fprintf(stderr, "%s: done.\n", progname); exit(0); userr: - 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", + fprintf(stderr, "Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-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: @@ -778,22 +819,22 @@ AddDirect(float *parr) /* Output a sun to indicated file if appropriate for this time step */ void -OutputSun(int id, FILE *fp) +OutputSun(int id, int goodsun, FILE *fp, FILE *mfp) { 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); + srad = goodsun ? dir_illum/(WHTEFFICACY * PI*srad*srad) : 0; 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); + + if (mfp != NULL) /* saving modifier IDs? */ + fprintf(mfp, "solar%d\n", id); } /* Initialize Reinhart sky patch positions (GW) */