21 |
|
#include <string.h> |
22 |
|
#include <math.h> |
23 |
|
#include <stdlib.h> |
24 |
+ |
#include <ctype.h> |
25 |
|
|
26 |
|
#include "rtio.h" |
27 |
|
#include "fvect.h" |
43 |
|
/* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */ |
44 |
|
double sky_brightness(); |
45 |
|
double sky_clearness(); |
46 |
+ |
void computesky(); |
47 |
|
|
48 |
|
/* calculation of the direct and diffuse components from the Perez parametrization */ |
49 |
|
double diffus_irradiance_from_sky_brightness(); |
73 |
|
double integ_lv(float *lv,float *theta); |
74 |
|
float *theta_ordered(char *filename); |
75 |
|
float *phi_ordered(char *filename); |
76 |
+ |
void skip_comments(FILE *fp); |
77 |
|
|
78 |
|
|
79 |
|
|
162 |
|
if (argc < 4) |
163 |
|
userror("arg count"); |
164 |
|
if (!strcmp(argv[1], "-ang")) { |
165 |
< |
altitude = atof(argv[2]) * (M_PI/180); |
166 |
< |
azimuth = atof(argv[3]) * (M_PI/180); |
165 |
> |
altitude = atof(argv[2]) * (PI/180); |
166 |
> |
azimuth = atof(argv[3]) * (PI/180); |
167 |
|
month = 0; |
168 |
|
} else { |
169 |
|
month = atoi(argv[1]); |
203 |
|
gprefl = atof(argv[++i]); |
204 |
|
break; |
205 |
|
case 'a': |
206 |
< |
s_latitude = atof(argv[++i]) * (M_PI/180); |
206 |
> |
s_latitude = atof(argv[++i]) * (PI/180); |
207 |
|
break; |
208 |
|
case 'o': |
209 |
< |
s_longitude = atof(argv[++i]) * (M_PI/180); |
209 |
> |
s_longitude = atof(argv[++i]) * (PI/180); |
210 |
|
break; |
211 |
|
case 'm': |
212 |
< |
s_meridian = atof(argv[++i]) * (M_PI/180); |
212 |
> |
s_meridian = atof(argv[++i]) * (PI/180); |
213 |
|
break; |
214 |
|
|
215 |
|
|
251 |
|
else |
252 |
|
userror("bad option"); |
253 |
|
|
254 |
< |
if (fabs(s_meridian-s_longitude) > 30*M_PI/180) |
254 |
> |
if (fabs(s_meridian-s_longitude) > 30*PI/180) |
255 |
|
fprintf(stderr, |
256 |
|
"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
257 |
< |
progname, (s_longitude-s_meridian)*12/M_PI); |
257 |
> |
progname, (s_longitude-s_meridian)*12/PI); |
258 |
|
|
259 |
|
|
260 |
|
/* allocation dynamique de memoire pour les pointeurs */ |
274 |
|
} |
275 |
|
|
276 |
|
|
277 |
+ |
void |
278 |
|
computesky() /* compute sky parameters */ |
279 |
|
{ |
280 |
|
|
281 |
|
/* new variables */ |
282 |
< |
int j, i; |
282 |
> |
int j; |
283 |
|
float *lv_mod; /* 145 luminance values*/ |
284 |
|
/* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */ |
285 |
|
float *theta_o, *phi_o, *coeff_perez; |
286 |
|
double dzeta, gamma; |
283 |
– |
double diffusion; |
287 |
|
double normfactor; |
288 |
|
|
289 |
|
|
306 |
|
daynumber = (double)jdate(month, day); |
307 |
|
|
308 |
|
} |
309 |
< |
if (!cloudy && altitude > 87.*M_PI/180.) { |
309 |
> |
if (!cloudy && altitude > 87.*PI/180.) { |
310 |
|
fprintf(stderr, |
311 |
|
"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n", |
312 |
|
progname); |
313 |
|
printf( |
314 |
|
"# warning - sun too close to zenith, reducing altitude to 87 degrees\n"); |
315 |
< |
altitude = 87.*M_PI/180.; |
315 |
> |
altitude = 87.*PI/180.; |
316 |
|
} |
317 |
|
sundir[0] = -sin(azimuth)*cos(altitude); |
318 |
|
sundir[1] = -cos(azimuth)*cos(altitude); |
320 |
|
|
321 |
|
|
322 |
|
/* calculation for the new functions */ |
323 |
< |
sunzenith = 90 - altitude*180/M_PI; |
323 |
> |
sunzenith = 90 - altitude*180/PI; |
324 |
|
|
325 |
|
|
326 |
|
|
459 |
|
|
460 |
|
/* calculation for the solar source */ |
461 |
|
if (output==0) |
462 |
< |
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY; |
462 |
> |
solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)))/WHTEFFICACY; |
463 |
|
|
464 |
|
else if (output==1) |
465 |
< |
solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
465 |
> |
solarradiance = directirradiance/(2*PI*(1-cos(half_sun_angle*PI/180))); |
466 |
|
|
467 |
|
else |
468 |
< |
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
468 |
> |
solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180))); |
469 |
|
|
470 |
|
|
471 |
|
|
473 |
|
/* Compute the ground radiance */ |
474 |
|
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
475 |
|
zenithbr*=diffnormalization; |
476 |
+ |
/* |
477 |
|
fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr); |
478 |
< |
|
478 |
> |
*/ |
479 |
> |
|
480 |
|
if (skyclearness==1) |
481 |
|
normfactor = 0.777778; |
482 |
|
|
483 |
|
if (skyclearness>=6) |
484 |
|
{ |
485 |
< |
F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); |
486 |
< |
normfactor = normsc()/F2/M_PI; |
485 |
> |
F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); |
486 |
> |
normfactor = normsc()/F2/PI; |
487 |
|
} |
488 |
|
|
489 |
|
if ( (skyclearness>1) && (skyclearness<6) ) |
490 |
|
{ |
491 |
|
S_INTER=1; |
492 |
< |
F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude)); |
493 |
< |
normfactor = normsc()/F2/M_PI; |
492 |
> |
F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(PI/2.0-altitude)*(.4441+1.48*altitude)); |
493 |
> |
normfactor = normsc()/F2/PI; |
494 |
|
} |
495 |
|
|
496 |
|
groundbr = zenithbr*normfactor; |
497 |
|
printf("# Ground ambient level: %.1f\n", groundbr); |
498 |
|
|
499 |
|
if (dosun&&(skyclearness>1)) |
500 |
< |
groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; |
500 |
> |
groundbr += 6.8e-5/PI*solarradiance*sundir[2]; |
501 |
|
|
502 |
|
groundbr *= gprefl; |
503 |
|
|
551 |
|
printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr); |
552 |
|
else |
553 |
|
printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity); |
554 |
< |
printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI)); |
555 |
< |
printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI)); |
556 |
< |
printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI)); |
554 |
> |
printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI)); |
555 |
> |
printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI)); |
556 |
> |
printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI)); |
557 |
|
} |
558 |
|
|
559 |
|
|
588 |
|
register int i; |
589 |
|
/* polynomial approximation */ |
590 |
|
nf = nfc[S_INTER]; |
591 |
< |
x = (altitude - M_PI/4.0)/(M_PI/4.0); |
591 |
> |
x = (altitude - PI/4.0)/(PI/4.0); |
592 |
|
nsc = nf[i=4]; |
593 |
|
while (i--) |
594 |
|
nsc = nsc*x + nf[i]; |
613 |
|
|
614 |
|
|
615 |
|
|
616 |
+ |
void |
617 |
+ |
skip_comments(FILE *fp) /* skip comments in file */ |
618 |
+ |
{ |
619 |
+ |
int c; |
620 |
+ |
|
621 |
+ |
while ((c = getc(fp)) != EOF) |
622 |
+ |
if (c == '#') { |
623 |
+ |
while ((c = getc(fp)) != EOF) |
624 |
+ |
if (c == '\n') |
625 |
+ |
break; |
626 |
+ |
} else if (!isspace(c)) { |
627 |
+ |
ungetc(c, fp); |
628 |
+ |
break; |
629 |
+ |
} |
630 |
+ |
} |
631 |
|
|
632 |
|
|
633 |
|
|
614 |
– |
|
615 |
– |
|
616 |
– |
|
617 |
– |
|
618 |
– |
|
634 |
|
/* Perez models */ |
635 |
|
|
636 |
|
/* Perez global horizontal luminous efficacy model */ |
707 |
|
} |
708 |
|
|
709 |
|
value = a[category_number] + b[category_number]*atm_preci_water + |
710 |
< |
c[category_number]*cos(sunzenith*M_PI/180) + d[category_number]*log(skybrightness); |
710 |
> |
c[category_number]*cos(sunzenith*PI/180) + d[category_number]*log(skybrightness); |
711 |
|
|
712 |
|
return(value); |
713 |
|
} |
785 |
|
category_number = i; |
786 |
|
} |
787 |
|
|
788 |
< |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) + |
788 |
> |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*PI/180) + |
789 |
|
d[category_number]*log(skybrightness); |
790 |
|
|
791 |
|
return(value); |
866 |
|
category_number = i; |
867 |
|
} |
868 |
|
|
869 |
< |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*exp(5.73*sunzenith*M_PI/180 - 5) + d[category_number]*skybrightness; |
869 |
> |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*exp(5.73*sunzenith*PI/180 - 5) + d[category_number]*skybrightness; |
870 |
|
|
871 |
|
if (value < 0) value = 0; |
872 |
|
|
926 |
|
{ |
927 |
|
double value; |
928 |
|
|
929 |
< |
value = ( (diffusirradiance + directirradiance)/(diffusirradiance) + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180 ) / (1 + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ; |
929 |
> |
value = ( (diffusirradiance + directirradiance)/(diffusirradiance) + 1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180 ) / (1 + 1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) ; |
930 |
|
|
931 |
|
return(value); |
932 |
|
} |
950 |
|
double value; |
951 |
|
|
952 |
|
value = diffus_irradiance_from_sky_brightness(); |
953 |
< |
value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ); |
953 |
> |
value = value * ( (skyclearness-1) * (1+1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) ); |
954 |
|
|
955 |
|
return(value); |
956 |
|
} |
1013 |
|
{ |
1014 |
|
/*printf("file %s open\n", filename);*/ |
1015 |
|
} |
1016 |
+ |
|
1017 |
+ |
skip_comments(fcoeff_perez); |
1018 |
|
|
1002 |
– |
rewind(fcoeff_perez); /* on se place en debut de fichier */ |
1003 |
– |
|
1019 |
|
for (i=0;i<8;i++) |
1020 |
|
for (j=0;j<20;j++) |
1021 |
|
{ |
1148 |
|
/* degrees into radians */ |
1149 |
|
double radians(double degres) |
1150 |
|
{ |
1151 |
< |
return degres*M_PI/180.0; |
1151 |
> |
return degres*PI/180.0; |
1152 |
|
} |
1153 |
|
|
1154 |
|
/* radian into degrees */ |
1155 |
|
double degres(double radians) |
1156 |
|
{ |
1157 |
< |
return radians/M_PI*180.0; |
1157 |
> |
return radians/PI*180.0; |
1158 |
|
} |
1159 |
|
|
1160 |
|
/* calculation of the angles dzeta and gamma */ |
1198 |
|
fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename); |
1199 |
|
exit(1); |
1200 |
|
} |
1201 |
+ |
|
1202 |
+ |
skip_comments(file_in); |
1203 |
|
|
1187 |
– |
rewind(file_in); |
1188 |
– |
|
1204 |
|
if ( (ptr = malloc(145*sizeof(float))) == NULL ) |
1205 |
|
{ |
1206 |
|
fprintf(stderr,"Out of memory in function theta_ordered\n"); |
1244 |
|
fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename); |
1245 |
|
exit(1); |
1246 |
|
} |
1247 |
+ |
|
1248 |
+ |
skip_comments(file_in); |
1249 |
|
|
1233 |
– |
rewind(file_in); |
1234 |
– |
|
1250 |
|
if ( (ptr = malloc(145*sizeof(float))) == NULL ) |
1251 |
|
{ |
1252 |
|
fprintf(stderr,"Out of memory in function phi_ordered"); |
1288 |
|
for (i=0;i<145;i++) |
1289 |
|
buffer += (*(lv+i))*cos(radians(*(theta+i))); |
1290 |
|
|
1291 |
< |
return buffer*2*M_PI/144; |
1291 |
> |
return buffer*2*PI/144; |
1292 |
|
|
1293 |
|
} |
1294 |
|
|
1304 |
|
double day_angle; |
1305 |
|
double E0; |
1306 |
|
|
1307 |
< |
day_angle = 2*M_PI*(daynumber - 1)/365; |
1307 |
> |
day_angle = 2*PI*(daynumber - 1)/365; |
1308 |
|
E0 = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+ |
1309 |
|
0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle); |
1310 |
|
|
1324 |
|
exit(1); |
1325 |
|
} |
1326 |
|
|
1327 |
< |
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1327 |
> |
m = 1/( cos(sunzenith*PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1328 |
|
return(m); |
1329 |
|
} |
1330 |
|
|
1339 |
|
if (sun_zenith == 0) |
1340 |
|
puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan"); |
1341 |
|
|
1342 |
< |
angle = acos( cos(sun_zenith*M_PI/180)*cos(direction_zenith*M_PI/180) + sin(sun_zenith*M_PI/180)*sin(direction_zenith*M_PI/180)*cos((sun_azimut-direction_azimut)*M_PI/180) ); |
1343 |
< |
angle = angle*180/M_PI; |
1342 |
> |
angle = acos( cos(sun_zenith*PI/180)*cos(direction_zenith*PI/180) + sin(sun_zenith*PI/180)*sin(direction_zenith*PI/180)*cos((sun_azimut-direction_azimut)*PI/180) ); |
1343 |
> |
angle = angle*180/PI; |
1344 |
|
return(angle); |
1345 |
|
} |
1346 |
|
|