--- ray/src/gen/gendaymtx.c 2013/01/18 19:56:03 2.2 +++ ray/src/gen/gendaymtx.c 2013/10/04 04:31:18 2.13 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: gendaymtx.c,v 2.2 2013/01/18 19:56:03 greg Exp $"; +static const char RCSid[] = "$Id: gendaymtx.c,v 2.13 2013/10/04 04:31:18 greg Exp $"; #endif /* * gendaymtx.c @@ -108,6 +108,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 +209,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 */ @@ -246,6 +247,10 @@ static const ModelCoeff DirectLumEff[8] = { 101.18, 1.58, -1.10, -8.29 } }; +#ifndef NSUNPATCH +#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); @@ -256,7 +261,8 @@ extern double s_latitude; extern double s_longitude; extern double s_meridian; -double grefl = 0.2; /* diffuse ground reflectance */ +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? */ @@ -264,10 +270,10 @@ int outfmt = 'a'; /* output format */ int rhsubdiv = 1; /* Reinhart sky subdivisions */ -float skycolor[3] = {.96, 1.004, 1.118}; /* sky coloration */ +COLOR skycolor = {.96, 1.004, 1.118}; /* sky coloration */ +COLOR suncolor = {1., 1., 1.}; /* sun color */ +COLOR grefl = {.2, .2, .2}; /* ground reflectance */ -int do_sun = 1; /* output direct solar contribution? */ - int nskypatch; /* number of Reinhart patches */ float *rh_palt; /* sky patch altitudes (radians) */ float *rh_pazi; /* sky patch azimuths (radians) */ @@ -290,6 +296,7 @@ int main(int argc, char *argv[]) { char buf[256]; + 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 */ @@ -305,13 +312,15 @@ main(int argc, char *argv[]) /* get options */ for (i = 1; i < argc && argv[i][0] == '-'; i++) switch (argv[i][1]) { - case 'g': - grefl = atof(argv[++i]); + case 'g': /* ground reflectance */ + grefl[0] = atof(argv[++i]); + grefl[1] = atof(argv[++i]); + grefl[2] = atof(argv[++i]); break; - case 'v': + case 'v': /* verbose progress reports */ verbose++; break; - case 'o': + case 'o': /* output format */ switch (argv[i][2]) { case 'f': case 'd': @@ -322,23 +331,47 @@ main(int argc, char *argv[]) goto userr; } break; - case 'm': + 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': + case 'c': /* sky color */ skycolor[0] = atof(argv[++i]); skycolor[1] = atof(argv[++i]); skycolor[2] = atof(argv[++i]); break; - case 'd': - do_sun = 1; + case 'd': /* solar (direct) only */ skycolor[0] = skycolor[1] = skycolor[2] = 0; + if (suncolor[1] <= 1e-4) + suncolor[0] = suncolor[1] = suncolor[2] = 1; break; - case 's': - do_sun = 0; + case 's': /* sky only (no direct) */ + suncolor[0] = suncolor[1] = suncolor[2] = 0; if (skycolor[1] <= 1e-4) 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 = 6.797e-05; + break; default: goto userr; } @@ -393,6 +426,9 @@ main(int argc, char *argv[]) progname, s_latitude, s_longitude); fprintf(stderr, "%s: %d sky patches per time step\n", progname, nskypatch); + if (rotation != 0) + fprintf(stderr, "%s: rotating output %.0f degrees\n", + progname, rotation); } /* convert quantities to radians */ s_latitude = DegToRad(s_latitude); @@ -416,7 +452,7 @@ main(int argc, char *argv[]) sda = sdec(julian_date); sta = stadj(julian_date); altitude = salt(sda, hr+sta); - azimuth = sazi(sda, hr+sta) + PI; + azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); /* convert measured values */ if (dir_is_horiz && altitude > 0.) dir /= sin(altitude); @@ -429,8 +465,7 @@ main(int argc, char *argv[]) } /* compute sky patch values */ ComputeSky(mtx_data+mtx_offset); - if (do_sun) - AddDirect(mtx_data+mtx_offset); + AddDirect(mtx_data+mtx_offset); } /* check for junk at end */ while ((i = fgetc(stdin)) != EOF) @@ -455,7 +490,7 @@ main(int argc, char *argv[]) switch (outfmt) { case 'a': for (j = 0; j < ntsteps; j++) { - printf("%.3e %.3e %.3e\n", mtx_data[mtx_offset], + printf("%.3g %.3g %.3g\n", mtx_data[mtx_offset], mtx_data[mtx_offset+1], mtx_data[mtx_offset+2]); mtx_offset += 3*nskypatch; @@ -490,7 +525,7 @@ main(int argc, char *argv[]) fprintf(stderr, "%s: done.\n", progname); exit(0); userr: - fprintf(stderr, "Usage: %s [-v][-d|-s][-m N][-g refl][-c r g b][-o{f|d}] [tape.wea]\n", + fprintf(stderr, "Usage: %s [-v][-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: @@ -511,19 +546,20 @@ ComputeSky(float *parr) { int index; /* Category index */ double norm_diff_illum; /* Normalized diffuse illuimnance */ - double zlumin; /* Zenith luminance */ int i; /* Calculate atmospheric precipitable water content */ apwc = CalcPrecipWater(dew_point); - /* Limit solar altitude to keep circumsolar off zenith */ - if (altitude > DegToRad(87.0)) - altitude = DegToRad(87.0); + /* Calculate sun zenith angle (don't let it dip below horizon) */ + /* Also limit minimum angle to keep circumsolar off zenith */ + if (altitude <= 0.0) + sun_zenith = DegToRad(90.0); + else if (altitude >= DegToRad(87.0)) + sun_zenith = DegToRad(3.0); + else + sun_zenith = DegToRad(90.0) - altitude; - /* Calculate sun zenith angle */ - sun_zenith = DegToRad(90.0) - altitude; - /* Compute the inputs for the calculation of the sky distribution */ if (input == 0) /* XXX never used */ @@ -542,6 +578,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); @@ -553,15 +597,21 @@ 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; } /* Compute ground radiance (include solar contribution if any) */ - parr[0] = diff_illum * (1./PI/WHTEFFICACY); + parr[0] = diff_illum; if (altitude > 0) - parr[0] += dir_illum * sin(altitude) * (1./PI/WHTEFFICACY); - parr[2] = parr[1] = parr[0]; + parr[0] += dir_illum * sin(altitude); + parr[2] = parr[1] = parr[0] *= (1./PI/WHTEFFICACY); + multcolor(parr, grefl); /* Calculate Perez sky model parameters */ CalcPerezParam(sun_zenith, sky_clearness, sky_brightness, index); @@ -572,18 +622,18 @@ 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; - /* Calculate relative zenith luminance */ - zlumin = CalcRelLuminance(sun_zenith, 0.0); - - /* Calculate absolute zenith illuminance */ - zlumin *= norm_diff_illum; - /* Apply to sky patches to get absolute radiance values */ for (i = 1; i < nskypatch; i++) { - scalecolor(parr+3*i, zlumin*(1./WHTEFFICACY)); + scalecolor(parr+3*i, norm_diff_illum*(1./WHTEFFICACY)); multcolor(parr+3*i, skycolor); } } @@ -593,15 +643,19 @@ void AddDirect(float *parr) { FVECT svec; - double near_dprod[4]; - int near_patch[4]; - double wta[4], wtot; + double near_dprod[NSUNPATCH]; + int near_patch[NSUNPATCH]; + double wta[NSUNPATCH], wtot; int i, j, p; - if (!do_sun || dir_illum < 1e-4) + if (dir_illum <= 1e-4 || bright(suncolor) <= 1e-4) return; - /* identify 4 closest patches */ - for (i = 4; 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++) { @@ -609,9 +663,9 @@ AddDirect(float *parr) double dprod; rh_vector(pvec, p); dprod = DOT(pvec, svec); - for (i = 0; i < 4; i++) + for (i = 0; i < nsuns; i++) if (dprod > near_dprod[i]) { - for (j = 4; --j > i; ) { + for (j = nsuns; --j > i; ) { near_dprod[j] = near_dprod[j-1]; near_patch[j] = near_patch[j-1]; } @@ -621,16 +675,18 @@ AddDirect(float *parr) } } wtot = 0; /* weight by proximity */ - for (i = 4; i--; ) + for (i = nsuns; i--; ) wtot += wta[i] = 1./(1.002 - near_dprod[i]); /* add to nearest patch radiances */ - for (i = 4; 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]]); - *pdest++ += val_add; - *pdest++ += val_add; - *pdest++ += val_add; + 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]; } } @@ -665,7 +721,8 @@ rh_init(void) for (i = 0; i < NROW*rhsubdiv; i++) { const float ralt = alpha*(i + .5); const int ninrow = tnaz[i/rhsubdiv]*rhsubdiv; - const float dom = (sin(alpha*(i+1)) - sin(alpha*i))/ninrow; + const float dom = 2.*PI*(sin(alpha*(i+1)) - sin(alpha*i)) / + (double)ninrow; for (j = 0; j < ninrow; j++) { rh_palt[p] = ralt; rh_pazi[p] = 2.*PI * j / (double)ninrow; @@ -770,7 +827,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); @@ -801,7 +858,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 */ @@ -827,7 +884,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) || @@ -851,7 +908,7 @@ int CalcSkyParamFromIllum() sky_clearness = 12.0; /* Limit sky brightness */ - if (sky_brightness < 0.05) + if (sky_brightness < 0.01) sky_brightness = 0.01; } @@ -937,9 +994,9 @@ double CalcRelHorzIllum( float *parr ) double rh_illum = 0.0; /* Relative horizontal illuminance */ for (i = 1; i < nskypatch; i++) - rh_illum += parr[3*i+1] * rh_cos(i); + rh_illum += parr[3*i+1] * rh_cos(i) * rh_dom[i]; - return rh_illum * (2.0 * PI / (nskypatch-1)); + return rh_illum; } /* Calculate earth orbit eccentricity correction factor */