89 |
|
#include "color.h" |
90 |
|
#include "sun.h" |
91 |
|
|
92 |
< |
char *progname; /* Program name */ |
93 |
< |
char errmsg[128]; /* Error message buffer */ |
92 |
> |
char *progname; /* Program name */ |
93 |
|
const double DC_SolarConstantE = 1367.0; /* Solar constant W/m^2 */ |
94 |
|
const double DC_SolarConstantL = 127.5; /* Solar constant klux */ |
95 |
|
|
134 |
|
/* Radiuans into degrees */ |
135 |
|
#define RadToDeg(rad) ((rad)*(180./PI)) |
136 |
|
|
138 |
– |
|
137 |
|
/* Perez sky model coefficients */ |
138 |
|
|
139 |
|
/* Reference: Perez, R., R. Seals, and J. Michalsky, 1993. "All- */ |
310 |
|
double rotation = 0; /* site rotation (degrees) */ |
311 |
|
double elevation; /* site elevation (meters) */ |
312 |
|
int leap_day = 0; /* add leap day? */ |
313 |
+ |
int sun_hours_only = 0; /* only output sun hours? */ |
314 |
|
int dir_is_horiz; /* direct is meas. on horizontal? */ |
315 |
|
FILE *sunsfp = NULL; /* output file for individual suns */ |
316 |
|
FILE *modsfp = NULL; /* modifier output file */ |
405 |
|
case 's': /* sky only (no direct) */ |
406 |
|
suncolor[0] = suncolor[1] = suncolor[2] = 0; |
407 |
|
break; |
408 |
+ |
case 'u': /* solar hours only */ |
409 |
+ |
sun_hours_only = 1; |
410 |
+ |
break; |
411 |
|
case 'r': /* rotate distribution */ |
412 |
|
if (argv[i][2] && argv[i][2] != 'z') |
413 |
|
goto userr; |
485 |
|
fprintf(stderr, "%s: %d sky patches\n", |
486 |
|
progname, nskypatch); |
487 |
|
if (sunsfp) |
488 |
< |
fprintf(stderr, "%s: outputting suns to file\n", |
489 |
< |
progname); |
488 |
> |
fprintf(stderr, "%s: outputting suns to %s\n", |
489 |
> |
progname, sunsfp==stdout ? "stdout" : "file"); |
490 |
|
if (rotation != 0) |
491 |
|
fprintf(stderr, "%s: rotating output %.0f degrees\n", |
492 |
|
progname, rotation); |
500 |
|
/* process each time step in tape */ |
501 |
|
while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) { |
502 |
|
double sda, sta; |
503 |
< |
|
502 |
< |
mtx_offset = 3*nskypatch*nstored; |
503 |
< |
nstored += !avgSky | !nstored; |
504 |
< |
/* make space for next row */ |
505 |
< |
if (nstored > tstorage) { |
506 |
< |
tstorage += (tstorage>>1) + nstored + 7; |
507 |
< |
mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch); |
508 |
< |
} |
509 |
< |
ntsteps++; /* keep count of time steps */ |
503 |
> |
int sun_in_sky; |
504 |
|
/* compute solar position */ |
505 |
|
if ((mo == 2) & (da == 29)) { |
506 |
|
julian_date = 60; |
510 |
|
sda = sdec(julian_date); |
511 |
|
sta = stadj(julian_date); |
512 |
|
altitude = salt(sda, hr+sta); |
513 |
+ |
sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG/2.)); |
514 |
+ |
if (sun_hours_only && !sun_in_sky) |
515 |
+ |
continue; /* skipping nighttime points */ |
516 |
|
azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); |
517 |
|
|
518 |
+ |
mtx_offset = 3*nskypatch*nstored; |
519 |
+ |
nstored += !avgSky | !nstored; |
520 |
+ |
/* make space for next row */ |
521 |
+ |
if (nstored > tstorage) { |
522 |
+ |
tstorage += (tstorage>>1) + nstored + 7; |
523 |
+ |
mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch); |
524 |
+ |
} |
525 |
+ |
ntsteps++; /* keep count of time steps */ |
526 |
+ |
|
527 |
|
if (dir+dif <= 1e-4) { /* effectively nighttime? */ |
528 |
|
if (!avgSky | !mtx_offset) |
529 |
|
memset(mtx_data+mtx_offset, 0, |
530 |
|
sizeof(float)*3*nskypatch); |
531 |
< |
if (sunsfp) /* output black sun */ |
531 |
> |
/* output black sun? */ |
532 |
> |
if (sunsfp && sun_in_sky) |
533 |
|
OutputSun(solar_minute(julian_date,hr), 0, |
534 |
|
sunsfp, modsfp); |
535 |
|
continue; |
536 |
|
} |
537 |
+ |
if (!sun_in_sky && dir > (input==1 ? 20. : 20.*WHTEFFICACY)) |
538 |
+ |
fprintf(stderr, |
539 |
+ |
"%s: warning - unusually bright at %.1f on %d-%d\n", |
540 |
+ |
progname, hr, mo, da); |
541 |
|
/* convert measured values */ |
542 |
< |
if (dir_is_horiz && altitude > 0.) |
542 |
> |
if (dir_is_horiz && altitude > FTINY) |
543 |
|
dir /= sin(altitude); |
544 |
|
if (input == 1) { |
545 |
|
dir_irrad = dir; |
550 |
|
} |
551 |
|
/* compute sky patch values */ |
552 |
|
ComputeSky(mtx_data+mtx_offset); |
553 |
< |
|
554 |
< |
if (sunsfp) /* output sun if requested */ |
553 |
> |
/* output sun if requested */ |
554 |
> |
if (sunsfp && sun_in_sky) |
555 |
|
OutputSun(solar_minute(julian_date,hr), 1, |
556 |
|
sunsfp, modsfp); |
557 |
|
|
653 |
|
fprintf(stderr, "%s: done.\n", progname); |
654 |
|
exit(0); |
655 |
|
userr: |
656 |
< |
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", |
656 |
> |
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", |
657 |
|
progname); |
658 |
|
exit(1); |
659 |
|
fmterr: |
825 |
|
FVECT sv; |
826 |
|
|
827 |
|
srad = DegToRad(SUN_ANG_DEG/2.); |
817 |
– |
|
818 |
– |
if (altitude < -srad) /* well below horizon? */ |
819 |
– |
return; |
820 |
– |
|
828 |
|
srad = goodsun ? dir_illum/(WHTEFFICACY * PI*srad*srad) : 0; |
829 |
|
vector(sv, altitude, azimuth); |
830 |
|
fprintf(fp, "\nvoid light solar%d\n0\n0\n", id); |