--- ray/src/gen/gendaymtx.c 2013/01/18 01:12:59 2.1 +++ ray/src/gen/gendaymtx.c 2013/01/20 01:16:15 2.4 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: gendaymtx.c,v 2.1 2013/01/18 01:12:59 greg Exp $"; +static const char RCSid[] = "$Id: gendaymtx.c,v 2.4 2013/01/20 01:16:15 greg Exp $"; #endif /* * gendaymtx.c @@ -246,6 +246,10 @@ static const ModelCoeff DirectLumEff[8] = { 101.18, 1.58, -1.10, -8.29 } }; +#ifndef NSUNPATCH +#define NSUNPATCH 4 /* # patches to spread sun into */ +#endif + extern int jdate(int month, int day); extern double stadj(int jd); extern double sdec(int jd); @@ -256,18 +260,16 @@ extern double s_latitude; extern double s_longitude; extern double s_meridian; -double grefl = 0.2; /* diffuse ground reflectance */ - int verbose = 0; /* progress reports to stderr? */ 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) */ @@ -305,13 +307,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,20 +326,21 @@ main(int argc, char *argv[]) goto userr; } break; - case 'm': + 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; @@ -358,7 +363,7 @@ main(int argc, char *argv[]) progname); } /* read weather tape header */ - if (scanf("place %[^\n]\n", buf) != 1) + if (scanf("place %[^\r\n] ", buf) != 1) goto fmterr; if (scanf("latitude %lf\n", &s_latitude) != 1) goto fmterr; @@ -394,6 +399,10 @@ main(int argc, char *argv[]) fprintf(stderr, "%s: %d sky patches per time step\n", progname, nskypatch); } + /* convert quantities to radians */ + s_latitude = DegToRad(s_latitude); + s_longitude = DegToRad(s_longitude); + s_meridian = DegToRad(s_meridian); /* process each time step in tape */ while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) { double sda, sta; @@ -412,7 +421,7 @@ main(int argc, char *argv[]) sda = sdec(julian_date); sta = stadj(julian_date); altitude = salt(sda, hr+sta); - azimuth = sazi(sda, hr+sta); + azimuth = sazi(sda, hr+sta) + PI; /* convert measured values */ if (dir_is_horiz && altitude > 0.) dir /= sin(altitude); @@ -425,8 +434,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) @@ -451,12 +459,13 @@ 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; } - fputc('\n', stdout); + if (ntsteps > 1) + fputc('\n', stdout); break; case 'f': for (j = 0; j < ntsteps; j++) { @@ -508,11 +517,6 @@ ComputeSky(float *parr) double norm_diff_illum; /* Normalized diffuse illuimnance */ double zlumin; /* Zenith luminance */ int i; - - if (bright(skycolor) <= 1e-4) { /* 0 sky component? */ - memset(parr, 0, sizeof(float)*3*nskypatch); - return; - } /* Calculate atmospheric precipitable water content */ apwc = CalcPrecipWater(dew_point); @@ -553,11 +557,16 @@ ComputeSky(float *parr) index = CalcSkyParamFromIllum(); } + 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); @@ -589,15 +598,15 @@ 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 NSUNPATCH closest patches */ + for (i = NSUNPATCH; i--; ) near_dprod[i] = -1.; vector(svec, altitude, azimuth); for (p = 1; p < nskypatch; p++) { @@ -605,9 +614,9 @@ AddDirect(float *parr) double dprod; rh_vector(pvec, p); dprod = DOT(pvec, svec); - for (i = 0; i < 4; i++) + for (i = 0; i < NSUNPATCH; i++) if (dprod > near_dprod[i]) { - for (j = 4; --j > i; ) { + for (j = NSUNPATCH; --j > i; ) { near_dprod[j] = near_dprod[j-1]; near_patch[j] = near_patch[j-1]; } @@ -617,12 +626,17 @@ AddDirect(float *parr) } } wtot = 0; /* weight by proximity */ - for (i = 4; i--; ) + for (i = NSUNPATCH; i--; ) wtot += wta[i] = 1./(1.002 - near_dprod[i]); /* add to nearest patch radiances */ - for (i = 4; i--; ) - parr[near_patch[i]] += wta[i] * dir_illum / - (wtot * rh_dom[near_patch[i]]); + for (i = NSUNPATCH; 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*suncolor[0]; + *pdest++ += val_add*suncolor[1]; + *pdest++ += val_add*suncolor[2]; + } } /* Initialize Reinhart sky patch positions (GW) */ @@ -656,7 +670,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;