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) */ |