| 128 |
|
extern void CalcSkyPatchLumin( float *parr ); |
| 129 |
|
extern void ComputeSky( float *parr ); |
| 130 |
|
|
| 131 |
+ |
|
| 132 |
+ |
extern double solar_sunset(int month, int day); |
| 133 |
+ |
extern double solar_sunrise(int month, int day); |
| 134 |
+ |
|
| 135 |
|
/* Degrees into radians */ |
| 136 |
|
#define DegToRad(deg) ((deg)*(PI/180.)) |
| 137 |
|
|
| 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]; |
| 335 |
|
/* get options */ |
| 431 |
|
case 'A': /* compute average sky */ |
| 432 |
|
avgSky = 1; |
| 433 |
|
break; |
| 434 |
+ |
case 'i': |
| 435 |
+ |
timeinterval = atof(argv[++i]); |
| 436 |
+ |
break; |
| 437 |
|
default: |
| 438 |
|
goto userr; |
| 439 |
|
} |
| 507 |
|
mtx_data = resize_dmatrix(mtx_data, tstorage=2, nskypatch); |
| 508 |
|
/* process each time step in tape */ |
| 509 |
|
while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) { |
| 510 |
< |
double sda, sta; |
| 510 |
> |
double sda, sta, st; |
| 511 |
|
int sun_in_sky; |
| 512 |
|
/* compute solar position */ |
| 513 |
|
if ((mo == 2) & (da == 29)) { |
| 517 |
|
julian_date = jdate(mo, da) + leap_day; |
| 518 |
|
sda = sdec(julian_date); |
| 519 |
|
sta = stadj(julian_date); |
| 520 |
< |
altitude = salt(sda, hr+sta); |
| 520 |
> |
st = hr + sta; |
| 521 |
> |
|
| 522 |
> |
if (timeinterval > 0) { |
| 523 |
> |
if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120) |
| 524 |
> |
st = (st + timeinterval/120 + solar_sunrise(mo, da))/2; |
| 525 |
> |
else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120) |
| 526 |
> |
st = (st - timeinterval/120 + solar_sunset(mo, da))/2; |
| 527 |
> |
} |
| 528 |
> |
altitude = salt(sda, st); |
| 529 |
|
sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG/2.)); |
| 530 |
|
if (sun_hours_only && !sun_in_sky) |
| 531 |
|
continue; /* skipping nighttime points */ |
| 532 |
< |
azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation); |
| 532 |
> |
azimuth = sazi(sda, st) + PI - DegToRad(rotation); |
| 533 |
|
|
| 534 |
|
mtx_offset = 3*nskypatch*nstored; |
| 535 |
|
nstored += !avgSky | !nstored; |
| 680 |
|
exit(1); |
| 681 |
|
} |
| 682 |
|
|
| 683 |
+ |
|
| 684 |
|
/* Return maximum of two doubles */ |
| 685 |
|
double dmax( double a, double b ) |
| 686 |
|
{ return (a > b) ? a : b; } |
| 726 |
|
/* Limit sky clearness */ |
| 727 |
|
if (sky_clearness > 11.9) |
| 728 |
|
sky_clearness = 11.9; |
| 729 |
+ |
else if (sky_clearness < 1.0) |
| 730 |
+ |
sky_clearness = 1.0; |
| 731 |
|
|
| 732 |
|
/* Limit sky brightness */ |
| 733 |
|
if (sky_brightness < 0.01) |
| 734 |
|
sky_brightness = 0.01; |
| 735 |
+ |
else if (sky_brightness > 0.6) |
| 736 |
+ |
sky_brightness = 0.6; |
| 737 |
|
|
| 738 |
|
/* Calculate illuminance */ |
| 739 |
|
index = GetCategoryIndex(); |
| 786 |
|
} |
| 787 |
|
} |
| 788 |
|
|
| 789 |
+ |
|
| 790 |
+ |
double |
| 791 |
+ |
solar_sunset(int month, int day) |
| 792 |
+ |
{ |
| 793 |
+ |
float W; |
| 794 |
+ |
W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day)))); |
| 795 |
+ |
return(12 + (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15)); |
| 796 |
+ |
} |
| 797 |
+ |
|
| 798 |
+ |
|
| 799 |
+ |
double |
| 800 |
+ |
solar_sunrise(int month, int day) |
| 801 |
+ |
{ |
| 802 |
+ |
float W; |
| 803 |
+ |
W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day)))); |
| 804 |
+ |
return(12 - (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15)); |
| 805 |
+ |
} |
| 806 |
+ |
|
| 807 |
+ |
|
| 808 |
|
/* Add in solar direct to nearest sky patches (GW) */ |
| 809 |
|
void |
| 810 |
|
AddDirect(float *parr) |
| 1074 |
|
if (sky_brightness < 0.01) |
| 1075 |
|
sky_brightness = 0.01; |
| 1076 |
|
|
| 1077 |
+ |
if (sky_clearness < 1.0000) |
| 1078 |
+ |
{ |
| 1079 |
+ |
sky_clearness = 1.0000; |
| 1080 |
+ |
} |
| 1081 |
+ |
|
| 1082 |
+ |
if (sky_brightness > 0.6) |
| 1083 |
+ |
{ |
| 1084 |
+ |
sky_brightness = 0.6; |
| 1085 |
+ |
} |
| 1086 |
+ |
|
| 1087 |
|
while (((fabs(diff_irrad - test1) > 10.0) || |
| 1088 |
|
(fabs(dir_irrad - test2) > 10.0)) && !(counter == 5)) |
| 1089 |
|
{ |
| 1109 |
|
/* Limit sky brightness */ |
| 1110 |
|
if (sky_brightness < 0.01) |
| 1111 |
|
sky_brightness = 0.01; |
| 1112 |
+ |
|
| 1113 |
+ |
if (sky_clearness < 1.0000) |
| 1114 |
+ |
{ |
| 1115 |
+ |
sky_clearness = 1.0000; |
| 1116 |
+ |
} |
| 1117 |
+ |
|
| 1118 |
+ |
if (sky_brightness > 0.6) |
| 1119 |
+ |
{ |
| 1120 |
+ |
sky_brightness = 0.6; |
| 1121 |
+ |
} |
| 1122 |
|
} |
| 1123 |
|
|
| 1124 |
|
return GetCategoryIndex(); |