ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.8
Committed: Wed Oct 5 17:33:36 2011 UTC (12 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R1
Changes since 2.7: +2 -2 lines
Log Message:
Forgot to include rtmath.h

File Contents

# Content
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