--- ray/src/gen/gendaymtx.c 2016/08/20 15:42:38 2.25 +++ ray/src/gen/gendaymtx.c 2025/02/27 19:00:00 2.42 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: gendaymtx.c,v 2.25 2016/08/20 15:42:38 greg Exp $"; +static const char RCSid[] = "$Id: gendaymtx.c,v 2.42 2025/02/27 19:00:00 greg Exp $"; #endif /* * gendaymtx.c @@ -81,19 +81,16 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.25 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" +#include "loadEPW.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 */ @@ -132,13 +129,16 @@ extern void CalcPerezParam( double, double, double, in extern void CalcSkyPatchLumin( float *parr ); extern void ComputeSky( float *parr ); + +extern double solar_sunset(int month, int day); +extern double solar_sunrise(int month, int day); + /* Degrees into radians */ #define DegToRad(deg) ((deg)*(PI/180.)) /* Radiuans into degrees */ #define RadToDeg(rad) ((rad)*(180./PI)) - /* Perez sky model coefficients */ /* Reference: Perez, R., R. Seals, and J. Michalsky, 1993. "All- */ @@ -255,15 +255,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? */ @@ -283,17 +275,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, int goodsun, FILE *fp, FILE *mfp); extern void AddDirect(float *parr); @@ -315,21 +310,27 @@ getfmtname(int fmt) int main(int argc, char *argv[]) { - char buf[256]; + EPWheader *epw = NULL; /* EPW/WEA input file */ + EPWrecord erec; /* current EPW/WEA input record */ + float dpthist[2]; /* previous dew point temps */ + double dir, dif; 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 ntsteps = 0; /* number of rows in matrix */ - int step_alloc = 0; + 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 */ int mtx_offset; int i, j; + double timeinterval = 0; progname = argv[0]; /* get options */ @@ -375,25 +376,44 @@ 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 '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; + /* 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 '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); } 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 'u': /* solar hours only */ + sun_hours_only = 1; + break; case 'r': /* rotate distribution */ if (argv[i][2] && argv[i][2] != 'z') goto userr; @@ -404,24 +424,28 @@ 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; break; + case 'A': /* compute average sky */ + avgSky = 1; + break; + case 'i': + timeinterval = atof(argv[++i]); + break; default: goto userr; } - if (i < argc-1) + if ((i < argc-1) | (i > argc)) 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]); + epw = EPWopen(argv[i]); + if (epw == NULL) 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", @@ -431,28 +455,21 @@ main(int argc, char *argv[]) progname); } /* read weather tape header */ - if (scanf("place %[^\r\n] ", buf) != 1) - goto fmterr; - if (scanf("latitude %lf\n", &s_latitude) != 1) - goto fmterr; - if (scanf("longitude %lf\n", &s_longitude) != 1) - goto fmterr; - if (scanf("time_zone %lf\n", &s_meridian) != 1) - goto fmterr; - if (scanf("site_elevation %lf\n", &elevation) != 1) - goto fmterr; - if (scanf("weather_data_file_units %d\n", &input) != 1) - goto fmterr; - switch (input) { /* translate units */ - case 1: + s_latitude = epw->loc.latitude; + s_longitude = -epw->loc.longitude; + s_meridian = -15.*epw->loc.timezone; + elevation = epw->loc.elevation; + switch (epw->isWEA) { /* translate units */ + case WEAnot: + case WEAradnorm: input = 1; /* radiometric quantities */ dir_is_horiz = 0; /* direct is perpendicular meas. */ break; - case 2: + case WEAradhoriz: input = 1; /* radiometric quantities */ dir_is_horiz = 1; /* solar measured horizontally */ break; - case 3: + case WEAphotnorm: input = 2; /* photometric quantities */ dir_is_horiz = 0; /* direct is perpendicular meas. */ break; @@ -461,11 +478,16 @@ main(int argc, char *argv[]) } rh_init(); /* initialize sky patches */ if (verbose) { - fprintf(stderr, "%s: location '%s'\n", progname, buf); + fprintf(stderr, "%s: location '%s, %s'\n", progname, + epw->loc.city, epw->loc.country); 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 %s\n", + progname, sunsfp==stdout ? "stdout" : "file"); if (rotation != 0) fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation); @@ -474,30 +496,92 @@ 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); + dpthist[0] = -100; /* 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++; - if (ntsteps > step_alloc) { - step_alloc += (step_alloc>>1) + ntsteps + 7; - mtx_data = resize_dmatrix(mtx_data, step_alloc, nskypatch); - } - if (dif <= 1e-4) { - 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); + while ((j = EPWread(epw, &erec)) > 0) { + const int mo = erec.date.month+1; + const int da = erec.date.day; + const double hr = erec.date.hour; + double sda, sta, st; + int sun_in_sky; + /* 3-hour dew point temp */ + if (EPWisset(&erec,dptemp)) { + if (dpthist[0] < -99) + dpthist[0] = dpthist[1] = erec.dptemp; + dew_point = (1./3.)*(dpthist[0] + dpthist[1] + erec.dptemp); + dpthist[0] = dpthist[1]; dpthist[1] = erec.dptemp; + } else + dpthist[0] = -100; /* compute solar position */ - julian_date = jdate(mo, da); + 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); - azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); + st = hr + sta; + + if (timeinterval > 0) { + if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120) + st = (st + timeinterval/120 + solar_sunrise(mo, da))/2; + else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120) + st = (st - timeinterval/120 + solar_sunset(mo, da))/2; + } + altitude = salt(sda, st); + sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG/2.)); + if (sun_hours_only & !sun_in_sky) + continue; /* skipping nighttime points */ + azimuth = sazi(sda, st) + PI - DegToRad(rotation); + + switch (epw->isWEA) { /* translate units */ + case WEAnot: + case WEAradnorm: + if (!EPWisset(&erec,dirirrad) | + !EPWisset(&erec,horizdiffirrad)) { + fprintf(stderr, "%s: missing required irradiances at line %d\n", + progname, epw->lino); + exit(1); + } + dir = erec.dirirrad; + dif = erec.horizdiffirrad; + break; + case WEAradhoriz: + dir = erec.globhorizirrad - erec.horizdiffirrad; + dif = erec.horizdiffirrad; + break; + case WEAphotnorm: + dir = erec.dirillum; + dif = erec.diffillum; + break; + } + 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 (dir+dif <= 1e-4) { /* effectively nighttime? */ + if (!avgSky | !mtx_offset) + 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; + } + 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; @@ -508,18 +592,39 @@ 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 (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); + /* note whether leap-day was given */ } - /* check for junk at end */ - while ((i = fgetc(stdin)) != EOF) - if (!isspace(i)) { - fprintf(stderr, "%s: warning - unexpected data past EOT: ", - progname); - buf[0] = i; buf[1] = '\0'; - fgets(buf+1, sizeof(buf)-1, stdin); - fputs(buf, stderr); fputc('\n', stderr); - break; - } + if (j != EOF) { + fprintf(stderr, "%s: error on input\n", progname); + exit(1); + } + EPWclose(epw); epw = NULL; + 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; /* write out matrix */ if (outfmt != 'a') SET_FILE_BINARY(stdout); @@ -528,15 +633,17 @@ main(int argc, char *argv[]) #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", ntsteps); + printf("NCOLS=%d\n", nstored); printf("NCOMP=3\n"); + if ((outfmt == 'f') | (outfmt == 'd')) + fputendian(stdout); fputformat((char *)getfmtname(outfmt), stdout); putchar('\n'); } @@ -545,24 +652,24 @@ main(int argc, char *argv[]) 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++) { + 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]; @@ -575,23 +682,25 @@ 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][-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][-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: - 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); exit(1); } + /* Return maximum of two doubles */ double dmax( double a, double b ) { return (a > b) ? a : b; } @@ -603,7 +712,7 @@ ComputeSky(float *parr) int index; /* Category index */ double norm_diff_illum; /* Normalized diffuse illuimnance */ int i; - + /* Calculate atmospheric precipitable water content */ apwc = CalcPrecipWater(dew_point); @@ -637,10 +746,14 @@ ComputeSky(float *parr) /* Limit sky clearness */ if (sky_clearness > 11.9) sky_clearness = 11.9; + else if (sky_clearness < 1.0) + sky_clearness = 1.0; /* Limit sky brightness */ if (sky_brightness < 0.01) sky_brightness = 0.01; + else if (sky_brightness > 0.6) + sky_brightness = 0.6; /* Calculate illuminance */ index = GetCategoryIndex(); @@ -657,11 +770,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) @@ -669,6 +777,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); @@ -694,6 +806,25 @@ ComputeSky(float *parr) } } + +double +solar_sunset(int month, int day) +{ + float W; + W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day)))); + return(12 + (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15)); +} + + +double +solar_sunrise(int month, int day) +{ + float W; + W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day)))); + return(12 - (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15)); +} + + /* Add in solar direct to nearest sky patches (GW) */ void AddDirect(float *parr) @@ -746,6 +877,26 @@ AddDirect(float *parr) } } +/* Output a sun to indicated file if appropriate for this time step */ +void +OutputSun(int id, int goodsun, FILE *fp, FILE *mfp) +{ + double srad; + FVECT sv; + + srad = DegToRad(SUN_ANG_DEG/2.); + 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) */ int rh_init(void) @@ -943,6 +1094,16 @@ int CalcSkyParamFromIllum() if (sky_brightness < 0.01) sky_brightness = 0.01; + if (sky_clearness < 1.0000) + { + sky_clearness = 1.0000; + } + + if (sky_brightness > 0.6) + { + sky_brightness = 0.6; + } + while (((fabs(diff_irrad - test1) > 10.0) || (fabs(dir_irrad - test2) > 10.0)) && !(counter == 5)) { @@ -953,7 +1114,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(); @@ -966,6 +1129,16 @@ int CalcSkyParamFromIllum() /* Limit sky brightness */ if (sky_brightness < 0.01) sky_brightness = 0.01; + + if (sky_clearness < 1.0000) + { + sky_clearness = 1.0000; + } + + if (sky_brightness > 0.6) + { + sky_brightness = 0.6; + } } return GetCategoryIndex();