ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.3
Committed: Sat Jun 20 21:34:34 2009 UTC (13 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0
Changes since 2.2: +4 -2 lines
Log Message:
Eliminated debug message

File Contents

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