| 88 |
|
#include "rtio.h" |
| 89 |
|
#include "color.h" |
| 90 |
|
#include "sun.h" |
| 91 |
+ |
#include "loadEPW.h" |
| 92 |
|
|
| 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 |
|
|
| 309 |
|
int |
| 310 |
|
main(int argc, char *argv[]) |
| 311 |
|
{ |
| 312 |
< |
char buf[256]; |
| 312 |
> |
EPWheader *epw = NULL; /* EPW/WEA input file */ |
| 313 |
> |
EPWrecord erec; /* current EPW/WEA input record */ |
| 314 |
> |
float dpthist[2]; /* previous dew point temps */ |
| 315 |
> |
double dir, dif; |
| 316 |
|
int doheader = 1; /* output header? */ |
| 317 |
|
double rotation = 0; /* site rotation (degrees) */ |
| 318 |
|
double elevation; /* site elevation (meters) */ |
| 327 |
|
int tstorage = 0; /* number of allocated time steps */ |
| 328 |
|
int nstored = 0; /* number of time steps in matrix */ |
| 329 |
|
int last_monthly = 0; /* month of last report */ |
| 327 |
– |
int mo, da; /* month (1-12) and day (1-31) */ |
| 328 |
– |
double hr; /* hour (local standard time) */ |
| 329 |
– |
double dir, dif; /* direct and diffuse values */ |
| 330 |
|
int mtx_offset; |
| 331 |
|
int i, j; |
| 332 |
|
double timeinterval = 0; |
| 333 |
|
|
| 334 |
< |
progname = argv[0]; |
| 334 |
> |
fixargv0(argv[0]); |
| 335 |
|
/* get options */ |
| 336 |
|
for (i = 1; i < argc && argv[i][0] == '-'; i++) |
| 337 |
|
switch (argv[i][1]) { |
| 437 |
|
default: |
| 438 |
|
goto userr; |
| 439 |
|
} |
| 440 |
< |
if (i < argc-1) |
| 440 |
> |
if ((i < argc-1) | (i > argc)) |
| 441 |
|
goto userr; |
| 442 |
< |
if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) { |
| 443 |
< |
fprintf(stderr, "%s: cannot open '%s' for input\n", |
| 444 |
< |
progname, argv[i]); |
| 442 |
> |
epw = EPWopen(argv[i]); |
| 443 |
> |
if (epw == NULL) |
| 444 |
|
exit(1); |
| 446 |
– |
} |
| 445 |
|
if ((modsfp != NULL) & (sunsfp == NULL)) |
| 446 |
|
fprintf(stderr, "%s: warning -M output will be empty without -D\n", |
| 447 |
|
progname); |
| 454 |
|
progname); |
| 455 |
|
} |
| 456 |
|
/* read weather tape header */ |
| 457 |
< |
if (scanf("place %[^\r\n] ", buf) != 1) |
| 458 |
< |
goto fmterr; |
| 459 |
< |
if (scanf("latitude %lf\n", &s_latitude) != 1) |
| 460 |
< |
goto fmterr; |
| 461 |
< |
if (scanf("longitude %lf\n", &s_longitude) != 1) |
| 462 |
< |
goto fmterr; |
| 463 |
< |
if (scanf("time_zone %lf\n", &s_meridian) != 1) |
| 466 |
< |
goto fmterr; |
| 467 |
< |
if (scanf("site_elevation %lf\n", &elevation) != 1) |
| 468 |
< |
goto fmterr; |
| 469 |
< |
if (scanf("weather_data_file_units %d\n", &input) != 1) |
| 470 |
< |
goto fmterr; |
| 471 |
< |
switch (input) { /* translate units */ |
| 472 |
< |
case 1: |
| 457 |
> |
s_latitude = epw->loc.latitude; |
| 458 |
> |
s_longitude = -epw->loc.longitude; |
| 459 |
> |
s_meridian = -15.*epw->loc.timezone; |
| 460 |
> |
elevation = epw->loc.elevation; |
| 461 |
> |
switch (epw->isWEA) { /* translate units */ |
| 462 |
> |
case WEAnot: |
| 463 |
> |
case WEAradnorm: |
| 464 |
|
input = 1; /* radiometric quantities */ |
| 465 |
|
dir_is_horiz = 0; /* direct is perpendicular meas. */ |
| 466 |
|
break; |
| 467 |
< |
case 2: |
| 467 |
> |
case WEAradhoriz: |
| 468 |
|
input = 1; /* radiometric quantities */ |
| 469 |
|
dir_is_horiz = 1; /* solar measured horizontally */ |
| 470 |
|
break; |
| 471 |
< |
case 3: |
| 471 |
> |
case WEAphotnorm: |
| 472 |
|
input = 2; /* photometric quantities */ |
| 473 |
|
dir_is_horiz = 0; /* direct is perpendicular meas. */ |
| 474 |
|
break; |
| 477 |
|
} |
| 478 |
|
rh_init(); /* initialize sky patches */ |
| 479 |
|
if (verbose) { |
| 480 |
< |
fprintf(stderr, "%s: location '%s'\n", progname, buf); |
| 480 |
> |
fprintf(stderr, "%s: location '%s, %s'\n", progname, |
| 481 |
> |
epw->loc.city, epw->loc.country); |
| 482 |
|
fprintf(stderr, "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n", |
| 483 |
|
progname, s_latitude, s_longitude); |
| 484 |
|
if (avgSky >= 0) |
| 497 |
|
s_meridian = DegToRad(s_meridian); |
| 498 |
|
/* initial allocation */ |
| 499 |
|
mtx_data = resize_dmatrix(mtx_data, tstorage=2, nskypatch); |
| 500 |
+ |
dpthist[0] = -100; |
| 501 |
|
/* process each time step in tape */ |
| 502 |
< |
while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) { |
| 502 |
> |
while ((j = EPWread(epw, &erec)) > 0) { |
| 503 |
> |
const int mo = erec.date.month+1; |
| 504 |
> |
const int da = erec.date.day; |
| 505 |
> |
const double hr = erec.date.hour; |
| 506 |
|
double sda, sta, st; |
| 507 |
|
int sun_in_sky; |
| 508 |
+ |
/* 3-hour dew point temp */ |
| 509 |
+ |
if (EPWisset(&erec,dptemp)) { |
| 510 |
+ |
if (dpthist[0] < -99) |
| 511 |
+ |
dpthist[0] = dpthist[1] = erec.dptemp; |
| 512 |
+ |
dew_point = (1./3.)*(dpthist[0] + dpthist[1] + erec.dptemp); |
| 513 |
+ |
dpthist[0] = dpthist[1]; dpthist[1] = erec.dptemp; |
| 514 |
+ |
} else |
| 515 |
+ |
dpthist[0] = -100; |
| 516 |
|
/* compute solar position */ |
| 517 |
|
if ((mo == 2) & (da == 29)) { |
| 518 |
|
julian_date = 60; |
| 531 |
|
} |
| 532 |
|
altitude = salt(sda, st); |
| 533 |
|
sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG/2.)); |
| 534 |
< |
if (sun_hours_only && !sun_in_sky) |
| 534 |
> |
if (sun_hours_only & !sun_in_sky) |
| 535 |
|
continue; /* skipping nighttime points */ |
| 536 |
|
azimuth = sazi(sda, st) + PI - DegToRad(rotation); |
| 537 |
|
|
| 538 |
+ |
switch (epw->isWEA) { /* translate units */ |
| 539 |
+ |
case WEAnot: |
| 540 |
+ |
case WEAradnorm: |
| 541 |
+ |
if (!EPWisset(&erec,dirirrad) | |
| 542 |
+ |
!EPWisset(&erec,horizdiffirrad)) { |
| 543 |
+ |
fprintf(stderr, "%s: missing required irradiances at line %d\n", |
| 544 |
+ |
progname, epw->lino); |
| 545 |
+ |
exit(1); |
| 546 |
+ |
} |
| 547 |
+ |
dir = erec.dirirrad; |
| 548 |
+ |
dif = erec.horizdiffirrad; |
| 549 |
+ |
break; |
| 550 |
+ |
case WEAradhoriz: |
| 551 |
+ |
dir = erec.globhorizirrad - erec.horizdiffirrad; |
| 552 |
+ |
dif = erec.horizdiffirrad; |
| 553 |
+ |
break; |
| 554 |
+ |
case WEAphotnorm: |
| 555 |
+ |
dir = erec.dirillum; |
| 556 |
+ |
dif = erec.diffillum; |
| 557 |
+ |
break; |
| 558 |
+ |
} |
| 559 |
|
mtx_offset = 3*nskypatch*nstored; |
| 560 |
|
nstored += !avgSky | !nstored; |
| 561 |
|
/* make space for next row */ |
| 609 |
|
progname, last_monthly=mo); |
| 610 |
|
/* note whether leap-day was given */ |
| 611 |
|
} |
| 612 |
+ |
if (j != EOF) { |
| 613 |
+ |
fprintf(stderr, "%s: error on input\n", progname); |
| 614 |
+ |
exit(1); |
| 615 |
+ |
} |
| 616 |
+ |
EPWclose(epw); epw = NULL; |
| 617 |
|
if (!ntsteps) { |
| 618 |
|
fprintf(stderr, "%s: no valid time steps on input\n", progname); |
| 619 |
|
exit(1); |
| 620 |
|
} |
| 591 |
– |
/* check for junk at end */ |
| 592 |
– |
while ((i = fgetc(stdin)) != EOF) |
| 593 |
– |
if (!isspace(i)) { |
| 594 |
– |
fprintf(stderr, "%s: warning - unexpected data past EOT: ", |
| 595 |
– |
progname); |
| 596 |
– |
buf[0] = i; buf[1] = '\0'; |
| 597 |
– |
fgets(buf+1, sizeof(buf)-1, stdin); |
| 598 |
– |
fputs(buf, stderr); fputc('\n', stderr); |
| 599 |
– |
break; |
| 600 |
– |
} |
| 601 |
– |
|
| 621 |
|
if (avgSky < 0) /* no matrix output? */ |
| 622 |
|
goto alldone; |
| 623 |
|
|
| 711 |
|
int index; /* Category index */ |
| 712 |
|
double norm_diff_illum; /* Normalized diffuse illuimnance */ |
| 713 |
|
int i; |
| 714 |
< |
|
| 714 |
> |
|
| 715 |
|
/* Calculate atmospheric precipitable water content */ |
| 716 |
|
apwc = CalcPrecipWater(dew_point); |
| 717 |
|
|
