| 272 |
|
float *rh_pazi; /* sky patch azimuths (radians) */ |
| 273 |
|
float *rh_dom; /* sky patch solid angle (sr) */ |
| 274 |
|
|
| 275 |
< |
#define vector(v,alt,azi) ( (v)[1] = tcos(alt), \ |
| 276 |
< |
(v)[0] = (v)[1]*tsin(azi), \ |
| 277 |
< |
(v)[1] *= tcos(azi), \ |
| 278 |
< |
(v)[2] = tsin(alt) ) |
| 275 |
> |
#define vector(v,alt,azi) ( (v)[1] = cos(alt), \ |
| 276 |
> |
(v)[0] = (v)[1]*sin(azi), \ |
| 277 |
> |
(v)[1] *= cos(azi), \ |
| 278 |
> |
(v)[2] = sin(alt) ) |
| 279 |
|
|
| 280 |
|
#define rh_vector(v,i) vector(v,rh_palt[i],rh_pazi[i]) |
| 281 |
|
|
| 282 |
|
#define rh_cos(i) tsin(rh_palt[i]) |
| 283 |
|
|
| 284 |
+ |
#define solar_minute(jd,hr) ((24*60)*((jd)-1)+(int)((hr)*60.+.5)) |
| 285 |
+ |
|
| 286 |
|
extern int rh_init(void); |
| 287 |
|
extern float * resize_dmatrix(float *mtx_data, int nsteps, int npatch); |
| 288 |
< |
extern void OutputSun(int id, FILE *fp); |
| 288 |
> |
extern void OutputSun(int id, int goodsun, FILE *fp, FILE *mfp); |
| 289 |
|
extern void AddDirect(float *parr); |
| 290 |
|
|
| 291 |
|
|
| 311 |
|
int doheader = 1; /* output header? */ |
| 312 |
|
double rotation = 0; /* site rotation (degrees) */ |
| 313 |
|
double elevation; /* site elevation (meters) */ |
| 314 |
+ |
int leap_day = 0; /* add leap day? */ |
| 315 |
|
int dir_is_horiz; /* direct is meas. on horizontal? */ |
| 316 |
|
FILE *sunsfp = NULL; /* output file for individual suns */ |
| 317 |
+ |
FILE *modsfp = NULL; /* modifier output file */ |
| 318 |
|
float *mtx_data = NULL; /* our matrix data */ |
| 319 |
|
int avgSky = 0; /* compute average sky r.t. matrix? */ |
| 320 |
|
int ntsteps = 0; /* number of time steps */ |
| 375 |
|
skycolor[1] = atof(argv[++i]); |
| 376 |
|
skycolor[2] = atof(argv[++i]); |
| 377 |
|
break; |
| 378 |
+ |
case 'D': /* output suns to file */ |
| 379 |
+ |
if (strcmp(argv[++i], "-")) { |
| 380 |
+ |
sunsfp = fopen(argv[i], "w"); |
| 381 |
+ |
if (sunsfp == NULL) { |
| 382 |
+ |
fprintf(stderr, |
| 383 |
+ |
"%s: cannot open '%s' for output\n", |
| 384 |
+ |
progname, argv[i]); |
| 385 |
+ |
exit(1); |
| 386 |
+ |
} |
| 387 |
+ |
break; /* still may output matrix */ |
| 388 |
+ |
} |
| 389 |
+ |
sunsfp = stdout; /* sending to stdout, so... */ |
| 390 |
+ |
/* fall through */ |
| 391 |
|
case 'n': /* no matrix output */ |
| 392 |
|
avgSky = -1; |
| 393 |
|
rhsubdiv = 1; |
| 396 |
|
skycolor[0] = skycolor[1] = skycolor[2] = 0; |
| 397 |
|
grefl[0] = grefl[1] = grefl[2] = 0; |
| 398 |
|
break; |
| 399 |
< |
case 'D': /* output suns to file */ |
| 400 |
< |
sunsfp = fopen(argv[++i], "w"); |
| 384 |
< |
if (!sunsfp) { |
| 399 |
> |
case 'M': /* send sun modifiers to file */ |
| 400 |
> |
if ((modsfp = fopen(argv[++i], "w")) == NULL) { |
| 401 |
|
fprintf(stderr, "%s: cannot open '%s' for output\n", |
| 402 |
|
progname, argv[i]); |
| 403 |
|
exit(1); |
| 434 |
|
progname, argv[i]); |
| 435 |
|
exit(1); |
| 436 |
|
} |
| 437 |
+ |
if ((modsfp != NULL) & (sunsfp == NULL)) |
| 438 |
+ |
fprintf(stderr, "%s: warning -M output will be empty without -D\n", |
| 439 |
+ |
progname); |
| 440 |
|
if (verbose) { |
| 441 |
|
if (i == argc-1) |
| 442 |
|
fprintf(stderr, "%s: reading weather tape '%s'\n", |
| 507 |
|
mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch); |
| 508 |
|
} |
| 509 |
|
ntsteps++; /* keep count of time steps */ |
| 491 |
– |
if (dif <= 1e-4) { |
| 492 |
– |
if (!avgSky | !mtx_offset) |
| 493 |
– |
memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch); |
| 494 |
– |
continue; |
| 495 |
– |
} |
| 510 |
|
/* compute solar position */ |
| 511 |
< |
julian_date = jdate(mo, da); |
| 511 |
> |
if ((mo == 2) & (da == 29)) { |
| 512 |
> |
julian_date = 60; |
| 513 |
> |
leap_day = 1; |
| 514 |
> |
} else |
| 515 |
> |
julian_date = jdate(mo, da) + leap_day; |
| 516 |
|
sda = sdec(julian_date); |
| 517 |
|
sta = stadj(julian_date); |
| 518 |
|
altitude = salt(sda, hr+sta); |
| 519 |
|
azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); |
| 520 |
+ |
|
| 521 |
+ |
if (dir+dif <= 1e-4) { /* effectively nighttime? */ |
| 522 |
+ |
if (!avgSky | !mtx_offset) |
| 523 |
+ |
memset(mtx_data+mtx_offset, 0, |
| 524 |
+ |
sizeof(float)*3*nskypatch); |
| 525 |
+ |
if (sunsfp) /* output black sun */ |
| 526 |
+ |
OutputSun(solar_minute(julian_date,hr), 0, |
| 527 |
+ |
sunsfp, modsfp); |
| 528 |
+ |
continue; |
| 529 |
+ |
} |
| 530 |
|
/* convert measured values */ |
| 531 |
|
if (dir_is_horiz && altitude > 0.) |
| 532 |
|
dir /= sin(altitude); |
| 540 |
|
/* compute sky patch values */ |
| 541 |
|
ComputeSky(mtx_data+mtx_offset); |
| 542 |
|
|
| 543 |
< |
if (sunsfp) /* output sun if indicated */ |
| 544 |
< |
OutputSun(ntsteps, sunsfp); |
| 543 |
> |
if (sunsfp) /* output sun if requested */ |
| 544 |
> |
OutputSun(solar_minute(julian_date,hr), 1, |
| 545 |
> |
sunsfp, modsfp); |
| 546 |
|
|
| 547 |
|
if (avgSky < 0) /* no matrix? */ |
| 548 |
|
continue; |
| 555 |
|
if (verbose && mo != last_monthly) |
| 556 |
|
fprintf(stderr, "%s: stepping through month %d...\n", |
| 557 |
|
progname, last_monthly=mo); |
| 558 |
+ |
/* note whether leap-day was given */ |
| 559 |
|
} |
| 560 |
|
if (!ntsteps) { |
| 561 |
|
fprintf(stderr, "%s: no valid time steps on input\n", progname); |
| 808 |
|
|
| 809 |
|
/* Output a sun to indicated file if appropriate for this time step */ |
| 810 |
|
void |
| 811 |
< |
OutputSun(int id, FILE *fp) |
| 811 |
> |
OutputSun(int id, int goodsun, FILE *fp, FILE *mfp) |
| 812 |
|
{ |
| 813 |
|
double srad; |
| 814 |
|
FVECT sv; |
| 815 |
|
|
| 816 |
< |
if ((altitude <= 0) | (dir_illum <= 1e-4)) |
| 816 |
> |
srad = DegToRad(SUN_ANG_DEG/2.); |
| 817 |
> |
|
| 818 |
> |
if (altitude < -srad) /* well below horizon? */ |
| 819 |
|
return; |
| 820 |
|
|
| 821 |
< |
srad = DegToRad(SUN_ANG_DEG/2.); |
| 790 |
< |
srad = dir_illum/(WHTEFFICACY * PI*srad*srad); |
| 821 |
> |
srad = goodsun ? dir_illum/(WHTEFFICACY * PI*srad*srad) : 0; |
| 822 |
|
vector(sv, altitude, azimuth); |
| 823 |
|
fprintf(fp, "\nvoid light solar%d\n0\n0\n", id); |
| 824 |
|
fprintf(fp, "3 %.3e %.3e %.3e\n", srad*suncolor[0], |
| 825 |
|
srad*suncolor[1], srad*suncolor[2]); |
| 826 |
|
fprintf(fp, "\nsolar%d source sun%d\n0\n0\n", id, id); |
| 827 |
|
fprintf(fp, "4 %.6f %.6f %.6f %.4f\n", sv[0], sv[1], sv[2], SUN_ANG_DEG); |
| 828 |
+ |
|
| 829 |
+ |
if (mfp != NULL) /* saving modifier IDs? */ |
| 830 |
+ |
fprintf(mfp, "solar%d\n", id); |
| 831 |
|
} |
| 832 |
|
|
| 833 |
|
/* Initialize Reinhart sky patch positions (GW) */ |