ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.2
Committed: Mon Jun 15 22:27:21 2009 UTC (14 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +24 -12 lines
Log Message:
Added missing library data files used by gendaylit

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.2 static const char RCSid[] = "$Id: gendaylit.c,v 2.1 2009/06/06 20:18:32 greg Exp $";
3 greg 2.1 #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 greg 2.2 #include <ctype.h>
25 greg 2.1
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 greg 2.2 void skip_comments(FILE *fp);
76 greg 2.1
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     fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
476    
477     if (skyclearness==1)
478     normfactor = 0.777778;
479    
480     if (skyclearness>=6)
481     {
482     F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
483     normfactor = normsc()/F2/M_PI;
484     }
485    
486     if ( (skyclearness>1) && (skyclearness<6) )
487     {
488     S_INTER=1;
489     F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
490     normfactor = normsc()/F2/M_PI;
491     }
492    
493     groundbr = zenithbr*normfactor;
494     printf("# Ground ambient level: %.1f\n", groundbr);
495    
496     if (dosun&&(skyclearness>1))
497     groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];
498    
499     groundbr *= gprefl;
500    
501    
502    
503     return;
504     }
505    
506    
507    
508    
509    
510    
511    
512     printsky() /* print out sky */
513     {
514     if (dosun&&(skyclearness>1))
515     {
516     printf("\nvoid light solar\n");
517     printf("0\n0\n");
518     printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
519     printf("\nsolar source sun\n");
520     printf("0\n0\n");
521     printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
522     }
523    
524     if (dosun&&(skyclearness==1))
525     {
526     printf("\nvoid light solar\n");
527     printf("0\n0\n");
528     printf("3 0.0 0.0 0.0\n");
529     printf("\nsolar source sun\n");
530     printf("0\n0\n");
531     printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
532     }
533    
534    
535     printf("\nvoid brightfunc skyfunc\n");
536     printf("2 skybright perezlum.cal\n");
537     printf("0\n");
538     printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
539     *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
540     sundir[0], sundir[1], sundir[2]);
541     }
542    
543    
544     printdefaults() /* print default values */
545     {
546     printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
547     if (zenithbr > 0.0)
548     printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
549     else
550     printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity);
551     printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI));
552     printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI));
553     printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI));
554     }
555    
556    
557     userror(msg) /* print usage error and quit */
558     char *msg;
559     {
560     if (msg != NULL)
561     fprintf(stderr, "%s: Use error - %s\n", progname, msg);
562     fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
563     fprintf(stderr, "or : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
564     fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n");
565     fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
566     fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
567     fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
568     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");
569     exit(1);
570     }
571    
572    
573    
574     double
575     normsc() /* compute normalization factor (E0*F2/L0) */
576     {
577     static double nfc[2][5] = {
578     /* clear sky approx. */
579     {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
580     /* intermediate sky approx. */
581     {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
582     };
583     register double *nf;
584     double x, nsc;
585     register int i;
586     /* polynomial approximation */
587     nf = nfc[S_INTER];
588     x = (altitude - M_PI/4.0)/(M_PI/4.0);
589     nsc = nf[i=4];
590     while (i--)
591     nsc = nsc*x + nf[i];
592    
593     return(nsc);
594     }
595    
596    
597    
598     printhead(ac, av) /* print command header */
599     register int ac;
600     register char **av;
601     {
602     putchar('#');
603     while (ac--) {
604     putchar(' ');
605     fputs(*av++, stdout);
606     }
607     putchar('\n');
608     }
609    
610    
611    
612    
613 greg 2.2 void
614     skip_comments(FILE *fp) /* skip comments in file */
615     {
616     int c;
617    
618     while ((c = getc(fp)) != EOF)
619     if (c == '#') {
620     while ((c = getc(fp)) != EOF)
621     if (c == '\n')
622     break;
623     } else if (!isspace(c)) {
624     ungetc(c, fp);
625     break;
626     }
627     }
628 greg 2.1
629    
630    
631     /* Perez models */
632    
633     /* Perez global horizontal luminous efficacy model */
634     double glob_h_effi_PEREZ()
635     {
636    
637     double value;
638     double category_bounds[10], a[10], b[10], c[10], d[10];
639     int category_total_number, category_number, i;
640    
641    
642     if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
643     fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
644    
645     /* initialize category bounds (clearness index bounds) */
646    
647     category_total_number = 8;
648    
649     category_bounds[1] = 1;
650     category_bounds[2] = 1.065;
651     category_bounds[3] = 1.230;
652     category_bounds[4] = 1.500;
653     category_bounds[5] = 1.950;
654     category_bounds[6] = 2.800;
655     category_bounds[7] = 4.500;
656     category_bounds[8] = 6.200;
657     category_bounds[9] = 12.01;
658    
659    
660     /* initialize model coefficients */
661     a[1] = 96.63;
662     a[2] = 107.54;
663     a[3] = 98.73;
664     a[4] = 92.72;
665     a[5] = 86.73;
666     a[6] = 88.34;
667     a[7] = 78.63;
668     a[8] = 99.65;
669    
670     b[1] = -0.47;
671     b[2] = 0.79;
672     b[3] = 0.70;
673     b[4] = 0.56;
674     b[5] = 0.98;
675     b[6] = 1.39;
676     b[7] = 1.47;
677     b[8] = 1.86;
678    
679     c[1] = 11.50;
680     c[2] = 1.79;
681     c[3] = 4.40;
682     c[4] = 8.36;
683     c[5] = 7.10;
684     c[6] = 6.06;
685     c[7] = 4.93;
686     c[8] = -4.46;
687    
688     d[1] = -9.16;
689     d[2] = -1.19;
690     d[3] = -6.95;
691     d[4] = -8.31;
692     d[5] = -10.94;
693     d[6] = -7.60;
694     d[7] = -11.37;
695     d[8] = -3.15;
696    
697    
698    
699    
700     for (i=1; i<=category_total_number; i++)
701     {
702     if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
703     category_number = i;
704     }
705    
706     value = a[category_number] + b[category_number]*atm_preci_water +
707     c[category_number]*cos(sunzenith*M_PI/180) + d[category_number]*log(skybrightness);
708    
709     return(value);
710     }
711    
712    
713     /* global horizontal diffuse efficacy model, according to PEREZ */
714     double glob_h_diffuse_effi_PEREZ()
715     {
716     double value;
717     double category_bounds[10], a[10], b[10], c[10], d[10];
718     int category_total_number, category_number, i;
719    
720    
721     if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
722     fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
723    
724     /* initialize category bounds (clearness index bounds) */
725    
726     category_total_number = 8;
727    
728     category_bounds[1] = 1;
729     category_bounds[2] = 1.065;
730     category_bounds[3] = 1.230;
731     category_bounds[4] = 1.500;
732     category_bounds[5] = 1.950;
733     category_bounds[6] = 2.800;
734     category_bounds[7] = 4.500;
735     category_bounds[8] = 6.200;
736     category_bounds[9] = 12.01;
737    
738    
739     /* initialize model coefficients */
740     a[1] = 97.24;
741     a[2] = 107.22;
742     a[3] = 104.97;
743     a[4] = 102.39;
744     a[5] = 100.71;
745     a[6] = 106.42;
746     a[7] = 141.88;
747     a[8] = 152.23;
748    
749     b[1] = -0.46;
750     b[2] = 1.15;
751     b[3] = 2.96;
752     b[4] = 5.59;
753     b[5] = 5.94;
754     b[6] = 3.83;
755     b[7] = 1.90;
756     b[8] = 0.35;
757    
758     c[1] = 12.00;
759     c[2] = 0.59;
760     c[3] = -5.53;
761     c[4] = -13.95;
762     c[5] = -22.75;
763     c[6] = -36.15;
764     c[7] = -53.24;
765     c[8] = -45.27;
766    
767     d[1] = -8.91;
768     d[2] = -3.95;
769     d[3] = -8.77;
770     d[4] = -13.90;
771     d[5] = -23.74;
772     d[6] = -28.83;
773     d[7] = -14.03;
774     d[8] = -7.98;
775    
776    
777    
778    
779     for (i=1; i<=category_total_number; i++)
780     {
781     if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
782     category_number = i;
783     }
784    
785     value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) +
786     d[category_number]*log(skybrightness);
787    
788     return(value);
789     }
790    
791    
792     /* direct normal efficacy model, according to PEREZ */
793    
794     double direct_n_effi_PEREZ()
795    
796     {
797     double value;
798     double category_bounds[10], a[10], b[10], c[10], d[10];
799     int category_total_number, category_number, i;
800    
801    
802     if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
803     fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
804    
805    
806     /* initialize category bounds (clearness index bounds) */
807    
808     category_total_number = 8;
809    
810     category_bounds[1] = 1;
811     category_bounds[2] = 1.065;
812     category_bounds[3] = 1.230;
813     category_bounds[4] = 1.500;
814     category_bounds[5] = 1.950;
815     category_bounds[6] = 2.800;
816     category_bounds[7] = 4.500;
817     category_bounds[8] = 6.200;
818     category_bounds[9] = 12.1;
819    
820    
821     /* initialize model coefficients */
822     a[1] = 57.20;
823     a[2] = 98.99;
824     a[3] = 109.83;
825     a[4] = 110.34;
826     a[5] = 106.36;
827     a[6] = 107.19;
828     a[7] = 105.75;
829     a[8] = 101.18;
830    
831     b[1] = -4.55;
832     b[2] = -3.46;
833     b[3] = -4.90;
834     b[4] = -5.84;
835     b[5] = -3.97;
836     b[6] = -1.25;
837     b[7] = 0.77;
838     b[8] = 1.58;
839    
840     c[1] = -2.98;
841     c[2] = -1.21;
842     c[3] = -1.71;
843     c[4] = -1.99;
844     c[5] = -1.75;
845     c[6] = -1.51;
846     c[7] = -1.26;
847     c[8] = -1.10;
848    
849     d[1] = 117.12;
850     d[2] = 12.38;
851     d[3] = -8.81;
852     d[4] = -4.56;
853     d[5] = -6.16;
854     d[6] = -26.73;
855     d[7] = -34.44;
856     d[8] = -8.29;
857    
858    
859    
860     for (i=1; i<=category_total_number; i++)
861     {
862     if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
863     category_number = i;
864     }
865    
866     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;
867    
868     if (value < 0) value = 0;
869    
870     return(value);
871     }
872    
873    
874     /*check the range of epsilon and delta indexes of the perez parametrization*/
875     void check_parametrization()
876     {
877     if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
878     {
879     fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
880     exit(1);
881     }
882     else return;
883     }
884    
885    
886     /* likelihood of the direct and diffuse components */
887     void check_illuminances()
888     {
889     if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
890     {
891     fprintf(stderr,"direct or diffuse illuminances out of range\n");
892     exit(1);
893     }
894     return;
895     }
896    
897    
898     void check_irradiances()
899     {
900     if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
901     {
902     fprintf(stderr,"direct or diffuse irradiances out of range\n");
903     exit(1);
904     }
905     return;
906     }
907    
908    
909    
910     /* Perez sky's brightness */
911     double sky_brightness()
912     {
913     double value;
914    
915     value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
916    
917     return(value);
918     }
919    
920    
921     /* Perez sky's clearness */
922     double sky_clearness()
923     {
924     double value;
925    
926     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) ;
927    
928     return(value);
929     }
930    
931    
932    
933     /* diffus horizontal irradiance from Perez sky's brightness */
934     double diffus_irradiance_from_sky_brightness()
935     {
936     double value;
937    
938     value = skybrightness / air_mass() * ( solar_constant_e*get_eccentricity());
939    
940     return(value);
941     }
942    
943    
944     /* direct normal irradiance from Perez sky's clearness */
945     double direct_irradiance_from_sky_clearness()
946     {
947     double value;
948    
949     value = diffus_irradiance_from_sky_brightness();
950     value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
951    
952     return(value);
953     }
954    
955    
956     void illu_to_irra_index(void)
957     {
958     double test1=0.1, test2=0.1;
959     int counter=0;
960    
961     diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
962     directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
963     skyclearness = sky_clearness();
964     skybrightness = sky_brightness();
965     if (skyclearness>12) skyclearness=12;
966     if (skybrightness<0.05) skybrightness=0.01;
967    
968    
969     while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
970     || skyclearness>skyclearinf || skyclearness<skyclearsup
971     || skybrightness>skybriginf || skybrightness<skybrigsup )
972     && !(counter==5) )
973     {
974     /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
975    
976     test1=diffusirradiance;
977     test2=directirradiance;
978     counter++;
979    
980     diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
981     directirradiance = directilluminance/direct_n_effi_PEREZ();
982     /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
983    
984     skybrightness = sky_brightness();
985     skyclearness = sky_clearness();
986     if (skyclearness>12) skyclearness=12;
987     if (skybrightness<0.05) skybrightness=0.01;
988    
989     /*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
990    
991     }
992    
993    
994     return;
995     }
996    
997    
998     int lect_coeff_perez(char *filename,float **coeff_perez)
999     {
1000     FILE *fcoeff_perez;
1001     float temp;
1002     int i,j;
1003    
1004     if ((fcoeff_perez = frlibopen(filename)) == NULL)
1005     {
1006     fprintf(stderr,"file %s cannot be opened\n", filename);
1007     return 1; /* il y a un probleme de fichier */
1008     }
1009     else
1010     {
1011     /*printf("file %s open\n", filename);*/
1012     }
1013 greg 2.2
1014     skip_comments(fcoeff_perez);
1015 greg 2.1
1016     for (i=0;i<8;i++)
1017     for (j=0;j<20;j++)
1018     {
1019     fscanf(fcoeff_perez,"%f",&temp);
1020     *(*coeff_perez+i*20+j) = temp;
1021     }
1022     fclose(fcoeff_perez);
1023    
1024     return 0; /* tout est OK */
1025     }
1026    
1027    
1028    
1029     /* sky luminance perez model */
1030     double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1031     double epsilon,double Delta,float *coeff_perez)
1032     {
1033     float x[5][4];
1034     int i,j,num_lin;
1035     double c_perez[5];
1036    
1037     if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) )
1038     {
1039     fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1040     exit(1);
1041     }
1042    
1043     /* correction de modele de Perez solar energy ...*/
1044     if ( (epsilon > 1.065) && (epsilon < 2.8) )
1045     {
1046     if ( Delta < 0.2 ) Delta = 0.2;
1047     }
1048    
1049     if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1050     if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1051     if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1052     if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1053     if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1054     if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1055     if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1056     if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1057    
1058     for (i=0;i<5;i++)
1059     for (j=0;j<4;j++)
1060     {
1061     x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1062     /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1063     }
1064    
1065    
1066     if (num_lin)
1067     {
1068     for (i=0;i<5;i++)
1069     c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
1070     }
1071     else
1072     {
1073     c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
1074     c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
1075     c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
1076     c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
1077     c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
1078     }
1079    
1080    
1081     return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) *
1082     (1 + c_perez[2]*exp(c_perez[3]*gamma) +
1083     c_perez[4]*cos(gamma)*cos(gamma) );
1084     }
1085    
1086    
1087    
1088     /* coefficients for the sky luminance perez model */
1089     void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1090     {
1091     float x[5][4];
1092     int i,j,num_lin;
1093    
1094     if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) )
1095     {
1096     fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1097     exit(1);
1098     }
1099    
1100     /* correction du modele de Perez solar energy ...*/
1101     if ( (epsilon > 1.065) && (epsilon < 2.8) )
1102     {
1103     if ( Delta < 0.2 ) Delta = 0.2;
1104     }
1105    
1106     if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1107     if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1108     if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
1109     if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
1110     if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
1111     if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
1112     if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
1113     if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
1114    
1115     for (i=0;i<5;i++)
1116     for (j=0;j<4;j++)
1117     {
1118     x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1119     /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1120     }
1121    
1122    
1123     if (num_lin)
1124     {
1125     for (i=0;i<5;i++)
1126     *(c_perez+i) = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
1127    
1128     }
1129     else
1130     {
1131     *(c_perez+0) = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
1132     *(c_perez+1) = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
1133     *(c_perez+4) = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
1134     *(c_perez+2) = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
1135     *(c_perez+3) = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
1136    
1137    
1138     }
1139    
1140    
1141     return;
1142     }
1143    
1144    
1145     /* degrees into radians */
1146     double radians(double degres)
1147     {
1148     return degres*M_PI/180.0;
1149     }
1150    
1151     /* radian into degrees */
1152     double degres(double radians)
1153     {
1154     return radians/M_PI*180.0;
1155     }
1156    
1157     /* calculation of the angles dzeta and gamma */
1158     void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1159     {
1160     *dzeta = theta; /* dzeta = phi */
1161     if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) )
1162     *gamma = 0;
1163     else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
1164     {
1165     printf("error in calculation of gamma (angle between point and sun");
1166     exit(3);
1167     }
1168     else
1169     *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
1170     }
1171    
1172    
1173     /********************************************************************************/
1174     /* Fonction: theta_ordered */
1175     /* */
1176     /* In: char *filename */
1177     /* */
1178     /* Out: float * */
1179     /* */
1180     /* Update: 29/08/93 */
1181     /* */
1182     /* Rem: theta en degres */
1183     /* */
1184     /* But: fournit les valeurs de theta du fichier d'entree a la memoire */
1185     /* */
1186     /********************************************************************************/
1187     float *theta_ordered(char *filename)
1188     {
1189     int i;
1190     float buffer,*ptr;
1191     FILE *file_in;
1192    
1193     if ( (file_in = frlibopen(filename)) == NULL )
1194     {
1195     fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
1196     exit(1);
1197     }
1198 greg 2.2
1199     skip_comments(file_in);
1200 greg 2.1
1201     if ( (ptr = malloc(145*sizeof(float))) == NULL )
1202     {
1203     fprintf(stderr,"Out of memory in function theta_ordered\n");
1204     exit(1);
1205     }
1206    
1207     for (i=0;i<145;i++)
1208     {
1209     fscanf(file_in,"%f",&buffer);
1210     *(ptr+i) = buffer;
1211     fscanf(file_in,"%f",&buffer);
1212     }
1213    
1214     fclose(file_in);
1215     return ptr;
1216     }
1217    
1218    
1219     /********************************************************************************/
1220     /* Fonction: phi_ordered */
1221     /* */
1222     /* In: char *filename */
1223     /* */
1224     /* Out: float * */
1225     /* */
1226     /* Update: 29/08/93 */
1227     /* */
1228     /* Rem: valeurs de Phi en DEGRES */
1229     /* */
1230     /* But: mettre les angles contenus dans le fichier d'entree dans la memoire */
1231     /* */
1232     /********************************************************************************/
1233     float *phi_ordered(char *filename)
1234     {
1235     int i;
1236     float buffer,*ptr;
1237     FILE *file_in;
1238    
1239     if ( (file_in = frlibopen(filename)) == NULL )
1240     {
1241     fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
1242     exit(1);
1243     }
1244 greg 2.2
1245     skip_comments(file_in);
1246 greg 2.1
1247     if ( (ptr = malloc(145*sizeof(float))) == NULL )
1248     {
1249     fprintf(stderr,"Out of memory in function phi_ordered");
1250     exit(1);
1251     }
1252    
1253     for (i=0;i<145;i++)
1254     {
1255     fscanf(file_in,"%f",&buffer);
1256     fscanf(file_in,"%f",&buffer);
1257     *(ptr+i) = buffer;
1258     }
1259    
1260     fclose(file_in);
1261     return ptr;
1262     }
1263    
1264    
1265     /********************************************************************************/
1266     /* Fonction: integ_lv */
1267     /* */
1268     /* In: float *lv,*theta */
1269     /* int sun_pos */
1270     /* */
1271     /* Out: double */
1272     /* */
1273     /* Update: 29/08/93 */
1274     /* */
1275     /* Rem: */
1276     /* */
1277     /* But: calcul l'integrale de luminance relative sans la dir. du soleil */
1278     /* */
1279     /********************************************************************************/
1280     double integ_lv(float *lv,float *theta)
1281     {
1282     int i;
1283     double buffer=0.0;
1284    
1285     for (i=0;i<145;i++)
1286     buffer += (*(lv+i))*cos(radians(*(theta+i)));
1287    
1288     return buffer*2*M_PI/144;
1289    
1290     }
1291    
1292    
1293    
1294    
1295    
1296    
1297     /* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */
1298    
1299     double get_eccentricity()
1300     {
1301     double day_angle;
1302     double E0;
1303    
1304     day_angle = 2*M_PI*(daynumber - 1)/365;
1305     E0 = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
1306     0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1307    
1308     return (E0);
1309    
1310     }
1311    
1312    
1313     /* enter sunzenith angle (degrees) return relative air mass (double) */
1314     double air_mass()
1315     {
1316     double m;
1317    
1318     if (sunzenith>90)
1319     {
1320     fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1321     exit(1);
1322     }
1323    
1324     m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1325     return(m);
1326     }
1327    
1328    
1329     double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
1330    
1331     {
1332    
1333     double angle;
1334    
1335    
1336     if (sun_zenith == 0)
1337     puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
1338    
1339     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) );
1340     angle = angle*180/M_PI;
1341     return(angle);
1342     }
1343    
1344    
1345    
1346    
1347    
1348    
1349    
1350