1 |
#ifndef lint |
2 |
static const char RCSid[] = "$Id: gendaylit.c,v 2.7 2011/10/05 17:20:55 greg Exp $"; |
3 |
#endif |
4 |
/* Copyright (c) 1994 *Fraunhofer Institut for Solar Energy Systems |
5 |
* Oltmannstr 5, D-79100 Freiburg, Germany |
6 |
* *Agence de l'Environnement et de la Maitrise de l'Energie |
7 |
* Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France |
8 |
* *BOUYGUES |
9 |
* 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France |
10 |
*/ |
11 |
|
12 |
|
13 |
|
14 |
/* |
15 |
* gendaylit.c program to generate the angular distribution of the daylight. |
16 |
* Our zenith is along the Z-axis, the X-axis |
17 |
* points east, and the Y-axis points north. |
18 |
*/ |
19 |
|
20 |
#include <stdio.h> |
21 |
#include <string.h> |
22 |
#include <math.h> |
23 |
#include <stdlib.h> |
24 |
#include <ctype.h> |
25 |
|
26 |
#include "rtmath.h" |
27 |
#include "rtio.h" |
28 |
#include "color.h" |
29 |
#include "paths.h" |
30 |
|
31 |
extern int jdate(int month, int day); |
32 |
extern double stadj(int jd); |
33 |
extern double sdec(int jd); |
34 |
extern double salt(double sd, double st); |
35 |
extern double sazi(double sd, double st); |
36 |
|
37 |
double normsc(); |
38 |
|
39 |
#define DATFILE "coeff_perez.dat" |
40 |
|
41 |
|
42 |
|
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(); |
50 |
double direct_irradiance_from_sky_clearness(); |
51 |
|
52 |
|
53 |
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
54 |
double glob_h_effi_PEREZ(); |
55 |
double glob_h_diffuse_effi_PEREZ(); |
56 |
double direct_n_effi_PEREZ(); |
57 |
/*likelihood check of the epsilon, delta, direct and diffuse components*/ |
58 |
void check_parametrization(); |
59 |
void check_irradiances(); |
60 |
void check_illuminances(); |
61 |
void illu_to_irra_index(); |
62 |
|
63 |
|
64 |
/* Perez sky luminance model */ |
65 |
int lect_coeff_perez(char *filename,float **coeff_perez); |
66 |
double calc_rel_lum_perez(double dzeta,double gamma,double Z, |
67 |
double epsilon,double Delta,float *coeff_perez); |
68 |
/* coefficients for the sky luminance perez model */ |
69 |
void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez); |
70 |
double radians(double degres); |
71 |
double degres(double radians); |
72 |
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z); |
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 |
|
80 |
/* astronomy and geometry*/ |
81 |
double get_eccentricity(); |
82 |
double air_mass(); |
83 |
double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut); |
84 |
|
85 |
|
86 |
/* date*/ |
87 |
int jdate(int month, int day); |
88 |
|
89 |
|
90 |
|
91 |
|
92 |
|
93 |
/* sun calculation constants */ |
94 |
extern double s_latitude; |
95 |
extern double s_longitude; |
96 |
extern double s_meridian; |
97 |
|
98 |
const double AU = 149597890E3; |
99 |
const double solar_constant_e = 1367; /* solar constant W/m^2 */ |
100 |
const double solar_constant_l = 127.5; /* solar constant klux */ |
101 |
|
102 |
const double half_sun_angle = 0.2665; |
103 |
const double half_direct_angle = 2.85; |
104 |
|
105 |
const double skyclearinf = 1.000; /* limitations for the variation of the Perez parameters */ |
106 |
const double skyclearsup = 12.1; |
107 |
const double skybriginf = 0.01; |
108 |
const double skybrigsup = 0.6; |
109 |
|
110 |
|
111 |
|
112 |
/* required values */ |
113 |
int month, day; /* date */ |
114 |
double hour; /* time */ |
115 |
int tsolar; /* 0=standard, 1=solar */ |
116 |
double altitude, azimuth; /* or solar angles */ |
117 |
|
118 |
|
119 |
|
120 |
/* definition of the sky conditions through the Perez parametrization */ |
121 |
double skyclearness, skybrightness; |
122 |
double solarradiance; /*radiance of the sun disk and of the circumsolar area*/ |
123 |
double diffusilluminance, directilluminance, diffusirradiance, directirradiance; |
124 |
double sunzenith, daynumber=150, atm_preci_water=2; |
125 |
|
126 |
double diffnormalization, dirnormalization; |
127 |
double *c_perez; |
128 |
|
129 |
int output=0; /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/ |
130 |
int input=0; /*define the input for the calulation*/ |
131 |
|
132 |
/* default values */ |
133 |
int cloudy = 0; /* 1=standard, 2=uniform */ |
134 |
int dosun = 1; |
135 |
double zenithbr = -1.0; |
136 |
double betaturbidity = 0.1; |
137 |
double gprefl = 0.2; |
138 |
int S_INTER=0; |
139 |
|
140 |
/* computed values */ |
141 |
double sundir[3]; |
142 |
double groundbr; |
143 |
double F2; |
144 |
double solarbr = 0.0; |
145 |
int u_solar = 0; /* -1=irradiance, 1=radiance */ |
146 |
|
147 |
char *progname; |
148 |
char errmsg[128]; |
149 |
|
150 |
|
151 |
main(argc, argv) |
152 |
int argc; |
153 |
char *argv[]; |
154 |
{ |
155 |
int i; |
156 |
|
157 |
progname = argv[0]; |
158 |
if (argc == 2 && !strcmp(argv[1], "-defaults")) { |
159 |
printdefaults(); |
160 |
exit(0); |
161 |
} |
162 |
if (argc < 4) |
163 |
userror("arg count"); |
164 |
if (!strcmp(argv[1], "-ang")) { |
165 |
altitude = atof(argv[2]) * (PI/180); |
166 |
azimuth = atof(argv[3]) * (PI/180); |
167 |
month = 0; |
168 |
} else { |
169 |
month = atoi(argv[1]); |
170 |
if (month < 1 || month > 12) |
171 |
userror("bad month"); |
172 |
day = atoi(argv[2]); |
173 |
if (day < 1 || day > 31) |
174 |
userror("bad day"); |
175 |
hour = atof(argv[3]); |
176 |
if (hour < 0 || hour >= 24) |
177 |
userror("bad hour"); |
178 |
tsolar = argv[3][0] == '+'; |
179 |
} |
180 |
for (i = 4; i < argc; i++) |
181 |
if (argv[i][0] == '-' || argv[i][0] == '+') |
182 |
switch (argv[i][1]) { |
183 |
case 's': |
184 |
cloudy = 0; |
185 |
dosun = argv[i][0] == '+'; |
186 |
break; |
187 |
case 'r': |
188 |
case 'R': |
189 |
u_solar = argv[i][1] == 'R' ? -1 : 1; |
190 |
solarbr = atof(argv[++i]); |
191 |
break; |
192 |
case 'c': |
193 |
cloudy = argv[i][0] == '+' ? 2 : 1; |
194 |
dosun = 0; |
195 |
break; |
196 |
case 't': |
197 |
betaturbidity = atof(argv[++i]); |
198 |
break; |
199 |
case 'b': |
200 |
zenithbr = atof(argv[++i]); |
201 |
break; |
202 |
case 'g': |
203 |
gprefl = atof(argv[++i]); |
204 |
break; |
205 |
case 'a': |
206 |
s_latitude = atof(argv[++i]) * (PI/180); |
207 |
break; |
208 |
case 'o': |
209 |
s_longitude = atof(argv[++i]) * (PI/180); |
210 |
break; |
211 |
case 'm': |
212 |
s_meridian = atof(argv[++i]) * (PI/180); |
213 |
break; |
214 |
|
215 |
|
216 |
case 'O': |
217 |
output = atof(argv[++i]); /*define the unit of the output of the program : |
218 |
sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) |
219 |
default is set to 0*/ |
220 |
break; |
221 |
|
222 |
case 'P': |
223 |
input = 0; /* Perez parameters: epsilon, delta */ |
224 |
skyclearness = atof(argv[++i]); |
225 |
skybrightness = atof(argv[++i]); |
226 |
break; |
227 |
|
228 |
case 'W': /* direct normal Irradiance [W/m^2] */ |
229 |
input = 1; /* diffuse horizontal Irrad. [W/m^2] */ |
230 |
directirradiance = atof(argv[++i]); |
231 |
diffusirradiance = atof(argv[++i]); |
232 |
break; |
233 |
|
234 |
case 'L': /* direct normal Illuminance [Lux] */ |
235 |
input = 2; /* diffuse horizontal Ill. [Lux] */ |
236 |
directilluminance = atof(argv[++i]); |
237 |
diffusilluminance = atof(argv[++i]); |
238 |
break; |
239 |
|
240 |
case 'G': /* direct horizontal Irradiance [W/m^2] */ |
241 |
input = 3; /* diffuse horizontal Irrad. [W/m^2] */ |
242 |
directirradiance = atof(argv[++i]); |
243 |
diffusirradiance = atof(argv[++i]); |
244 |
break; |
245 |
|
246 |
|
247 |
default: |
248 |
sprintf(errmsg, "unknown option: %s", argv[i]); |
249 |
userror(errmsg); |
250 |
} |
251 |
else |
252 |
userror("bad option"); |
253 |
|
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/PI); |
258 |
|
259 |
|
260 |
/* allocation dynamique de memoire pour les pointeurs */ |
261 |
if ( (c_perez = malloc(5*sizeof(double))) == NULL ) |
262 |
{ |
263 |
fprintf(stderr,"Out of memory error in function main !"); |
264 |
exit(1); |
265 |
} |
266 |
|
267 |
|
268 |
printhead(argc, argv); |
269 |
|
270 |
computesky(); |
271 |
printsky(); |
272 |
|
273 |
exit(0); |
274 |
} |
275 |
|
276 |
|
277 |
void |
278 |
computesky() /* compute sky parameters */ |
279 |
{ |
280 |
|
281 |
/* new variables */ |
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; |
287 |
double normfactor; |
288 |
|
289 |
|
290 |
|
291 |
/* compute solar direction */ |
292 |
|
293 |
if (month) { /* from date and time */ |
294 |
int jd; |
295 |
double sd, st; |
296 |
|
297 |
jd = jdate(month, day); /* Julian date */ |
298 |
sd = sdec(jd); /* solar declination */ |
299 |
if (tsolar) /* solar time */ |
300 |
st = hour; |
301 |
else |
302 |
st = hour + stadj(jd); |
303 |
altitude = salt(sd, st); |
304 |
azimuth = sazi(sd, st); |
305 |
|
306 |
daynumber = (double)jdate(month, day); |
307 |
|
308 |
} |
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.*PI/180.; |
316 |
} |
317 |
sundir[0] = -sin(azimuth)*cos(altitude); |
318 |
sundir[1] = -cos(azimuth)*cos(altitude); |
319 |
sundir[2] = sin(altitude); |
320 |
|
321 |
|
322 |
/* calculation for the new functions */ |
323 |
sunzenith = 90 - altitude*180/PI; |
324 |
|
325 |
|
326 |
|
327 |
/* compute the inputs for the calculation of the light distribution over the sky*/ |
328 |
if (input==0) |
329 |
{ |
330 |
check_parametrization(); |
331 |
diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ |
332 |
directirradiance = direct_irradiance_from_sky_clearness(); |
333 |
check_irradiances(); |
334 |
|
335 |
if (output==0 || output==2) |
336 |
{ |
337 |
diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
338 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
339 |
check_illuminances(); |
340 |
} |
341 |
} |
342 |
|
343 |
|
344 |
else if (input==1) |
345 |
{ |
346 |
check_irradiances(); |
347 |
skybrightness = sky_brightness(); |
348 |
skyclearness = sky_clearness(); |
349 |
check_parametrization(); |
350 |
|
351 |
if (output==0 || output==2) |
352 |
{ |
353 |
diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
354 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
355 |
check_illuminances(); |
356 |
} |
357 |
|
358 |
} |
359 |
|
360 |
|
361 |
else if (input==2) |
362 |
{ |
363 |
check_illuminances(); |
364 |
illu_to_irra_index(); |
365 |
check_parametrization(); |
366 |
} |
367 |
|
368 |
|
369 |
else if (input==3) |
370 |
{ |
371 |
if (altitude<=0) |
372 |
{ |
373 |
fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n"); |
374 |
exit(1); |
375 |
} |
376 |
|
377 |
directirradiance=directirradiance/sin(altitude); |
378 |
check_irradiances(); |
379 |
skybrightness = sky_brightness(); |
380 |
skyclearness = sky_clearness(); |
381 |
check_parametrization(); |
382 |
|
383 |
if (output==0 || output==2) |
384 |
{ |
385 |
diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
386 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
387 |
check_illuminances(); |
388 |
} |
389 |
|
390 |
} |
391 |
|
392 |
|
393 |
else {fprintf(stderr,"error in giving the input arguments"); exit(1);} |
394 |
|
395 |
|
396 |
|
397 |
/* normalization factor for the relative sky luminance distribution, diffuse part*/ |
398 |
|
399 |
/* allocation dynamique de memoire pour les pointeurs */ |
400 |
if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL ) |
401 |
{ |
402 |
fprintf(stderr,"Out of memory error in function main !"); |
403 |
exit(1); |
404 |
} |
405 |
|
406 |
/* read the coefficients for the Perez sky luminance model */ |
407 |
if (lect_coeff_perez(DATFILE, &coeff_perez) > 0) |
408 |
{ |
409 |
fprintf(stderr,"lect_coeff_perez does not work\n"); |
410 |
exit(2); |
411 |
} |
412 |
|
413 |
if ( (lv_mod = malloc(145*sizeof(float))) == NULL) |
414 |
{ |
415 |
fprintf(stderr,"Out of memory in function main"); |
416 |
exit(1); |
417 |
} |
418 |
|
419 |
/* read the angles */ |
420 |
theta_o = theta_ordered("defangle.dat"); |
421 |
phi_o = phi_ordered("defangle.dat"); |
422 |
|
423 |
/* parameters for the perez model */ |
424 |
coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez); |
425 |
|
426 |
/*calculation of the modelled luminance */ |
427 |
for (j=0;j<145;j++) |
428 |
{ |
429 |
theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith)); |
430 |
*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
431 |
/*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/ |
432 |
} |
433 |
|
434 |
/* integration of luminance for the normalization factor, diffuse part of the sky*/ |
435 |
diffnormalization = integ_lv(lv_mod, theta_o); |
436 |
/*printf("perez integration %lf\n", diffnormalization);*/ |
437 |
|
438 |
|
439 |
|
440 |
|
441 |
/*normalization coefficient in lumen or in watt*/ |
442 |
if (output==0) |
443 |
{ |
444 |
diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY; |
445 |
} |
446 |
else if (output==1) |
447 |
{ |
448 |
diffnormalization = diffusirradiance/diffnormalization; |
449 |
} |
450 |
else if (output==2) |
451 |
{ |
452 |
diffnormalization = diffusilluminance/diffnormalization; |
453 |
} |
454 |
|
455 |
else {fprintf(stderr,"output argument : wrong number"); exit(1);} |
456 |
|
457 |
|
458 |
|
459 |
|
460 |
/* calculation for the solar source */ |
461 |
if (output==0) |
462 |
solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)))/WHTEFFICACY; |
463 |
|
464 |
else if (output==1) |
465 |
solarradiance = directirradiance/(2*PI*(1-cos(half_sun_angle*PI/180))); |
466 |
|
467 |
else |
468 |
solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180))); |
469 |
|
470 |
|
471 |
|
472 |
|
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 |
*/ |
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*(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(-(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/PI*solarradiance*sundir[2]; |
501 |
|
502 |
groundbr *= gprefl; |
503 |
|
504 |
|
505 |
|
506 |
return; |
507 |
} |
508 |
|
509 |
|
510 |
|
511 |
|
512 |
|
513 |
|
514 |
|
515 |
printsky() /* print out sky */ |
516 |
{ |
517 |
if (dosun&&(skyclearness>1)) |
518 |
{ |
519 |
printf("\nvoid light solar\n"); |
520 |
printf("0\n0\n"); |
521 |
printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance); |
522 |
printf("\nsolar source sun\n"); |
523 |
printf("0\n0\n"); |
524 |
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
525 |
} |
526 |
|
527 |
if (dosun&&(skyclearness==1)) |
528 |
{ |
529 |
printf("\nvoid light solar\n"); |
530 |
printf("0\n0\n"); |
531 |
printf("3 0.0 0.0 0.0\n"); |
532 |
printf("\nsolar source sun\n"); |
533 |
printf("0\n0\n"); |
534 |
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
535 |
} |
536 |
|
537 |
|
538 |
printf("\nvoid brightfunc skyfunc\n"); |
539 |
printf("2 skybright perezlum.cal\n"); |
540 |
printf("0\n"); |
541 |
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
542 |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
543 |
sundir[0], sundir[1], sundir[2]); |
544 |
} |
545 |
|
546 |
|
547 |
printdefaults() /* print default values */ |
548 |
{ |
549 |
printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl); |
550 |
if (zenithbr > 0.0) |
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/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 |
|
560 |
userror(msg) /* print usage error and quit */ |
561 |
char *msg; |
562 |
{ |
563 |
if (msg != NULL) |
564 |
fprintf(stderr, "%s: Use error - %s\n", progname, msg); |
565 |
fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname); |
566 |
fprintf(stderr, "or : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname); |
567 |
fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n"); |
568 |
fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
569 |
fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n"); |
570 |
fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
571 |
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"); |
572 |
exit(1); |
573 |
} |
574 |
|
575 |
|
576 |
|
577 |
double |
578 |
normsc() /* compute normalization factor (E0*F2/L0) */ |
579 |
{ |
580 |
static double nfc[2][5] = { |
581 |
/* clear sky approx. */ |
582 |
{2.766521, 0.547665, -0.369832, 0.009237, 0.059229}, |
583 |
/* intermediate sky approx. */ |
584 |
{3.5556, -2.7152, -1.3081, 1.0660, 0.60227}, |
585 |
}; |
586 |
register double *nf; |
587 |
double x, nsc; |
588 |
register int i; |
589 |
/* polynomial approximation */ |
590 |
nf = nfc[S_INTER]; |
591 |
x = (altitude - PI/4.0)/(PI/4.0); |
592 |
nsc = nf[i=4]; |
593 |
while (i--) |
594 |
nsc = nsc*x + nf[i]; |
595 |
|
596 |
return(nsc); |
597 |
} |
598 |
|
599 |
|
600 |
|
601 |
printhead(ac, av) /* print command header */ |
602 |
register int ac; |
603 |
register char **av; |
604 |
{ |
605 |
putchar('#'); |
606 |
while (ac--) { |
607 |
putchar(' '); |
608 |
fputs(*av++, stdout); |
609 |
} |
610 |
putchar('\n'); |
611 |
} |
612 |
|
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 |
|
634 |
/* Perez models */ |
635 |
|
636 |
/* Perez global horizontal luminous efficacy model */ |
637 |
double glob_h_effi_PEREZ() |
638 |
{ |
639 |
|
640 |
double value; |
641 |
double category_bounds[10], a[10], b[10], c[10], d[10]; |
642 |
int category_total_number, category_number, i; |
643 |
|
644 |
|
645 |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
646 |
fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n"); |
647 |
|
648 |
/* initialize category bounds (clearness index bounds) */ |
649 |
|
650 |
category_total_number = 8; |
651 |
|
652 |
category_bounds[1] = 1; |
653 |
category_bounds[2] = 1.065; |
654 |
category_bounds[3] = 1.230; |
655 |
category_bounds[4] = 1.500; |
656 |
category_bounds[5] = 1.950; |
657 |
category_bounds[6] = 2.800; |
658 |
category_bounds[7] = 4.500; |
659 |
category_bounds[8] = 6.200; |
660 |
category_bounds[9] = 12.01; |
661 |
|
662 |
|
663 |
/* initialize model coefficients */ |
664 |
a[1] = 96.63; |
665 |
a[2] = 107.54; |
666 |
a[3] = 98.73; |
667 |
a[4] = 92.72; |
668 |
a[5] = 86.73; |
669 |
a[6] = 88.34; |
670 |
a[7] = 78.63; |
671 |
a[8] = 99.65; |
672 |
|
673 |
b[1] = -0.47; |
674 |
b[2] = 0.79; |
675 |
b[3] = 0.70; |
676 |
b[4] = 0.56; |
677 |
b[5] = 0.98; |
678 |
b[6] = 1.39; |
679 |
b[7] = 1.47; |
680 |
b[8] = 1.86; |
681 |
|
682 |
c[1] = 11.50; |
683 |
c[2] = 1.79; |
684 |
c[3] = 4.40; |
685 |
c[4] = 8.36; |
686 |
c[5] = 7.10; |
687 |
c[6] = 6.06; |
688 |
c[7] = 4.93; |
689 |
c[8] = -4.46; |
690 |
|
691 |
d[1] = -9.16; |
692 |
d[2] = -1.19; |
693 |
d[3] = -6.95; |
694 |
d[4] = -8.31; |
695 |
d[5] = -10.94; |
696 |
d[6] = -7.60; |
697 |
d[7] = -11.37; |
698 |
d[8] = -3.15; |
699 |
|
700 |
|
701 |
|
702 |
|
703 |
for (i=1; i<=category_total_number; i++) |
704 |
{ |
705 |
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
706 |
category_number = i; |
707 |
} |
708 |
|
709 |
value = a[category_number] + b[category_number]*atm_preci_water + |
710 |
c[category_number]*cos(sunzenith*PI/180) + d[category_number]*log(skybrightness); |
711 |
|
712 |
return(value); |
713 |
} |
714 |
|
715 |
|
716 |
/* global horizontal diffuse efficacy model, according to PEREZ */ |
717 |
double glob_h_diffuse_effi_PEREZ() |
718 |
{ |
719 |
double value; |
720 |
double category_bounds[10], a[10], b[10], c[10], d[10]; |
721 |
int category_total_number, category_number, i; |
722 |
|
723 |
|
724 |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
725 |
fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n"); |
726 |
|
727 |
/* initialize category bounds (clearness index bounds) */ |
728 |
|
729 |
category_total_number = 8; |
730 |
|
731 |
category_bounds[1] = 1; |
732 |
category_bounds[2] = 1.065; |
733 |
category_bounds[3] = 1.230; |
734 |
category_bounds[4] = 1.500; |
735 |
category_bounds[5] = 1.950; |
736 |
category_bounds[6] = 2.800; |
737 |
category_bounds[7] = 4.500; |
738 |
category_bounds[8] = 6.200; |
739 |
category_bounds[9] = 12.01; |
740 |
|
741 |
|
742 |
/* initialize model coefficients */ |
743 |
a[1] = 97.24; |
744 |
a[2] = 107.22; |
745 |
a[3] = 104.97; |
746 |
a[4] = 102.39; |
747 |
a[5] = 100.71; |
748 |
a[6] = 106.42; |
749 |
a[7] = 141.88; |
750 |
a[8] = 152.23; |
751 |
|
752 |
b[1] = -0.46; |
753 |
b[2] = 1.15; |
754 |
b[3] = 2.96; |
755 |
b[4] = 5.59; |
756 |
b[5] = 5.94; |
757 |
b[6] = 3.83; |
758 |
b[7] = 1.90; |
759 |
b[8] = 0.35; |
760 |
|
761 |
c[1] = 12.00; |
762 |
c[2] = 0.59; |
763 |
c[3] = -5.53; |
764 |
c[4] = -13.95; |
765 |
c[5] = -22.75; |
766 |
c[6] = -36.15; |
767 |
c[7] = -53.24; |
768 |
c[8] = -45.27; |
769 |
|
770 |
d[1] = -8.91; |
771 |
d[2] = -3.95; |
772 |
d[3] = -8.77; |
773 |
d[4] = -13.90; |
774 |
d[5] = -23.74; |
775 |
d[6] = -28.83; |
776 |
d[7] = -14.03; |
777 |
d[8] = -7.98; |
778 |
|
779 |
|
780 |
|
781 |
|
782 |
for (i=1; i<=category_total_number; i++) |
783 |
{ |
784 |
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
785 |
category_number = i; |
786 |
} |
787 |
|
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); |
792 |
} |
793 |
|
794 |
|
795 |
/* direct normal efficacy model, according to PEREZ */ |
796 |
|
797 |
double direct_n_effi_PEREZ() |
798 |
|
799 |
{ |
800 |
double value; |
801 |
double category_bounds[10], a[10], b[10], c[10], d[10]; |
802 |
int category_total_number, category_number, i; |
803 |
|
804 |
|
805 |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
806 |
fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n"); |
807 |
|
808 |
|
809 |
/* initialize category bounds (clearness index bounds) */ |
810 |
|
811 |
category_total_number = 8; |
812 |
|
813 |
category_bounds[1] = 1; |
814 |
category_bounds[2] = 1.065; |
815 |
category_bounds[3] = 1.230; |
816 |
category_bounds[4] = 1.500; |
817 |
category_bounds[5] = 1.950; |
818 |
category_bounds[6] = 2.800; |
819 |
category_bounds[7] = 4.500; |
820 |
category_bounds[8] = 6.200; |
821 |
category_bounds[9] = 12.1; |
822 |
|
823 |
|
824 |
/* initialize model coefficients */ |
825 |
a[1] = 57.20; |
826 |
a[2] = 98.99; |
827 |
a[3] = 109.83; |
828 |
a[4] = 110.34; |
829 |
a[5] = 106.36; |
830 |
a[6] = 107.19; |
831 |
a[7] = 105.75; |
832 |
a[8] = 101.18; |
833 |
|
834 |
b[1] = -4.55; |
835 |
b[2] = -3.46; |
836 |
b[3] = -4.90; |
837 |
b[4] = -5.84; |
838 |
b[5] = -3.97; |
839 |
b[6] = -1.25; |
840 |
b[7] = 0.77; |
841 |
b[8] = 1.58; |
842 |
|
843 |
c[1] = -2.98; |
844 |
c[2] = -1.21; |
845 |
c[3] = -1.71; |
846 |
c[4] = -1.99; |
847 |
c[5] = -1.75; |
848 |
c[6] = -1.51; |
849 |
c[7] = -1.26; |
850 |
c[8] = -1.10; |
851 |
|
852 |
d[1] = 117.12; |
853 |
d[2] = 12.38; |
854 |
d[3] = -8.81; |
855 |
d[4] = -4.56; |
856 |
d[5] = -6.16; |
857 |
d[6] = -26.73; |
858 |
d[7] = -34.44; |
859 |
d[8] = -8.29; |
860 |
|
861 |
|
862 |
|
863 |
for (i=1; i<=category_total_number; i++) |
864 |
{ |
865 |
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
866 |
category_number = i; |
867 |
} |
868 |
|
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 |
|
873 |
return(value); |
874 |
} |
875 |
|
876 |
|
877 |
/*check the range of epsilon and delta indexes of the perez parametrization*/ |
878 |
void check_parametrization() |
879 |
{ |
880 |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup) |
881 |
{ |
882 |
fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness); |
883 |
exit(1); |
884 |
} |
885 |
else return; |
886 |
} |
887 |
|
888 |
|
889 |
/* likelihood of the direct and diffuse components */ |
890 |
void check_illuminances() |
891 |
{ |
892 |
if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) )) |
893 |
{ |
894 |
fprintf(stderr,"direct or diffuse illuminances out of range\n"); |
895 |
exit(1); |
896 |
} |
897 |
return; |
898 |
} |
899 |
|
900 |
|
901 |
void check_irradiances() |
902 |
{ |
903 |
if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) )) |
904 |
{ |
905 |
fprintf(stderr,"direct or diffuse irradiances out of range\n"); |
906 |
exit(1); |
907 |
} |
908 |
return; |
909 |
} |
910 |
|
911 |
|
912 |
|
913 |
/* Perez sky's brightness */ |
914 |
double sky_brightness() |
915 |
{ |
916 |
double value; |
917 |
|
918 |
value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity()); |
919 |
|
920 |
return(value); |
921 |
} |
922 |
|
923 |
|
924 |
/* Perez sky's clearness */ |
925 |
double sky_clearness() |
926 |
{ |
927 |
double value; |
928 |
|
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 |
} |
933 |
|
934 |
|
935 |
|
936 |
/* diffus horizontal irradiance from Perez sky's brightness */ |
937 |
double diffus_irradiance_from_sky_brightness() |
938 |
{ |
939 |
double value; |
940 |
|
941 |
value = skybrightness / air_mass() * ( solar_constant_e*get_eccentricity()); |
942 |
|
943 |
return(value); |
944 |
} |
945 |
|
946 |
|
947 |
/* direct normal irradiance from Perez sky's clearness */ |
948 |
double direct_irradiance_from_sky_clearness() |
949 |
{ |
950 |
double value; |
951 |
|
952 |
value = diffus_irradiance_from_sky_brightness(); |
953 |
value = value * ( (skyclearness-1) * (1+1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) ); |
954 |
|
955 |
return(value); |
956 |
} |
957 |
|
958 |
|
959 |
void illu_to_irra_index(void) |
960 |
{ |
961 |
double test1=0.1, test2=0.1; |
962 |
int counter=0; |
963 |
|
964 |
diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000); |
965 |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000); |
966 |
skyclearness = sky_clearness(); |
967 |
skybrightness = sky_brightness(); |
968 |
if (skyclearness>12) skyclearness=12; |
969 |
if (skybrightness<0.05) skybrightness=0.01; |
970 |
|
971 |
|
972 |
while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10) |
973 |
|| skyclearness>skyclearinf || skyclearness<skyclearsup |
974 |
|| skybrightness>skybriginf || skybrightness<skybrigsup ) |
975 |
&& !(counter==5) ) |
976 |
{ |
977 |
/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/ |
978 |
|
979 |
test1=diffusirradiance; |
980 |
test2=directirradiance; |
981 |
counter++; |
982 |
|
983 |
diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ(); |
984 |
directirradiance = directilluminance/direct_n_effi_PEREZ(); |
985 |
/*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/ |
986 |
|
987 |
skybrightness = sky_brightness(); |
988 |
skyclearness = sky_clearness(); |
989 |
if (skyclearness>12) skyclearness=12; |
990 |
if (skybrightness<0.05) skybrightness=0.01; |
991 |
|
992 |
/*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/ |
993 |
|
994 |
} |
995 |
|
996 |
|
997 |
return; |
998 |
} |
999 |
|
1000 |
|
1001 |
int lect_coeff_perez(char *filename,float **coeff_perez) |
1002 |
{ |
1003 |
FILE *fcoeff_perez; |
1004 |
float temp; |
1005 |
int i,j; |
1006 |
|
1007 |
if ((fcoeff_perez = frlibopen(filename)) == NULL) |
1008 |
{ |
1009 |
fprintf(stderr,"file %s cannot be opened\n", filename); |
1010 |
return 1; /* il y a un probleme de fichier */ |
1011 |
} |
1012 |
else |
1013 |
{ |
1014 |
/*printf("file %s open\n", filename);*/ |
1015 |
} |
1016 |
|
1017 |
skip_comments(fcoeff_perez); |
1018 |
|
1019 |
for (i=0;i<8;i++) |
1020 |
for (j=0;j<20;j++) |
1021 |
{ |
1022 |
fscanf(fcoeff_perez,"%f",&temp); |
1023 |
*(*coeff_perez+i*20+j) = temp; |
1024 |
} |
1025 |
fclose(fcoeff_perez); |
1026 |
|
1027 |
return 0; /* tout est OK */ |
1028 |
} |
1029 |
|
1030 |
|
1031 |
|
1032 |
/* sky luminance perez model */ |
1033 |
double calc_rel_lum_perez(double dzeta,double gamma,double Z, |
1034 |
double epsilon,double Delta,float *coeff_perez) |
1035 |
{ |
1036 |
float x[5][4]; |
1037 |
int i,j,num_lin; |
1038 |
double c_perez[5]; |
1039 |
|
1040 |
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1041 |
{ |
1042 |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); |
1043 |
exit(1); |
1044 |
} |
1045 |
|
1046 |
/* correction de modele de Perez solar energy ...*/ |
1047 |
if ( (epsilon > 1.065) && (epsilon < 2.8) ) |
1048 |
{ |
1049 |
if ( Delta < 0.2 ) Delta = 0.2; |
1050 |
} |
1051 |
|
1052 |
if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0; |
1053 |
if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1; |
1054 |
if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2; |
1055 |
if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3; |
1056 |
if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4; |
1057 |
if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5; |
1058 |
if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6; |
1059 |
if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7; |
1060 |
|
1061 |
for (i=0;i<5;i++) |
1062 |
for (j=0;j<4;j++) |
1063 |
{ |
1064 |
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
1065 |
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
1066 |
} |
1067 |
|
1068 |
|
1069 |
if (num_lin) |
1070 |
{ |
1071 |
for (i=0;i<5;i++) |
1072 |
c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z); |
1073 |
} |
1074 |
else |
1075 |
{ |
1076 |
c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z); |
1077 |
c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z); |
1078 |
c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z); |
1079 |
c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3]; |
1080 |
c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3]; |
1081 |
} |
1082 |
|
1083 |
|
1084 |
return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) * |
1085 |
(1 + c_perez[2]*exp(c_perez[3]*gamma) + |
1086 |
c_perez[4]*cos(gamma)*cos(gamma) ); |
1087 |
} |
1088 |
|
1089 |
|
1090 |
|
1091 |
/* coefficients for the sky luminance perez model */ |
1092 |
void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez) |
1093 |
{ |
1094 |
float x[5][4]; |
1095 |
int i,j,num_lin; |
1096 |
|
1097 |
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
1098 |
{ |
1099 |
fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n"); |
1100 |
exit(1); |
1101 |
} |
1102 |
|
1103 |
/* correction du modele de Perez solar energy ...*/ |
1104 |
if ( (epsilon > 1.065) && (epsilon < 2.8) ) |
1105 |
{ |
1106 |
if ( Delta < 0.2 ) Delta = 0.2; |
1107 |
} |
1108 |
|
1109 |
if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0; |
1110 |
if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1; |
1111 |
if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2; |
1112 |
if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3; |
1113 |
if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4; |
1114 |
if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5; |
1115 |
if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6; |
1116 |
if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7; |
1117 |
|
1118 |
for (i=0;i<5;i++) |
1119 |
for (j=0;j<4;j++) |
1120 |
{ |
1121 |
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
1122 |
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
1123 |
} |
1124 |
|
1125 |
|
1126 |
if (num_lin) |
1127 |
{ |
1128 |
for (i=0;i<5;i++) |
1129 |
*(c_perez+i) = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z); |
1130 |
|
1131 |
} |
1132 |
else |
1133 |
{ |
1134 |
*(c_perez+0) = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z); |
1135 |
*(c_perez+1) = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z); |
1136 |
*(c_perez+4) = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z); |
1137 |
*(c_perez+2) = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3]; |
1138 |
*(c_perez+3) = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3]; |
1139 |
|
1140 |
|
1141 |
} |
1142 |
|
1143 |
|
1144 |
return; |
1145 |
} |
1146 |
|
1147 |
|
1148 |
/* degrees into radians */ |
1149 |
double radians(double degres) |
1150 |
{ |
1151 |
return degres*PI/180.0; |
1152 |
} |
1153 |
|
1154 |
/* radian into degrees */ |
1155 |
double degres(double radians) |
1156 |
{ |
1157 |
return radians/PI*180.0; |
1158 |
} |
1159 |
|
1160 |
/* calculation of the angles dzeta and gamma */ |
1161 |
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z) |
1162 |
{ |
1163 |
*dzeta = theta; /* dzeta = phi */ |
1164 |
if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) ) |
1165 |
*gamma = 0; |
1166 |
else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 ) |
1167 |
{ |
1168 |
printf("error in calculation of gamma (angle between point and sun"); |
1169 |
exit(3); |
1170 |
} |
1171 |
else |
1172 |
*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)); |
1173 |
} |
1174 |
|
1175 |
|
1176 |
/********************************************************************************/ |
1177 |
/* Fonction: theta_ordered */ |
1178 |
/* */ |
1179 |
/* In: char *filename */ |
1180 |
/* */ |
1181 |
/* Out: float * */ |
1182 |
/* */ |
1183 |
/* Update: 29/08/93 */ |
1184 |
/* */ |
1185 |
/* Rem: theta en degres */ |
1186 |
/* */ |
1187 |
/* But: fournit les valeurs de theta du fichier d'entree a la memoire */ |
1188 |
/* */ |
1189 |
/********************************************************************************/ |
1190 |
float *theta_ordered(char *filename) |
1191 |
{ |
1192 |
int i; |
1193 |
float buffer,*ptr; |
1194 |
FILE *file_in; |
1195 |
|
1196 |
if ( (file_in = frlibopen(filename)) == NULL ) |
1197 |
{ |
1198 |
fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename); |
1199 |
exit(1); |
1200 |
} |
1201 |
|
1202 |
skip_comments(file_in); |
1203 |
|
1204 |
if ( (ptr = malloc(145*sizeof(float))) == NULL ) |
1205 |
{ |
1206 |
fprintf(stderr,"Out of memory in function theta_ordered\n"); |
1207 |
exit(1); |
1208 |
} |
1209 |
|
1210 |
for (i=0;i<145;i++) |
1211 |
{ |
1212 |
fscanf(file_in,"%f",&buffer); |
1213 |
*(ptr+i) = buffer; |
1214 |
fscanf(file_in,"%f",&buffer); |
1215 |
} |
1216 |
|
1217 |
fclose(file_in); |
1218 |
return ptr; |
1219 |
} |
1220 |
|
1221 |
|
1222 |
/********************************************************************************/ |
1223 |
/* Fonction: phi_ordered */ |
1224 |
/* */ |
1225 |
/* In: char *filename */ |
1226 |
/* */ |
1227 |
/* Out: float * */ |
1228 |
/* */ |
1229 |
/* Update: 29/08/93 */ |
1230 |
/* */ |
1231 |
/* Rem: valeurs de Phi en DEGRES */ |
1232 |
/* */ |
1233 |
/* But: mettre les angles contenus dans le fichier d'entree dans la memoire */ |
1234 |
/* */ |
1235 |
/********************************************************************************/ |
1236 |
float *phi_ordered(char *filename) |
1237 |
{ |
1238 |
int i; |
1239 |
float buffer,*ptr; |
1240 |
FILE *file_in; |
1241 |
|
1242 |
if ( (file_in = frlibopen(filename)) == NULL ) |
1243 |
{ |
1244 |
fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename); |
1245 |
exit(1); |
1246 |
} |
1247 |
|
1248 |
skip_comments(file_in); |
1249 |
|
1250 |
if ( (ptr = malloc(145*sizeof(float))) == NULL ) |
1251 |
{ |
1252 |
fprintf(stderr,"Out of memory in function phi_ordered"); |
1253 |
exit(1); |
1254 |
} |
1255 |
|
1256 |
for (i=0;i<145;i++) |
1257 |
{ |
1258 |
fscanf(file_in,"%f",&buffer); |
1259 |
fscanf(file_in,"%f",&buffer); |
1260 |
*(ptr+i) = buffer; |
1261 |
} |
1262 |
|
1263 |
fclose(file_in); |
1264 |
return ptr; |
1265 |
} |
1266 |
|
1267 |
|
1268 |
/********************************************************************************/ |
1269 |
/* Fonction: integ_lv */ |
1270 |
/* */ |
1271 |
/* In: float *lv,*theta */ |
1272 |
/* int sun_pos */ |
1273 |
/* */ |
1274 |
/* Out: double */ |
1275 |
/* */ |
1276 |
/* Update: 29/08/93 */ |
1277 |
/* */ |
1278 |
/* Rem: */ |
1279 |
/* */ |
1280 |
/* But: calcul l'integrale de luminance relative sans la dir. du soleil */ |
1281 |
/* */ |
1282 |
/********************************************************************************/ |
1283 |
double integ_lv(float *lv,float *theta) |
1284 |
{ |
1285 |
int i; |
1286 |
double buffer=0.0; |
1287 |
|
1288 |
for (i=0;i<145;i++) |
1289 |
buffer += (*(lv+i))*cos(radians(*(theta+i))); |
1290 |
|
1291 |
return buffer*2*PI/144; |
1292 |
|
1293 |
} |
1294 |
|
1295 |
|
1296 |
|
1297 |
|
1298 |
|
1299 |
|
1300 |
/* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */ |
1301 |
|
1302 |
double get_eccentricity() |
1303 |
{ |
1304 |
double day_angle; |
1305 |
double E0; |
1306 |
|
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 |
|
1311 |
return (E0); |
1312 |
|
1313 |
} |
1314 |
|
1315 |
|
1316 |
/* enter sunzenith angle (degrees) return relative air mass (double) */ |
1317 |
double air_mass() |
1318 |
{ |
1319 |
double m; |
1320 |
|
1321 |
if (sunzenith>90) |
1322 |
{ |
1323 |
fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n"); |
1324 |
exit(1); |
1325 |
} |
1326 |
|
1327 |
m = 1/( cos(sunzenith*PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
1328 |
return(m); |
1329 |
} |
1330 |
|
1331 |
|
1332 |
double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut) |
1333 |
|
1334 |
{ |
1335 |
|
1336 |
double angle; |
1337 |
|
1338 |
|
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*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 |
|
1347 |
|
1348 |
|
1349 |
|
1350 |
|
1351 |
|
1352 |
|
1353 |
|