4 |
|
* Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France |
5 |
|
* *BOUYGUES |
6 |
|
* 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France |
7 |
+ |
* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 |
8 |
|
*/ |
9 |
|
|
10 |
< |
#define _USE_MATH_DEFINES |
10 |
> |
#define _USE_MATH_DEFINES |
11 |
|
#include <stdio.h> |
12 |
|
#include <string.h> |
13 |
|
#include <math.h> |
54 |
|
90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, |
55 |
|
80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, |
56 |
|
240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0}; |
57 |
+ |
/* default values for Berlin */ |
58 |
+ |
float locus[] = { |
59 |
+ |
-4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786}; |
60 |
|
|
61 |
|
|
62 |
|
|
63 |
< |
/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */ |
63 |
> |
/* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */ |
64 |
|
double sky_brightness(); |
65 |
|
double sky_clearness(); |
66 |
|
|
93 |
|
void check_sun_position(); |
94 |
|
void computesky(); |
95 |
|
void printhead(int ac, char** av); |
96 |
< |
void userror(char* msg); |
96 |
> |
void usage_error(char* msg); |
97 |
|
void printsky(); |
98 |
|
|
99 |
|
FILE * frlibopen(char* fname); |
102 |
|
double get_eccentricity(); |
103 |
|
double air_mass(); |
104 |
|
|
105 |
< |
double solar_sunset(); |
106 |
< |
double solar_sunrise(); |
105 |
> |
double solar_sunset(int month, int day); |
106 |
> |
double solar_sunrise(int month, int day); |
107 |
|
double stadj(); |
108 |
|
int jdate(int month, int day); |
109 |
|
|
106 |
– |
|
107 |
– |
|
108 |
– |
/* sun calculation constants */ |
109 |
– |
extern double s_latitude; |
110 |
– |
extern double s_longitude; |
111 |
– |
extern double s_meridian; |
112 |
– |
|
110 |
|
const double AU = 149597890E3; |
111 |
|
const double solar_constant_e = 1367; /* solar constant W/m^2 */ |
112 |
< |
const double solar_constant_l = 127.5; /* solar constant klux */ |
112 |
> |
const double solar_constant_l = 127500; /* solar constant lux */ |
113 |
|
|
114 |
|
const double half_sun_angle = 0.2665; |
115 |
|
const double half_direct_angle = 2.85; |
116 |
|
|
117 |
< |
const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */ |
118 |
< |
const double skyclearsup = 12.1; |
117 |
> |
const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */ |
118 |
> |
const double skyclearsup = 12.01; |
119 |
|
const double skybriginf = 0.01; |
120 |
|
const double skybrigsup = 0.6; |
121 |
|
|
144 |
|
int output=0; /* define the unit of the output (sky luminance or radiance): */ |
145 |
|
/* visible watt=0, solar watt=1, lumen=2 */ |
146 |
|
int input=0; /* define the input for the calulation */ |
147 |
< |
|
147 |
> |
int color_output=0; |
148 |
|
int suppress_warnings=0; |
149 |
|
|
150 |
|
/* default values */ |
155 |
|
double gprefl = 0.2; |
156 |
|
int S_INTER=0; |
157 |
|
|
158 |
+ |
|
159 |
|
/* computed values */ |
160 |
|
double sundir[3]; |
161 |
|
double groundbr = 0; |
167 |
|
char *progname; |
168 |
|
char errmsg[128]; |
169 |
|
|
170 |
+ |
double st; |
171 |
|
|
172 |
|
|
174 |
– |
|
173 |
|
int main(int argc, char** argv) |
174 |
|
{ |
175 |
|
int i; |
180 |
|
return 0; |
181 |
|
} |
182 |
|
if (argc < 4) |
183 |
< |
userror("arg count"); |
183 |
> |
usage_error("arg count"); |
184 |
|
if (!strcmp(argv[1], "-ang")) { |
185 |
|
altitude = atof(argv[2]) * (M_PI/180); |
186 |
|
azimuth = atof(argv[3]) * (M_PI/180); |
188 |
|
} else { |
189 |
|
month = atoi(argv[1]); |
190 |
|
if (month < 1 || month > 12) |
191 |
< |
userror("bad month"); |
191 |
> |
usage_error("bad month"); |
192 |
|
day = atoi(argv[2]); |
193 |
|
if (day < 1 || day > 31) |
194 |
< |
userror("bad day"); |
194 |
> |
usage_error("bad day"); |
195 |
|
hour = atof(argv[3]); |
196 |
|
if (hour < 0 || hour >= 24) |
197 |
< |
userror("bad hour"); |
197 |
> |
usage_error("bad hour"); |
198 |
|
tsolar = argv[3][0] == '+'; |
199 |
|
} |
200 |
|
for (i = 4; i < argc; i++) |
212 |
|
cloudy = argv[i][0] == '+' ? 2 : 1; |
213 |
|
dosun = 0; |
214 |
|
break; |
215 |
+ |
case 'C': |
216 |
+ |
if (argv[i][2] == 'I' && argv[i][3] == 'E' ) { |
217 |
+ |
locus[0] = -4.607e9; |
218 |
+ |
locus[1] = 2.9678e6; |
219 |
+ |
locus[2] = 0.09911e3; |
220 |
+ |
locus[3] = 0.244063; |
221 |
+ |
locus[4] = -2.0064e9; |
222 |
+ |
locus[5] = 1.9018e6; |
223 |
+ |
locus[6] = 0.24748e3; |
224 |
+ |
locus[7] = 0.23704; |
225 |
+ |
locus[8] = -3.0; |
226 |
+ |
locus[9] = 2.87; |
227 |
+ |
locus[10] = -0.275; |
228 |
+ |
}else{ color_output = 1; |
229 |
+ |
} |
230 |
+ |
break; |
231 |
+ |
case 'l': |
232 |
+ |
locus[0] = atof(argv[++i]); |
233 |
+ |
locus[1] = atof(argv[++i]); |
234 |
+ |
locus[2] = atof(argv[++i]); |
235 |
+ |
locus[3] = atof(argv[++i]); |
236 |
+ |
locus[4] = locus[0]; |
237 |
+ |
locus[5] = locus[1]; |
238 |
+ |
locus[6] = locus[2]; |
239 |
+ |
locus[7] = locus[3]; |
240 |
+ |
locus[8] = atof(argv[++i]); |
241 |
+ |
locus[9] = atof(argv[++i]); |
242 |
+ |
locus[10] = atof(argv[++i]); |
243 |
+ |
break; |
244 |
+ |
|
245 |
|
case 't': |
246 |
|
betaturbidity = atof(argv[++i]); |
247 |
|
break; |
265 |
|
break; |
266 |
|
|
267 |
|
case 'O': |
268 |
< |
output = atof(argv[++i]); /*define the unit of the output of the program : |
269 |
< |
sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */ |
268 |
> |
output = atof(argv[++i]); /*define the unit of the output of the program: |
269 |
> |
sky and sun luminance/radiance |
270 |
> |
(0==W visible, 1==W solar radiation, 2==lm) */ |
271 |
|
break; |
272 |
|
|
273 |
|
case 'P': |
299 |
|
globalirradiance = atof(argv[++i]); |
300 |
|
break; |
301 |
|
|
273 |
– |
/* |
274 |
– |
case 'l': |
275 |
– |
sunaltitude_border = atof(argv[++i]); |
276 |
– |
break; |
277 |
– |
*/ |
278 |
– |
|
302 |
|
case 'i': |
303 |
|
timeinterval = atof(argv[++i]); |
304 |
|
break; |
306 |
|
|
307 |
|
default: |
308 |
|
sprintf(errmsg, "unknown option: %s", argv[i]); |
309 |
< |
userror(errmsg); |
309 |
> |
usage_error(errmsg); |
310 |
|
} |
311 |
|
else |
312 |
< |
userror("bad option"); |
312 |
> |
usage_error("bad option"); |
313 |
|
|
314 |
< |
if (fabs(s_meridian-s_longitude) > 30*M_PI/180) |
315 |
< |
fprintf(stderr, |
293 |
< |
"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
314 |
> |
if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180) |
315 |
> |
fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
316 |
|
progname, (s_longitude-s_meridian)*12/M_PI); |
317 |
|
|
318 |
|
|
319 |
|
/* dynamic memory allocation for the pointers */ |
298 |
– |
|
320 |
|
if ( (c_perez = calloc(5, sizeof(double))) == NULL ) |
321 |
< |
{ |
301 |
< |
fprintf(stderr,"Out of memory error in function main"); |
302 |
< |
return 1; |
303 |
< |
} |
321 |
> |
{ fprintf(stderr,"Out of memory error in function main"); return 1; } |
322 |
|
|
323 |
+ |
|
324 |
|
printhead(argc, argv); |
325 |
|
computesky(); |
307 |
– |
|
308 |
– |
if(*(c_perez+1)>0) |
309 |
– |
{ |
310 |
– |
fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1)); |
311 |
– |
print_error_sky(); |
312 |
– |
exit(1); |
313 |
– |
} |
314 |
– |
|
326 |
|
printsky(); |
327 |
|
return 0; |
328 |
+ |
|
329 |
|
} |
330 |
|
|
331 |
|
|
332 |
|
|
333 |
|
|
334 |
|
|
323 |
– |
|
324 |
– |
|
325 |
– |
|
335 |
|
void computesky() |
336 |
|
{ |
337 |
|
|
348 |
|
|
349 |
|
if (month) { /* from date and time */ |
350 |
|
int jd; |
351 |
< |
double sd, st; |
351 |
> |
double sd; |
352 |
|
|
353 |
|
jd = jdate(month, day); /* Julian date */ |
354 |
|
sd = sdec(jd); /* solar declination */ |
356 |
|
st = hour; |
357 |
|
else |
358 |
|
st = hour + stadj(jd); |
350 |
– |
|
359 |
|
|
352 |
– |
if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) { |
353 |
– |
print_error_sky(); |
354 |
– |
exit(1); |
355 |
– |
} |
356 |
– |
|
360 |
|
|
361 |
|
if(timeinterval) { |
362 |
|
|
364 |
|
fprintf(stderr, "time interval negative\n"); |
365 |
|
exit(1); |
366 |
|
} |
367 |
< |
|
368 |
< |
if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) { |
369 |
< |
|
370 |
< |
fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour); |
371 |
< |
st= (st+timeinterval/120+solar_sunrise(month,day))/2; |
367 |
> |
|
368 |
> |
if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) { |
369 |
> |
st= (st+timeinterval/120+solar_sunrise(month,day))/2; |
370 |
> |
if(suppress_warnings==0) |
371 |
> |
{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } |
372 |
|
} |
373 |
|
|
374 |
< |
if(fabs(solar_sunset(month,day)-st)<timeinterval/60) { |
375 |
< |
fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour); |
376 |
< |
st= (st-timeinterval/120+solar_sunset(month,day))/2; |
374 |
> |
if(fabs(solar_sunset(month,day)-st)<timeinterval/120) { |
375 |
> |
st= (st-timeinterval/120+solar_sunset(month,day))/2; |
376 |
> |
if(suppress_warnings==0) |
377 |
> |
{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } |
378 |
|
} |
379 |
+ |
|
380 |
+ |
if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) { |
381 |
+ |
if(suppress_warnings==0) |
382 |
+ |
{ fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); } |
383 |
+ |
altitude = salt(sd, st); |
384 |
+ |
azimuth = sazi(sd, st); |
385 |
+ |
print_error_sky(); |
386 |
+ |
exit(0); |
387 |
+ |
} |
388 |
|
} |
389 |
+ |
else |
390 |
|
|
391 |
+ |
if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) { |
392 |
+ |
if(suppress_warnings==0) |
393 |
+ |
{ fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); } |
394 |
+ |
altitude = salt(sd, st); |
395 |
+ |
azimuth = sazi(sd, st); |
396 |
+ |
print_error_sky(); |
397 |
+ |
exit(0); |
398 |
+ |
} |
399 |
|
|
400 |
|
altitude = salt(sd, st); |
401 |
|
azimuth = sazi(sd, st); |
407 |
|
|
408 |
|
|
409 |
|
|
388 |
– |
/* if loop for the -l option. W.Sprenger (01/2013) */ |
389 |
– |
/* |
390 |
– |
if (altitude*180/M_PI < sunaltitude_border) { |
391 |
– |
|
392 |
– |
if (suppress_warnings==0) { |
393 |
– |
fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border); |
394 |
– |
} |
395 |
– |
print_error_sky(); |
396 |
– |
exit(1); |
397 |
– |
} |
398 |
– |
*/ |
399 |
– |
|
410 |
|
|
411 |
|
if (!cloudy && altitude > 87.*M_PI/180.) { |
412 |
|
|
418 |
|
altitude = 87.*M_PI/180.; |
419 |
|
} |
420 |
|
|
421 |
+ |
|
422 |
+ |
|
423 |
|
sundir[0] = -sin(azimuth)*cos(altitude); |
424 |
|
sundir[1] = -cos(azimuth)*cos(altitude); |
425 |
|
sundir[2] = sin(altitude); |
427 |
|
|
428 |
|
/* calculation for the new functions */ |
429 |
|
sunzenith = 90 - altitude*180/M_PI; |
430 |
< |
|
419 |
< |
|
430 |
> |
|
431 |
|
|
432 |
|
/* compute the inputs for the calculation of the light distribution over the sky*/ |
433 |
|
if (input==0) /* P */ |
451 |
|
check_irradiances(); |
452 |
|
skybrightness = sky_brightness(); |
453 |
|
skyclearness = sky_clearness(); |
454 |
+ |
|
455 |
|
check_parametrization(); |
456 |
< |
|
456 |
> |
|
457 |
|
if (output==0 || output==2) |
458 |
|
{ |
459 |
|
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
477 |
|
if (altitude<=0) |
478 |
|
{ |
479 |
|
if (suppress_warnings==0) |
480 |
< |
fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n"); |
480 |
> |
fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n"); |
481 |
|
directirradiance = 0; |
482 |
|
diffuseirradiance = 0; |
483 |
|
} else { |
517 |
|
{ |
518 |
|
if (suppress_warnings==0) |
519 |
|
fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n"); |
508 |
– |
|
520 |
|
globalirradiance=erbs_s0*0.999; |
521 |
|
} |
522 |
|
|
529 |
|
directirradiance=globalirradiance-diffuseirradiance; |
530 |
|
|
531 |
|
printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance); |
532 |
< |
printf("# WARNING: the -E option is only recommended for a rough estimation!"); |
532 |
> |
printf("# WARNING: the -E option is only recommended for a rough estimation!\n"); |
533 |
|
|
534 |
|
directirradiance=directirradiance/sin(altitude); |
535 |
|
|
552 |
|
|
553 |
|
|
554 |
|
|
555 |
< |
else {fprintf(stderr,"error at the input arguments"); exit(1);} |
555 |
> |
else { fprintf(stderr,"error at the input arguments"); exit(1); } |
556 |
|
|
557 |
|
|
558 |
|
|
574 |
|
|
575 |
|
|
576 |
|
|
566 |
– |
|
567 |
– |
|
577 |
|
/*calculation of the modelled luminance */ |
578 |
|
for (j=0;j<145;j++) |
579 |
|
{ |
619 |
|
else |
620 |
|
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
621 |
|
|
613 |
– |
|
622 |
|
|
623 |
|
|
624 |
< |
/* Compute the ground radiance */ |
625 |
< |
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
626 |
< |
zenithbr*=diffnormalization; |
624 |
> |
/* Compute the ground radiance */ |
625 |
> |
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
626 |
> |
zenithbr*=diffnormalization; |
627 |
|
|
628 |
< |
if (skyclearness==1) |
628 |
> |
if (skyclearness==1) |
629 |
|
normfactor = 0.777778; |
630 |
|
|
631 |
< |
if (skyclearness>=6) |
631 |
> |
if (skyclearness>=6) |
632 |
|
{ |
633 |
|
F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); |
634 |
|
normfactor = normsc()/F2/M_PI; |
635 |
|
} |
636 |
|
|
637 |
< |
if ( (skyclearness>1) && (skyclearness<6) ) |
637 |
> |
if ( (skyclearness>1) && (skyclearness<6) ) |
638 |
|
{ |
639 |
|
S_INTER=1; |
640 |
|
F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude)); |
641 |
|
normfactor = normsc()/F2/M_PI; |
642 |
|
} |
643 |
|
|
644 |
< |
groundbr = zenithbr*normfactor; |
644 |
> |
groundbr = zenithbr*normfactor; |
645 |
|
|
646 |
< |
if (dosun&&(skyclearness>1)) |
646 |
> |
if (dosun&&(skyclearness>1)) |
647 |
|
groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; |
648 |
|
|
649 |
< |
groundbr *= gprefl; |
649 |
> |
groundbr *= gprefl; |
650 |
|
|
651 |
|
|
652 |
+ |
|
653 |
+ |
if(*(c_perez+1)>0) |
654 |
+ |
{ |
655 |
+ |
if(suppress_warnings==0) |
656 |
+ |
{ fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));} |
657 |
+ |
print_error_sky(); |
658 |
+ |
exit(0); |
659 |
+ |
} |
660 |
|
|
661 |
+ |
|
662 |
|
return; |
663 |
|
} |
664 |
|
|
665 |
|
|
666 |
|
|
667 |
|
|
651 |
– |
void print_error_sky() |
652 |
– |
{ |
653 |
– |
sundir[0] = -sin(azimuth)*cos(altitude); |
654 |
– |
sundir[1] = -cos(azimuth)*cos(altitude); |
655 |
– |
sundir[2] = sin(altitude); |
656 |
– |
|
657 |
– |
printf("\nvoid brightfunc skyfunc\n"); |
658 |
– |
printf("2 skybright perezlum.cal\n"); |
659 |
– |
printf("0\n"); |
660 |
– |
printf("10 0.00 0.00 0.000 0.000 0.000 0.000 0.000 %f %f %f \n", sundir[0], sundir[1], sundir[2]); |
661 |
– |
} |
662 |
– |
|
668 |
|
|
664 |
– |
|
665 |
– |
|
666 |
– |
|
667 |
– |
|
668 |
– |
|
669 |
|
double solar_sunset(int month,int day) |
670 |
|
{ |
671 |
|
float W; |
675 |
|
} |
676 |
|
|
677 |
|
|
678 |
+ |
|
679 |
+ |
|
680 |
|
double solar_sunrise(int month,int day) |
681 |
|
{ |
682 |
|
float W; |
688 |
|
|
689 |
|
|
690 |
|
|
691 |
+ |
void printsky() |
692 |
+ |
{ |
693 |
+ |
|
694 |
+ |
printf("# Local solar time: %.2f\n", st); |
695 |
+ |
printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); |
696 |
|
|
697 |
|
|
691 |
– |
|
692 |
– |
|
693 |
– |
|
694 |
– |
|
695 |
– |
|
696 |
– |
void printsky() /* print out sky */ |
697 |
– |
{ |
698 |
|
if (dosun&&(skyclearness>1)) |
699 |
|
{ |
700 |
|
printf("\nvoid light solar\n"); |
711 |
|
printf("0\n0\n"); |
712 |
|
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
713 |
|
} |
714 |
+ |
/* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */ |
715 |
+ |
if (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 ) |
716 |
+ |
{ |
717 |
+ |
fprintf(stderr, " warning: sky clearness(epsilon)= %f \n",skyclearness); |
718 |
+ |
fprintf(stderr, " warning: intermediate sky!! \n"); |
719 |
+ |
fprintf(stderr, " warning: color model for intermediate sky pending \n"); |
720 |
+ |
fprintf(stderr, " warning: no color output ! \n"); |
721 |
+ |
color_output=0; |
722 |
+ |
} |
723 |
+ |
if (color_output==1) |
724 |
+ |
{ |
725 |
+ |
printf("\nvoid colorfunc skyfunc\n"); |
726 |
+ |
printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n"); |
727 |
+ |
printf("0\n"); |
728 |
+ |
printf("22 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", diffnormalization, groundbr, |
729 |
+ |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
730 |
+ |
sundir[0], sundir[1], sundir[2],skyclearness,locus[0],locus[1],locus[2],locus[3],locus[4],locus[5],locus[6],locus[7],locus[8],locus[9],locus[10]); |
731 |
+ |
}else{ |
732 |
+ |
printf("\nvoid brightfunc skyfunc\n"); |
733 |
+ |
printf("2 skybright perezlum.cal\n"); |
734 |
+ |
printf("0\n"); |
735 |
+ |
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
736 |
+ |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
737 |
+ |
sundir[0], sundir[1], sundir[2]); |
738 |
+ |
} |
739 |
|
|
740 |
+ |
} |
741 |
|
|
742 |
+ |
|
743 |
+ |
|
744 |
+ |
void print_error_sky() |
745 |
+ |
{ |
746 |
+ |
|
747 |
+ |
|
748 |
+ |
sundir[0] = -sin(azimuth)*cos(altitude); |
749 |
+ |
sundir[1] = -cos(azimuth)*cos(altitude); |
750 |
+ |
sundir[2] = sin(altitude); |
751 |
+ |
|
752 |
+ |
printf("# Local solar time: %.2f\n", st); |
753 |
+ |
printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); |
754 |
+ |
|
755 |
|
printf("\nvoid brightfunc skyfunc\n"); |
756 |
|
printf("2 skybright perezlum.cal\n"); |
757 |
|
printf("0\n"); |
758 |
< |
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
720 |
< |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
721 |
< |
sundir[0], sundir[1], sundir[2]); |
722 |
< |
|
758 |
> |
printf("10 0.00 0.00 0.000 0.000 0.000 0.000 0.000 %f %f %f \n", sundir[0], sundir[1], sundir[2]); |
759 |
|
} |
760 |
+ |
|
761 |
|
|
762 |
|
|
763 |
+ |
|
764 |
+ |
|
765 |
|
void printdefaults() /* print default values */ |
766 |
|
{ |
767 |
|
printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl); |
775 |
|
} |
776 |
|
|
777 |
|
|
778 |
< |
void userror(char* msg) /* print usage error and quit */ |
778 |
> |
|
779 |
> |
|
780 |
> |
void usage_error(char* msg) /* print usage error and quit */ |
781 |
|
{ |
782 |
|
if (msg != NULL) |
783 |
|
fprintf(stderr, "%s: Use error - %s\n\n", progname, msg); |
794 |
|
fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n"); |
795 |
|
fprintf(stderr, " Output specification with option:\n"); |
796 |
|
fprintf(stderr, " -O [0|1|2] (0=output in W/m^2/sr visible, 1=output in W/m^2/sr solar, 2=output in candela/m^2), default is 0 \n"); |
797 |
< |
fprintf(stderr, " gendaylit version 2.3 (2013/08/08) \n\n"); |
797 |
> |
fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n"); |
798 |
|
exit(1); |
799 |
|
} |
800 |
|
|
801 |
|
|
802 |
|
|
803 |
+ |
|
804 |
|
double normsc() /* compute normalization factor (E0*F2/L0) */ |
805 |
|
{ |
806 |
|
static double nfc[2][5] = { |
824 |
|
|
825 |
|
|
826 |
|
|
827 |
+ |
|
828 |
+ |
|
829 |
|
void printhead(int ac, char** av) /* print command header */ |
830 |
|
{ |
831 |
|
putchar('#'); |
839 |
|
|
840 |
|
|
841 |
|
|
842 |
+ |
|
843 |
+ |
|
844 |
|
/* Perez models */ |
845 |
|
|
846 |
|
/* Perez global horizontal luminous efficacy model */ |
912 |
|
|
913 |
|
|
914 |
|
|
869 |
– |
|
915 |
|
for (i=1; i<=category_total_number; i++) |
916 |
|
{ |
917 |
|
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
925 |
|
} |
926 |
|
|
927 |
|
|
928 |
+ |
|
929 |
+ |
|
930 |
|
/* global horizontal diffuse efficacy model, according to PEREZ */ |
931 |
|
double glob_h_diffuse_effi_PEREZ() |
932 |
|
{ |
934 |
|
double category_bounds[10], a[10], b[10], c[10], d[10]; |
935 |
|
int category_total_number, category_number, i; |
936 |
|
|
890 |
– |
|
891 |
– |
|
892 |
– |
|
937 |
|
check_parametrization(); |
938 |
|
|
939 |
|
|
940 |
|
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
941 |
< |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */ |
941 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */ |
942 |
|
|
899 |
– |
|
900 |
– |
|
943 |
|
/* initialize category bounds (clearness index bounds) */ |
944 |
|
|
945 |
|
category_total_number = 8; |
995 |
|
|
996 |
|
|
997 |
|
|
956 |
– |
|
998 |
|
category_number = -1; |
999 |
|
for (i=1; i<=category_total_number; i++) |
1000 |
|
{ |
1004 |
|
|
1005 |
|
if (category_number == -1) { |
1006 |
|
if (suppress_warnings==0) |
1007 |
< |
fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness); |
1007 |
> |
fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness); |
1008 |
|
print_error_sky(); |
1009 |
< |
exit(1); |
1009 |
> |
exit(0); |
1010 |
|
} |
1011 |
|
|
1012 |
|
|
1019 |
|
|
1020 |
|
|
1021 |
|
|
1022 |
+ |
|
1023 |
+ |
|
1024 |
+ |
|
1025 |
|
/* direct normal efficacy model, according to PEREZ */ |
1026 |
|
|
1027 |
|
double direct_n_effi_PEREZ() |
1032 |
|
int category_total_number, category_number, i; |
1033 |
|
|
1034 |
|
|
1035 |
< |
if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
1036 |
< |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n"); |
1035 |
> |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
1036 |
> |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/ |
1037 |
|
|
1038 |
|
|
1039 |
|
/* initialize category bounds (clearness index bounds) */ |
1107 |
|
/*check the range of epsilon and delta indexes of the perez parametrization*/ |
1108 |
|
void check_parametrization() |
1109 |
|
{ |
1110 |
+ |
|
1111 |
|
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) |
1112 |
|
{ |
1113 |
|
|
1114 |
|
/* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */ |
1115 |
+ |
|
1116 |
|
if (skyclearness<skyclearinf){ |
1117 |
< |
if (suppress_warnings==0) |
1118 |
< |
/* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */ |
1117 |
> |
/* if (suppress_warnings==0) |
1118 |
> |
fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */ |
1119 |
|
skyclearness=skyclearinf; |
1120 |
|
} |
1121 |
|
if (skyclearness>skyclearsup){ |
1122 |
< |
if (suppress_warnings==0) |
1123 |
< |
/* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ |
1124 |
< |
skyclearness=skyclearsup-0.1; |
1122 |
> |
/* if (suppress_warnings==0) |
1123 |
> |
fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ |
1124 |
> |
skyclearness=skyclearsup-0.001; |
1125 |
|
} |
1126 |
|
if (skybrightness<skybriginf){ |
1127 |
< |
if (suppress_warnings==0) |
1128 |
< |
/* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */ |
1127 |
> |
/* if (suppress_warnings==0) |
1128 |
> |
fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */ |
1129 |
|
skybrightness=skybriginf; |
1130 |
|
} |
1131 |
|
if (skybrightness>skybrigsup){ |
1132 |
< |
if (suppress_warnings==0) |
1133 |
< |
/* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ |
1132 |
> |
/* if (suppress_warnings==0) |
1133 |
> |
fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ |
1134 |
|
skybrightness=skybrigsup; |
1135 |
|
} |
1136 |
|
|
1139 |
|
} |
1140 |
|
|
1141 |
|
|
1142 |
+ |
|
1143 |
+ |
|
1144 |
+ |
|
1145 |
|
/* validity of the direct and diffuse components */ |
1146 |
|
void check_illuminances() |
1147 |
|
{ |
1148 |
|
if (directilluminance < 0) { |
1149 |
< |
fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n"); |
1149 |
> |
if(suppress_warnings==0) |
1150 |
> |
{ fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); } |
1151 |
|
directilluminance = 0.0; |
1152 |
|
} |
1153 |
|
if (diffuseilluminance < 0) { |
1154 |
< |
fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n"); |
1154 |
> |
if(suppress_warnings==0) |
1155 |
> |
{ fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); } |
1156 |
|
diffuseilluminance = 0.0; |
1157 |
|
} |
1158 |
< |
if (directilluminance > solar_constant_l*1000.0) { |
1159 |
< |
fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n"); |
1160 |
< |
exit(1); |
1158 |
> |
|
1159 |
> |
if (directilluminance+diffuseilluminance==0 && altitude > 0) { |
1160 |
> |
if(suppress_warnings==0) |
1161 |
> |
{ fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); } |
1162 |
> |
print_error_sky(); |
1163 |
> |
exit(0); |
1164 |
|
} |
1165 |
+ |
|
1166 |
+ |
if (directilluminance > solar_constant_l) { |
1167 |
+ |
if(suppress_warnings==0) |
1168 |
+ |
{ fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); } |
1169 |
+ |
print_error_sky(); |
1170 |
+ |
exit(0); |
1171 |
+ |
} |
1172 |
|
} |
1173 |
|
|
1174 |
|
|
1175 |
|
void check_irradiances() |
1176 |
|
{ |
1177 |
|
if (directirradiance < 0) { |
1178 |
< |
fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n"); |
1178 |
> |
if(suppress_warnings==0) |
1179 |
> |
{ fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); } |
1180 |
|
directirradiance = 0.0; |
1181 |
|
} |
1182 |
|
if (diffuseirradiance < 0) { |
1183 |
< |
fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n"); |
1183 |
> |
if(suppress_warnings==0) |
1184 |
> |
{ fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); } |
1185 |
|
diffuseirradiance = 0.0; |
1186 |
|
} |
1187 |
+ |
|
1188 |
+ |
if (directirradiance+diffuseirradiance==0 && altitude > 0) { |
1189 |
+ |
if(suppress_warnings==0) |
1190 |
+ |
{ fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); } |
1191 |
+ |
print_error_sky(); |
1192 |
+ |
exit(0); |
1193 |
+ |
} |
1194 |
+ |
|
1195 |
|
if (directirradiance > solar_constant_e) { |
1196 |
< |
fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n"); |
1197 |
< |
exit(1); |
1196 |
> |
if(suppress_warnings==0) |
1197 |
> |
{ fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); } |
1198 |
> |
print_error_sky(); |
1199 |
> |
exit(0); |
1200 |
|
} |
1201 |
|
} |
1202 |
|
|
1255 |
|
double test1=0.1, test2=0.1, d_eff; |
1256 |
|
int counter=0; |
1257 |
|
|
1258 |
< |
diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000); |
1259 |
< |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000); |
1258 |
> |
diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l); |
1259 |
> |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l); |
1260 |
|
skyclearness = sky_clearness(); |
1261 |
|
skybrightness = sky_brightness(); |
1262 |
|
check_parametrization(); |
1284 |
|
skybrightness = sky_brightness(); |
1285 |
|
skyclearness = sky_clearness(); |
1286 |
|
check_parametrization(); |
1214 |
– |
|
1215 |
– |
/*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/ |
1287 |
|
|
1288 |
|
} |
1289 |
|
|
1320 |
|
|
1321 |
|
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1322 |
|
{ |
1323 |
< |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n"); |
1323 |
> |
fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n"); |
1324 |
|
exit(1); |
1325 |
|
} |
1326 |
|
|
1371 |
|
|
1372 |
|
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1373 |
|
{ |
1374 |
< |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); |
1374 |
> |
fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n"); |
1375 |
|
exit(1); |
1376 |
|
} |
1377 |
|
|
1440 |
|
else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 ) |
1441 |
|
{ |
1442 |
|
printf("error in calculation of gamma (angle between point and sun"); |
1443 |
< |
exit(3); |
1443 |
> |
exit(1); |
1444 |
|
} |
1445 |
|
else |
1446 |
|
*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)); |
1484 |
|
double m; |
1485 |
|
if (sunzenith>90) |
1486 |
|
{ |
1487 |
< |
fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n"); |
1488 |
< |
exit(1); |
1487 |
> |
if(suppress_warnings==0) |
1488 |
> |
{ fprintf(stderr, "Warning: air mass has reached the maximal value\n"); } |
1489 |
> |
sunzenith=90; |
1490 |
|
} |
1491 |
|
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1492 |
|
return(m); |