ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/gen/gendaylit.c
(Generate patch)

Comparing src/gen/gendaylit.c (file contents):
Revision 2.8 by greg, Wed Oct 5 17:33:36 2011 UTC vs.
Revision 2.12 by greg, Fri Aug 9 16:44:19 2013 UTC

# Line 1 | Line 1
1 < #ifndef lint
2 < static const char RCSid[] = "$Id$";
3 < #endif
4 < /*        Copyright (c) 1994    *Fraunhofer Institut for Solar Energy Systems
5 < *                              Oltmannstr 5, D-79100 Freiburg, Germany
1 > /*      Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems
2 > *                              Heidenhofstr. 2, D-79110 Freiburg, Germany
3   *                              *Agence de l'Environnement et de la Maitrise de l'Energie
4   *                              Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France
5   *                              *BOUYGUES
6   *                              1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France
7   */
8  
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 –
9   #include  <stdio.h>
10   #include  <string.h>
11   #include  <math.h>
12   #include  <stdlib.h>
24 – #include  <ctype.h>
13  
26 – #include  "rtmath.h"
27 – #include  "rtio.h"
14   #include  "color.h"
15 + #include  "sun.h"
16   #include  "paths.h"
17  
18 < 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);
18 > #define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
19  
20   double  normsc();
21  
22 < #define DATFILE         "coeff_perez.dat"
22 > /*static        char *rcsid="$Header$";*/
23  
24 + float coeff_perez[] = {
25 +        1.3525,-0.2576,-0.2690,-1.4366,-0.7670,0.0007,1.2734,-0.1233,2.8000,0.6004,1.2375,1.000,1.8734,0.6297,
26 +        0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,-1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156,
27 +        6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285,-1.1000,-0.2515,
28 +        0.8952,0.0156,0.2782,-0.1812,-4.5000,1.1766,24.7219,-13.0812,-37.7000,34.8438,-5.0000,1.5218,3.9229,
29 +        -2.6204,-0.0156,0.1597,0.4199,-0.5562,-0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297,
30 +        33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329,-0.6000,
31 +        -0.3566,-2.5000,2.3250,0.2937,0.0496,-5.6812,1.8415,21.0000,-4.7656,-21.5906,7.2492,-3.5000,-0.1554,
32 +        1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294,-1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750,
33 +        14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844,-1.0000,
34 +        0.0211,0.5025,-0.5119,-0.3000,0.1922,0.7023,-1.6317,19.0000,-5.0000,1.2438,-1.9094,-4.0000,0.0250,0.3844,
35 +        0.2656,1.0468,-0.3788,-2.4517,1.4656,-1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625,
36 +        -14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636};
37  
38  
39 + float defangle_theta[] = {
40 +        84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84,
41 +        84, 84, 84, 84, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72,
42 +        72, 72, 72, 72, 72, 72, 72, 72, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60,
43 +        60, 60, 60, 60, 60, 60, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48,
44 +        48, 48, 48, 48, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 24, 24, 24, 24,
45 +        24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};
46 +
47 + float defangle_phi[] = {
48 +        0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264,
49 +        276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180,
50 +        192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105,
51 +        120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75,
52 +        90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
53 +        80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
54 +        240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
55 +
56 +
57 +
58   /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
59   double sky_brightness();
60   double sky_clearness();
46 – void computesky();
61  
62   /* calculation of the direct and diffuse components from the Perez parametrization */
63 < double  diffus_irradiance_from_sky_brightness();
63 > double  diffuse_irradiance_from_sky_brightness();
64   double  direct_irradiance_from_sky_clearness();
65  
66 + /* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */
67 + /* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
68  
53 – /* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
69   double  glob_h_effi_PEREZ();
70   double  glob_h_diffuse_effi_PEREZ();
71   double  direct_n_effi_PEREZ();
72 +
73   /*likelihood check of the epsilon, delta, direct and diffuse components*/
74   void    check_parametrization();
75   void    check_irradiances();
76   void    check_illuminances();
77   void    illu_to_irra_index();
78 + void    print_error_sky();
79  
80 <
81 < /* 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);
80 > double  calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
81 > void    coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
82   double  radians(double degres);
83   double  degres(double radians);
84   void    theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);
85   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);
86  
87 + void printdefaults();
88 + void check_sun_position();
89 + void computesky();
90 + void printhead(int ac, char** av);
91 + void userror(char* msg);
92 + void printsky();
93  
94 + FILE * frlibopen(char* fname);
95  
96   /* astronomy and geometry*/
97   double  get_eccentricity();
98   double  air_mass();
83 – double  get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
99  
100 <
101 < /* date*/
100 > double  solar_sunset();
101 > double  solar_sunrise();
102 > double  stadj();
103   int     jdate(int month, int day);
104  
105  
106  
91 –
92 –
107   /* sun calculation constants */
108 < extern double  s_latitude;
109 < extern double  s_longitude;
110 < extern double  s_meridian;
108 > extern double   s_latitude;
109 > extern double   s_longitude;
110 > extern double   s_meridian;
111  
112   const double    AU = 149597890E3;
113   const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
# Line 102 | Line 116 | const double   solar_constant_l = 127.5;   /* solar co
116   const double    half_sun_angle = 0.2665;
117   const double    half_direct_angle = 2.85;
118  
119 < const double    skyclearinf = 1.000;    /* limitations for the variation of the Perez parameters */
119 > const double    skyclearinf = 1.0;      /* limitations for the variation of the Perez parameters */
120   const double    skyclearsup = 12.1;
121   const double    skybriginf = 0.01;
122   const double    skybrigsup = 0.6;
# Line 118 | Line 132 | double  altitude, azimuth;                     /* or solar angles */
132  
133  
134   /* definition of the sky conditions through the Perez parametrization */
135 < double  skyclearness, skybrightness;
136 < double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
137 < double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
138 < double  sunzenith, daynumber=150, atm_preci_water=2;
135 > double  skyclearness = 0;
136 > double  skybrightness = 0;
137 > double  solarradiance;
138 > double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
139 > double  sunzenith, daynumber, atm_preci_water=2;
140  
141 < double  diffnormalization, dirnormalization;
141 > /*double  sunaltitude_border = 0;*/
142 > double  diffnormalization = 0;
143 > double  dirnormalization = 0;
144   double  *c_perez;
145  
146 < int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
147 < int     input=0;        /*define the input for the calulation*/
146 > int     output=0;       /* define the unit of the output (sky luminance or radiance): */
147 >                        /* visible watt=0, solar watt=1, lumen=2 */
148 > int     input=0;        /* define the input for the calulation */
149  
150 + int     suppress_warnings=0;
151 +
152          /* default values */
153 < int  cloudy = 0;                                /* 1=standard, 2=uniform */
154 < int  dosun = 1;
153 > int     cloudy = 0;                             /* 1=standard, 2=uniform */
154 > int     dosun = 1;
155   double  zenithbr = -1.0;
156   double  betaturbidity = 0.1;
157   double  gprefl = 0.2;
# Line 139 | Line 159 | int    S_INTER=0;
159  
160          /* computed values */
161   double  sundir[3];
162 < double  groundbr;
162 > double  groundbr = 0;
163   double  F2;
164   double  solarbr = 0.0;
165   int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
166 + float   timeinterval = 0;
167  
168 < char  *progname;
169 < char  errmsg[128];
168 > char    *progname;
169 > char    errmsg[128];
170  
171  
172 < main(argc, argv)
173 < int  argc;
174 < char  *argv[];
172 >
173 >
174 > int main(int argc, char** argv)
175   {
176          int  i;
177  
178          progname = argv[0];
179          if (argc == 2 && !strcmp(argv[1], "-defaults")) {
180                  printdefaults();
181 <                exit(0);
181 >                return 0;
182          }
183          if (argc < 4)
184                  userror("arg count");
185          if (!strcmp(argv[1], "-ang")) {
186 <                altitude = atof(argv[2]) * (PI/180);
187 <                azimuth = atof(argv[3]) * (PI/180);
186 >                altitude = atof(argv[2]) * (M_PI/180);
187 >                azimuth = atof(argv[3]) * (M_PI/180);
188                  month = 0;
189          } else {
190                  month = atoi(argv[1]);
# Line 184 | Line 205 | char  *argv[];
205                                  cloudy = 0;
206                                  dosun = argv[i][0] == '+';
207                                  break;
187 –                        case 'r':
208                          case 'R':
209                                  u_solar = argv[i][1] == 'R' ? -1 : 1;
210                                  solarbr = atof(argv[++i]);
# Line 196 | Line 216 | char  *argv[];
216                          case 't':
217                                  betaturbidity = atof(argv[++i]);
218                                  break;
219 +                        case 'w':
220 +                                suppress_warnings = 1;
221 +                                break;                  
222                          case 'b':
223                                  zenithbr = atof(argv[++i]);
224                                  break;
# Line 203 | Line 226 | char  *argv[];
226                                  gprefl = atof(argv[++i]);
227                                  break;
228                          case 'a':
229 <                                s_latitude = atof(argv[++i]) * (PI/180);
229 >                                s_latitude = atof(argv[++i]) * (M_PI/180);
230                                  break;
231                          case 'o':
232 <                                s_longitude = atof(argv[++i]) * (PI/180);
232 >                                s_longitude = atof(argv[++i]) * (M_PI/180);
233                                  break;
234                          case 'm':
235 <                                s_meridian = atof(argv[++i]) * (PI/180);
235 >                                s_meridian = atof(argv[++i]) * (M_PI/180);
236                                  break;
214 –
237                          
238                          case 'O':
239                                  output = atof(argv[++i]);       /*define the unit of the output of the program :
240 <                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm)
219 <                                                                default is set to 0*/
240 >                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
241                                  break;
242                                  
243                          case 'P':
# Line 228 | Line 249 | char  *argv[];
249                          case 'W':                                       /* direct normal Irradiance [W/m^2] */
250                                  input = 1;                              /* diffuse horizontal Irrad. [W/m^2] */
251                                  directirradiance = atof(argv[++i]);
252 <                                diffusirradiance = atof(argv[++i]);
252 >                                diffuseirradiance = atof(argv[++i]);
253                                  break;
254                                  
255                          case 'L':                                       /* direct normal Illuminance [Lux] */
256                                  input = 2;                              /* diffuse horizontal Ill. [Lux] */
257                                  directilluminance = atof(argv[++i]);
258 <                                diffusilluminance = atof(argv[++i]);
258 >                                diffuseilluminance = atof(argv[++i]);
259                                  break;
260                          
261                          case 'G':                                       /* direct horizontal Irradiance [W/m^2] */
262                                  input = 3;                              /* diffuse horizontal Irrad. [W/m^2] */
263                                  directirradiance = atof(argv[++i]);
264 <                                diffusirradiance = atof(argv[++i]);
264 >                                diffuseirradiance = atof(argv[++i]);
265                                  break;
266 +                        
267 +                        case 'E':                                       /* Erbs model based on the */
268 +                                input = 4;                              /* global-horizontal irradiance [W/m^2] */
269 +                                globalirradiance = atof(argv[++i]);
270 +                                break;
271 +                        
272 +                        /*
273 +                        case 'l':
274 +                                sunaltitude_border = atof(argv[++i]);
275 +                                break;
276 +                        */
277                                  
278 +                        case 'i':
279 +                                timeinterval = atof(argv[++i]);
280 +                                break;
281                          
282 +                        
283                          default:
284                                  sprintf(errmsg, "unknown option: %s", argv[i]);
285                                  userror(errmsg);
# Line 251 | Line 287 | char  *argv[];
287                  else
288                          userror("bad option");
289  
290 <        if (fabs(s_meridian-s_longitude) > 30*PI/180)
290 >        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
291                  fprintf(stderr,
292                      "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
293 <                    progname, (s_longitude-s_meridian)*12/PI);
293 >                    progname, (s_longitude-s_meridian)*12/M_PI);
294  
295  
296 <        /* allocation dynamique de memoire pour les pointeurs */
297 <        if ( (c_perez = malloc(5*sizeof(double))) == NULL )
296 >        /* dynamic memory allocation for the pointers */
297 >        
298 >        if ( (c_perez = calloc(5, sizeof(double))) == NULL )
299          {
300 <                fprintf(stderr,"Out of memory error in function main !");
301 <                exit(1);
300 >                fprintf(stderr,"Out of memory error in function main");
301 >                return 1;
302          }
303  
267 –
304          printhead(argc, argv);
269 –
305          computesky();
306 +        
307 +        if(*(c_perez+1)>0)
308 +        {
309 +          fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));      
310 +          print_error_sky();
311 +          exit(1);
312 +        }
313 +        
314          printsky();
315 <
273 <        exit(0);
315 >        return 0;
316   }
317  
318  
319 < void
320 < computesky()                    /* compute sky parameters */
319 >
320 >
321 >
322 >
323 >
324 >
325 > void computesky()
326   {
327  
281 –        /* new variables */
328          int     j;
329 <        float   *lv_mod;  /* 145 luminance values*/
330 <          /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
331 <        float   *theta_o, *phi_o, *coeff_perez;
329 >        
330 >        float   *lv_mod;  /* 145 luminance values */
331 >        float   *theta_o, *phi_o;
332          double  dzeta, gamma;
333          double  normfactor;
334 +        double  erbs_s0, erbs_kt;
335  
336  
290 –
337          /* compute solar direction */
338 <
338 >                
339          if (month) {                    /* from date and time */
340                  int  jd;
341                  double  sd, st;
# Line 300 | Line 346 | computesky()                   /* compute sky parameters */
346                          st = hour;
347                  else
348                          st = hour + stadj(jd);
349 +                                
350 +                
351 +                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
352 +                print_error_sky();
353 +                exit(1);
354 +                }
355 +                
356 +                                        
357 +                if(timeinterval) {
358 +                        
359 +                        if(timeinterval<0) {
360 +                        fprintf(stderr, "time interval negative\n");
361 +                        exit(1);
362 +                        }
363 +                                        
364 +                        if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) {                
365 +                        
366 +                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
367 +                        st= (st+timeinterval/120+solar_sunrise(month,day))/2;
368 +                        }
369 +                
370 +                        if(fabs(solar_sunset(month,day)-st)<timeinterval/60) {
371 +                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
372 +                        st= (st-timeinterval/120+solar_sunset(month,day))/2;
373 +                        }
374 +                }
375 +                
376 +                
377                  altitude = salt(sd, st);
378                  azimuth = sazi(sd, st);
379                  
380                  daynumber = (double)jdate(month, day);
381 +                
382 +        }
383 +        
384 +        
385  
386 +
387 +        /* if loop for the -l option. W.Sprenger (01/2013) */
388 +        /*
389 +        if (altitude*180/M_PI < sunaltitude_border) {
390 +        
391 +        if (suppress_warnings==0) {
392 +            fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
393 +                                   }
394 +        print_error_sky();
395 +        exit(1);
396          }
397 <        if (!cloudy && altitude > 87.*PI/180.) {
398 <                fprintf(stderr,
397 >        */
398 >
399 >                        
400 >        if (!cloudy && altitude > 87.*M_PI/180.) {
401 >                
402 >                if (suppress_warnings==0) {
403 >                    fprintf(stderr,
404                      "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
405                      progname);
406 <                printf(
407 <                    "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
315 <                altitude = 87.*PI/180.;
406 >                }
407 >                altitude = 87.*M_PI/180.;
408          }
409 +        
410          sundir[0] = -sin(azimuth)*cos(altitude);
411          sundir[1] = -cos(azimuth)*cos(altitude);
412          sundir[2] = sin(altitude);
413          
414                  
415          /* calculation for the new functions */
416 <        sunzenith = 90 - altitude*180/PI;
416 >        sunzenith = 90 - altitude*180/M_PI;
417          
418          
419  
420 < /* compute the inputs for the calculation of the light distribution over the sky*/
421 <        if (input==0)
420 >        /* compute the inputs for the calculation of the light distribution over the sky*/
421 >        if (input==0)           /* P */
422                  {
423                  check_parametrization();
424 <                diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
424 >                diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
425                  directirradiance = direct_irradiance_from_sky_clearness();
426                  check_irradiances();
427                  
428                  if (output==0 || output==2)
429                          {
430 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
430 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
431                          directilluminance = directirradiance*direct_n_effi_PEREZ();
432                          check_illuminances();
433                          }
434                  }
435          
436  
437 <        else if (input==1)
437 >        else if (input==1)      /* W */
438                  {
439                  check_irradiances();
440                  skybrightness = sky_brightness();
# Line 350 | Line 443 | computesky()                   /* compute sky parameters */
443  
444                  if (output==0 || output==2)
445                          {
446 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
446 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
447                          directilluminance = directirradiance*direct_n_effi_PEREZ();
448                          check_illuminances();
449                          }
# Line 358 | Line 451 | computesky()                   /* compute sky parameters */
451                  }
452                          
453          
454 <        else if (input==2)
454 >        else if (input==2)      /* L */
455                  {              
456                  check_illuminances();
457                  illu_to_irra_index();
# Line 366 | Line 459 | computesky()                   /* compute sky parameters */
459                  }
460                  
461  
462 <        else if (input==3)
462 >        else if (input==3)      /* G */
463                  {
464                          if (altitude<=0)
465                          {
466 <                        fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n");
467 <                        exit(1);
466 >                                if (suppress_warnings==0)
467 >                                     fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n");
468 >                                directirradiance = 0;
469 >                                diffuseirradiance = 0;
470 >                        } else {
471 >                        
472 >                                directirradiance=directirradiance/sin(altitude);
473                          }
474 <
377 <                directirradiance=directirradiance/sin(altitude);
474 >                                
475                  check_irradiances();
476                  skybrightness = sky_brightness();
477                  skyclearness =  sky_clearness();
# Line 382 | Line 479 | computesky()                   /* compute sky parameters */
479  
480                  if (output==0 || output==2)
481                          {
482 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
482 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
483                          directilluminance = directirradiance*direct_n_effi_PEREZ();
484                          check_illuminances();
485                          }
486  
487                  }
488  
392 –        
393 –        else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
489  
490 +        else if (input==4)      /* E */         /* Implementation of the Erbs model. W.Sprenger (04/13) */
491 +                {
492 +                        
493 +                        if (altitude<=0)
494 +                        {
495 +                                if (suppress_warnings==0 && globalirradiance > 50)
496 +                                        fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n");
497 +                                globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0;
498 +                        
499 +                        } else {
500 +                        
501 +                        erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude);
502 +                        
503 +                        if (globalirradiance>erbs_s0)
504 +                        {
505 +                                if (suppress_warnings==0)
506 +                                        fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
507 +                                
508 +                                globalirradiance=erbs_s0*0.999;                
509 +                        }
510 +                        
511 +                        erbs_kt=globalirradiance/erbs_s0;
512 +                        
513 +                        if (erbs_kt<=0.22)      diffuseirradiance=globalirradiance*(1-0.09*erbs_kt);
514 +                        else if (erbs_kt<=0.8)  diffuseirradiance=globalirradiance*(0.9511-0.1604*erbs_kt+4.388*pow(erbs_kt,2)-16.638*pow(erbs_kt,3)+12.336*pow(erbs_kt,4));
515 +                        else if (erbs_kt<1)     diffuseirradiance=globalirradiance*(0.165);
516 +                        
517 +                        directirradiance=globalirradiance-diffuseirradiance;
518 +                        
519 +                        printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
520 +                        printf("# WARNING: the -E option is only recommended for a rough estimation!");
521 +                        
522 +                        directirradiance=directirradiance/sin(altitude);
523 +                                                                                                                        
524 +                        }
525 +                        
526 +                check_irradiances();
527 +                skybrightness = sky_brightness();
528 +                skyclearness =  sky_clearness();
529 +                check_parametrization();
530  
531 +                if (output==0 || output==2)
532 +                        {
533 +                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
534 +                        directilluminance = directirradiance*direct_n_effi_PEREZ();
535 +                        check_illuminances();
536 +                        }
537 +
538 +                }
539 +                
540 +                
541 +                
542          
543 < /* normalization factor for the relative sky luminance distribution, diffuse part*/
543 >        else    {fprintf(stderr,"error at the input arguments"); exit(1);}
544  
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 –        }
545  
546 < /* read the coefficients for the Perez sky luminance model */
547 <        if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
548 <        {
409 <                fprintf(stderr,"lect_coeff_perez does not work\n");
410 <                exit(2);
411 <        }
412 <
546 >        
547 >        /* normalization factor for the relative sky luminance distribution, diffuse part*/
548 >        
549          if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
550          {
551                  fprintf(stderr,"Out of memory in function main");
# Line 417 | Line 553 | computesky()                   /* compute sky parameters */
553          }
554  
555          /* read the angles */
556 <        theta_o = theta_ordered("defangle.dat");
557 <        phi_o = phi_ordered("defangle.dat");
556 >        theta_o = defangle_theta;
557 >        phi_o = defangle_phi;
558 >        
559  
560 < /* parameters for the perez model */
560 >        /* parameters for the perez model */
561          coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
562  
563 < /*calculation of the modelled luminance */
563 >        
564 >        
565 >        
566 >        
567 >        /*calculation of the modelled luminance */
568          for (j=0;j<145;j++)
569          {
570                  theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
571 +                                
572                  *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
573 <                /*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
573 >                
574 >                /* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
575          }
576 <
576 >        
577          /* integration of luminance for the normalization factor, diffuse part of the sky*/
578 +        
579          diffnormalization = integ_lv(lv_mod, theta_o);
436 –        /*printf("perez integration %lf\n", diffnormalization);*/
580          
438 –        
581  
582  
583 < /*normalization coefficient in lumen or in watt*/
583 >        /*normalization coefficient in lumen or in watt*/
584          if (output==0)
585                  {
586 <                diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
586 >                diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
587                  }
588          else if (output==1)
589                  {
590 <                diffnormalization = diffusirradiance/diffnormalization;
590 >                diffnormalization = diffuseirradiance/diffnormalization;
591                  }
592          else if (output==2)
593                  {
594 <                diffnormalization = diffusilluminance/diffnormalization;
594 >                diffnormalization = diffuseilluminance/diffnormalization;
595                  }
596  
597 <        else    {fprintf(stderr,"output argument : wrong number"); exit(1);}
597 >        else    {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
598  
599  
600  
601  
602 < /* calculation for the solar source */  
602 >        /* calculation for the solar source */  
603          if (output==0)
604 <                solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)))/WHTEFFICACY;
604 >                solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
605                  
606          else if (output==1)
607 <                solarradiance = directirradiance/(2*PI*(1-cos(half_sun_angle*PI/180)));
607 >                solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
608          
609          else
610 <                solarradiance = directilluminance/(2*PI*(1-cos(half_sun_angle*PI/180)));
610 >                solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
611          
612                          
613  
# Line 473 | Line 615 | computesky()                   /* compute sky parameters */
615   /* Compute the ground radiance */
616   zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
617   zenithbr*=diffnormalization;
618 < /*
477 < fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
478 < */
479 <
618 >        
619   if (skyclearness==1)
620          normfactor = 0.777778;
621                  
622   if (skyclearness>=6)
623          {              
624 <        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
625 <        normfactor = normsc()/F2/PI;
624 >        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
625 >        normfactor = normsc()/F2/M_PI;
626          }
627  
628   if ( (skyclearness>1) && (skyclearness<6) )
629          {
630          S_INTER=1;
631 <        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(PI/2.0-altitude)*(.4441+1.48*altitude));
632 <        normfactor = normsc()/F2/PI;
631 >        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
632 >        normfactor = normsc()/F2/M_PI;
633          }
634  
635   groundbr = zenithbr*normfactor;
497 – printf("# Ground ambient level: %.1f\n", groundbr);
636  
637   if (dosun&&(skyclearness>1))
638 < groundbr += 6.8e-5/PI*solarradiance*sundir[2];          
638 >        groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
639  
640   groundbr *= gprefl;
641  
# Line 509 | Line 647 | return;
647  
648  
649  
650 + void print_error_sky()
651 + {
652 +        sundir[0] = -sin(azimuth)*cos(altitude);
653 +        sundir[1] = -cos(azimuth)*cos(altitude);
654 +        sundir[2] = sin(altitude);
655 +        
656 +        printf("\nvoid brightfunc skyfunc\n");
657 +        printf("2 skybright perezlum.cal\n");
658 +        printf("0\n");
659 +        printf("10 0.00 0.00  0.000 0.000 0.000 0.000 0.000  %f %f %f \n", sundir[0], sundir[1], sundir[2]);
660 + }
661 +        
662  
663  
664  
665 < printsky()                      /* print out sky */
665 >
666 >
667 >
668 > double solar_sunset(int month,int day)
669   {
670 +     float W;
671 +     extern double s_latitude;
672 +     W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
673 +     return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
674 + }
675 +
676 +
677 + double solar_sunrise(int month,int day)
678 + {
679 +     float W;
680 +     extern double s_latitude;
681 +     W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
682 +     return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
683 + }
684 +
685 +
686 +
687 +
688 +
689 +
690 +
691 +
692 +
693 +
694 +
695 + void printsky()                 /* print out sky */
696 + {
697          if (dosun&&(skyclearness>1))
698 <                {              
698 >        {              
699                  printf("\nvoid light solar\n");
700                  printf("0\n0\n");
701                  printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
702                  printf("\nsolar source sun\n");
703                  printf("0\n0\n");
704                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
705 <                }
526 <                
527 <        if (dosun&&(skyclearness==1))
528 <                {              
705 >        } else if (dosun) {
706                  printf("\nvoid light solar\n");
707                  printf("0\n0\n");
708                  printf("3 0.0 0.0 0.0\n");
709                  printf("\nsolar source sun\n");
710                  printf("0\n0\n");
711                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
712 <                }
712 >        }
713          
714  
715          printf("\nvoid brightfunc skyfunc\n");
716          printf("2 skybright perezlum.cal\n");
717          printf("0\n");
718          printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
719 <            *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
720 <            sundir[0], sundir[1], sundir[2]);
719 >                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
720 >                sundir[0], sundir[1], sundir[2]);
721 >        
722   }
723  
724  
725 < printdefaults()                 /* print default values */
725 > void printdefaults()                    /* print default values */
726   {
727          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
728          if (zenithbr > 0.0)
729                  printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
730          else
731                  printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity);
732 <        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/PI));
733 <        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/PI));
734 <        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/PI));
732 >        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI));
733 >        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI));
734 >        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI));
735   }
736  
737  
738 < userror(msg)                    /* print usage error and quit */
561 < char  *msg;
738 > void userror(char* msg)                 /* print usage error and quit */
739   {
740          if (msg != NULL)
741 <                fprintf(stderr, "%s: Use error - %s\n", progname, msg);
742 <        fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
743 <        fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
741 >                fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
742 >        fprintf(stderr, "Usage: %s      month day hour    [...]\n", progname);
743 >        fprintf(stderr, "   or: %s -ang altitude azimuth  [...]\n", progname);
744 >        fprintf(stderr, "               followed by:      -P          epsilon delta [options]\n");
745 >        fprintf(stderr, "                        or:      [-W|-L|-G]  direct_value diffuse_value [options]\n");
746 >        fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
747 >        fprintf(stderr, "       Description:\n");
748          fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
749          fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
750          fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
751          fprintf(stderr, "       -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
752 +        fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
753 +        fprintf(stderr, "       Output specification with option:\n");
754          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");
755 +        fprintf(stderr, "       gendaylit version 2.3 (2013/08/08)  \n\n");
756          exit(1);
757   }
758  
759  
760  
761 < double
578 < normsc()                        /* compute normalization factor (E0*F2/L0) */
761 > double normsc()           /* compute normalization factor (E0*F2/L0) */
762   {
763          static double  nfc[2][5] = {
764                                  /* clear sky approx. */
# Line 588 | Line 771 | normsc()                       /* compute normalization factor (E0*F2/L0)
771          register int  i;
772                                          /* polynomial approximation */
773          nf = nfc[S_INTER];
774 <        x = (altitude - PI/4.0)/(PI/4.0);
774 >        x = (altitude - M_PI/4.0)/(M_PI/4.0);
775          nsc = nf[i=4];
776          while (i--)
777                  nsc = nsc*x + nf[i];
# Line 598 | Line 781 | normsc()                       /* compute normalization factor (E0*F2/L0)
781  
782  
783  
784 < printhead(ac, av)               /* print command header */
602 < register int  ac;
603 < register char  **av;
784 > void printhead(int ac, char** av)               /* print command header */
785   {
786          putchar('#');
787          while (ac--) {
# Line 613 | Line 794 | register char  **av;
794  
795  
796  
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 –
797   /* Perez models */
798  
799   /* Perez global horizontal luminous efficacy model */
# Line 640 | Line 803 | double glob_h_effi_PEREZ()
803          double  value;
804          double  category_bounds[10], a[10], b[10], c[10], d[10];
805          int     category_total_number, category_number, i;
806 <
807 <
808 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
809 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
810 <
806 >        
807 >        check_parametrization();
808 >        
809 >        
810 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
811 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */
812 >    
813 >    
814          /* initialize category bounds (clearness index bounds) */
815  
816          category_total_number = 8;
# Line 707 | Line 873 | fprintf(stderr, "Warning : skyclearness or skybrightne
873          }
874  
875          value = a[category_number] + b[category_number]*atm_preci_water  +
876 <            c[category_number]*cos(sunzenith*PI/180) +  d[category_number]*log(skybrightness);
876 >            c[category_number]*cos(sunzenith*M_PI/180) +  d[category_number]*log(skybrightness);
877  
878          return(value);
879   }
# Line 720 | Line 886 | double glob_h_diffuse_effi_PEREZ()
886          double  category_bounds[10], a[10], b[10], c[10], d[10];
887          int     category_total_number, category_number, i;
888  
889 +
890 +
891          
892 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
893 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
892 >        check_parametrization();
893 >        
894 >        
895 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
896 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
897 >    
898 >
899  
900   /* initialize category bounds (clearness index bounds) */
901  
902          category_total_number = 8;
903  
904 + //XXX:  category_bounds > 0.1
905          category_bounds[1] = 1;
906          category_bounds[2] = 1.065;
907          category_bounds[3] = 1.230;
# Line 779 | Line 953 | fprintf(stderr, "Warning : skyclearness or skybrightne
953  
954  
955  
956 +        category_number = -1;
957          for (i=1; i<=category_total_number; i++)
958          {
959                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
960                          category_number = i;
961          }
962  
963 <        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*PI/180) +
963 >        if (category_number == -1) {
964 >                if (suppress_warnings==0)
965 >                    fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness);
966 >                print_error_sky();
967 >                exit(1);
968 >        }
969 >                
970 >
971 >        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) +
972              d[category_number]*log(skybrightness);
973  
974          return(value);
975 +
976   }
977  
978  
979 +
980   /* direct normal efficacy model, according to PEREZ */
981  
982   double direct_n_effi_PEREZ()
# Line 802 | Line 987 | double         category_bounds[10], a[10], b[10], c[10], d[10
987   int     category_total_number, category_number, i;
988  
989  
990 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
991 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
990 > if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
991 >   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");
992  
993  
994   /* initialize category bounds (clearness index bounds) */
# Line 866 | Line 1051 | if ( (skyclearness >= category_bounds[i]) && (skyclear
1051   category_number = i;
1052   }
1053  
1054 < value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*exp(5.73*sunzenith*PI/180 - 5) +  d[category_number]*skybrightness;
1054 > 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;
1055  
1056   if (value < 0) value = 0;
1057  
# Line 877 | Line 1062 | return(value);
1062   /*check the range of epsilon and delta indexes of the perez parametrization*/
1063   void check_parametrization()
1064   {
1065 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
1065 > if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1066                  {
1067 <                fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
1068 <                exit(1);        
1067 >
1068 > /*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1069 >                if (skyclearness<skyclearinf){
1070 >                        if (suppress_warnings==0)
1071 >                            /* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1072 >                        skyclearness=skyclearinf;
1073                  }
1074 +                if (skyclearness>skyclearsup){
1075 +                        if (suppress_warnings==0)
1076 +                            /* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1077 +                        skyclearness=skyclearsup-0.1;
1078 +                }
1079 +                if (skybrightness<skybriginf){
1080 +                        if (suppress_warnings==0)
1081 +                            /* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1082 +                        skybrightness=skybriginf;
1083 +                }
1084 +                if (skybrightness>skybrigsup){
1085 +                        if (suppress_warnings==0)
1086 +                            /* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1087 +                        skybrightness=skybrigsup;
1088 +                }
1089 +
1090 +        return; }
1091          else return;
1092   }
1093  
1094  
1095 < /* likelihood of the direct and diffuse components */
1095 > /* validity of the direct and diffuse components */
1096   void    check_illuminances()
1097   {
1098 <        if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
1099 <        {
1100 <        fprintf(stderr,"direct or diffuse illuminances out of range\n");
895 <        exit(1);
1098 >        if (directilluminance < 0) {
1099 >                fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
1100 >                directilluminance = 0.0;
1101          }
1102 < return;
1102 >        if (diffuseilluminance < 0) {
1103 >                fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n");
1104 >                diffuseilluminance = 0.0;
1105 >        }
1106 >        if (directilluminance > solar_constant_l*1000.0) {
1107 >                fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n");
1108 >                exit(1);
1109 >        }
1110   }
1111  
1112  
1113   void    check_irradiances()
1114   {
1115 <        if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
1116 <        {
1117 <        fprintf(stderr,"direct or diffuse irradiances out of range\n");
1118 <        exit(1);
1119 <        }      
1120 < return;
1115 >        if (directirradiance < 0) {
1116 >                fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n");
1117 >                directirradiance = 0.0;
1118 >        }
1119 >        if (diffuseirradiance < 0) {
1120 >                fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n");
1121 >                diffuseirradiance = 0.0;
1122 >        }
1123 >        if (directirradiance > solar_constant_e) {
1124 >                fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n");
1125 >                exit(1);
1126 >        }
1127   }
1128          
1129  
# Line 915 | Line 1133 | double sky_brightness()
1133   {
1134   double value;
1135  
1136 < value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1136 > value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1137  
1138   return(value);
1139   }
# Line 924 | Line 1142 | return(value);
1142   /* Perez sky's clearness */
1143   double sky_clearness()
1144   {
1145 < double value;
1145 >        double value;
1146  
1147 < 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) ;
1147 >        value = ( (diffuseirradiance + directirradiance)/(diffuseirradiance) + 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) ;
1148  
1149 < return(value);
1149 >        return(value);
1150   }
1151  
1152  
1153  
1154   /* diffus horizontal irradiance from Perez sky's brightness */
1155 < double diffus_irradiance_from_sky_brightness()
1155 > double diffuse_irradiance_from_sky_brightness()
1156   {
1157          double value;
1158  
# Line 949 | Line 1167 | double direct_irradiance_from_sky_clearness()
1167   {
1168          double value;
1169  
1170 <        value = diffus_irradiance_from_sky_brightness();
1171 <        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*PI/180*sunzenith*PI/180*sunzenith*PI/180) );
1170 >        value = diffuse_irradiance_from_sky_brightness();
1171 >        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
1172  
1173          return(value);
1174   }
1175  
1176  
1177 < void illu_to_irra_index(void)
1177 >
1178 >
1179 > void illu_to_irra_index()
1180   {
1181 < double  test1=0.1, test2=0.1;
1181 > double  test1=0.1, test2=0.1, d_eff;
1182   int     counter=0;      
1183  
1184 < diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
1184 > diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000);
1185   directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1186   skyclearness =  sky_clearness();
1187   skybrightness = sky_brightness();
1188 < 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);*/
1188 > check_parametrization();
1189  
1190 <        test1=diffusirradiance;
1190 >
1191 > while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1192 >                || (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1193 >                || (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1194 >                 && !(counter==9) )
1195 >        {
1196 >        
1197 >        test1=diffuseirradiance;
1198          test2=directirradiance;
1199          counter++;
1200          
1201 <        diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
1202 <        directirradiance = directilluminance/direct_n_effi_PEREZ();
985 <        /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1201 >        diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1202 >        d_eff = direct_n_effi_PEREZ();
1203          
1204 +        
1205 +        if (d_eff < 0.1)
1206 +                directirradiance = 0;
1207 +        else    
1208 +                directirradiance = directilluminance/d_eff;
1209 +        
1210          skybrightness = sky_brightness();
1211          skyclearness =  sky_clearness();
1212 <        if (skyclearness>12) skyclearness=12;
1213 <        if (skybrightness<0.05) skybrightness=0.01;
1214 <
1215 <                /*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
993 <
1212 >        check_parametrization();
1213 >        
1214 >        /*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/
1215 >                
1216          }
1217  
1218  
1219   return;
1220   }              
1221  
1222 <
1001 < int lect_coeff_perez(char *filename,float **coeff_perez)
1222 > static int get_numlin(float epsilon)
1223   {
1224 <        FILE *fcoeff_perez;
1225 <        float temp;
1226 <        int i,j;
1227 <
1228 <        if ((fcoeff_perez = frlibopen(filename)) == NULL)
1229 <        {
1230 <                fprintf(stderr,"file %s cannot be opened\n", filename);
1231 <                return 1; /* il y a un probleme de fichier */
1232 <        }
1233 <        else
1234 <        {
1235 <                /*printf("file %s  open\n", filename);*/
1236 <        }
1237 <        
1238 <        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 */
1224 >        if (epsilon < 1.065)
1225 >                return 0;
1226 >        else if (epsilon < 1.230)
1227 >                return 1;
1228 >        else if (epsilon < 1.500)
1229 >                return 2;
1230 >        else if (epsilon < 1.950)
1231 >                return 3;
1232 >        else if (epsilon < 2.800)
1233 >                return 4;
1234 >        else if (epsilon < 4.500)
1235 >                return 5;
1236 >        else if (epsilon < 6.200)
1237 >                return 6;
1238 >        return 7;
1239   }
1240  
1030 –
1031 –
1241   /* sky luminance perez model */
1242 < double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1034 < double epsilon,double Delta,float *coeff_perez)
1242 > double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1243   {
1244 +                        
1245          float x[5][4];
1246          int i,j,num_lin;
1247          double c_perez[5];
1248  
1249          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1250          {
1251 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1251 >                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n");
1252                  exit(1);
1253          }
1254  
# Line 1048 | Line 1257 | double epsilon,double Delta,float *coeff_perez)
1257          {
1258                  if ( Delta < 0.2 ) Delta = 0.2;
1259          }
1260 <
1261 <        if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1262 <        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1263 <        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 <
1260 >        
1261 >        
1262 >        num_lin = get_numlin(epsilon);
1263 >        
1264          for (i=0;i<5;i++)
1265                  for (j=0;j<4;j++)
1266                  {
1267                          x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1268 <                        /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1268 >                        /* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1269                  }
1270  
1271  
# Line 1089 | Line 1292 | double epsilon,double Delta,float *coeff_perez)
1292  
1293  
1294   /* coefficients for the sky luminance perez model */
1295 < void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1295 > void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1296   {
1297          float x[5][4];
1298          int i,j,num_lin;
# Line 1105 | Line 1308 | void coeff_lum_perez(double Z, double epsilon, double
1308          {
1309                  if ( Delta < 0.2 ) Delta = 0.2;
1310          }
1311 +        
1312 +        
1313 +        num_lin = get_numlin(epsilon);
1314  
1315 <        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;
1315 >        /*fprintf(stderr,"numlin %d\n", num_lin);*/
1316  
1317          for (i=0;i<5;i++)
1318                  for (j=0;j<4;j++)
# Line 1145 | Line 1344 | void coeff_lum_perez(double Z, double epsilon, double
1344   }
1345  
1346  
1347 +
1348   /* degrees into radians */
1349   double radians(double degres)
1350   {
1351 <        return degres*PI/180.0;
1351 >        return degres*M_PI/180.0;
1352   }
1353  
1354 +
1355   /* radian into degrees */
1356   double degres(double radians)
1357   {
1358 <        return radians/PI*180.0;
1358 >        return radians/M_PI*180.0;
1359   }
1360  
1361 +
1362   /* calculation of the angles dzeta and gamma */
1363   void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1364   {
# Line 1173 | Line 1375 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1375   }
1376  
1377  
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;
1378  
1379 <        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)
1379 > double integ_lv(float *lv,float *theta)
1380   {
1381          int i;
1382 <        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 <        }
1382 >        double buffer=0.0;
1383          
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 –
1384          for (i=0;i<145;i++)
1385          {
1386 <                fscanf(file_in,"%f",&buffer);
1259 <                fscanf(file_in,"%f",&buffer);
1260 <                *(ptr+i) = buffer;
1386 >                buffer += (*(lv+i))*cos(radians(*(theta+i)));
1387          }
1388 <
1389 <        fclose(file_in);
1264 <        return ptr;
1388 >                        
1389 >        return buffer*2*M_PI/144;
1390   }
1391  
1392  
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;
1393  
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 –
1394   /* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1395  
1396   double get_eccentricity()
# Line 1304 | Line 1398 | double get_eccentricity()
1398          double day_angle;
1399          double E0;
1400  
1401 <        day_angle  = 2*PI*(daynumber - 1)/365;
1401 >        day_angle  = 2*M_PI*(daynumber - 1)/365;
1402          E0         = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
1403              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1404  
1405          return (E0);
1312 –
1406   }
1407  
1408  
# Line 1317 | Line 1410 | double get_eccentricity()
1410   double  air_mass()
1411   {
1412   double  m;
1320 –
1413   if (sunzenith>90)
1414          {
1415 <        fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1415 >        fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1416          exit(1);
1417          }
1418 <        
1327 < m = 1/( cos(sunzenith*PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1418 > m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1419   return(m);
1420   }
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 –
1421  
1422  
1423  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines