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

Comparing ray/src/gen/gendaylit.c (file contents):
Revision 2.2 by greg, Mon Jun 15 22:27:21 2009 UTC vs.
Revision 2.16 by greg, Wed Jul 30 17:30:27 2014 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  
9 <
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 > #define  _USE_MATH_DEFINES
10   #include  <stdio.h>
11   #include  <string.h>
12   #include  <math.h>
13   #include  <stdlib.h>
24 – #include  <ctype.h>
14  
26 – #include  "rtio.h"
27 – #include  "fvect.h"
15   #include  "color.h"
16 + #include  "sun.h"
17   #include  "paths.h"
18  
19 < 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);
19 > #define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
20  
21   double  normsc();
22  
23 < #define DATFILE         "coeff_perez.dat"
23 > /*static        char *rcsid="$Header$";*/
24  
25 + float coeff_perez[] = {
26 +        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,
27 +        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,
28 +        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,
29 +        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,
30 +        -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,
31 +        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,
32 +        -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,
33 +        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,
34 +        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,
35 +        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,
36 +        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,
37 +        -14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636};
38  
39  
40 < /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
40 > float defangle_theta[] = {
41 >        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,
42 >        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,
43 >        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,
44 >        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,
45 >        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,
46 >        24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};
47 >
48 > float defangle_phi[] = {
49 >        0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264,
50 >        276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180,
51 >        192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105,
52 >        120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75,
53 >        90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
54 >        80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
55 >        240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
56 >
57 >
58 >
59 > /* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
60   double sky_brightness();
61   double sky_clearness();
62  
63   /* calculation of the direct and diffuse components from the Perez parametrization */
64 < double  diffus_irradiance_from_sky_brightness();
64 > double  diffuse_irradiance_from_sky_brightness();
65   double  direct_irradiance_from_sky_clearness();
66  
67 + /* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */
68 + /* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
69  
52 – /* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
70   double  glob_h_effi_PEREZ();
71   double  glob_h_diffuse_effi_PEREZ();
72   double  direct_n_effi_PEREZ();
73 +
74   /*likelihood check of the epsilon, delta, direct and diffuse components*/
75   void    check_parametrization();
76   void    check_irradiances();
77   void    check_illuminances();
78   void    illu_to_irra_index();
79 + void    print_error_sky();
80  
81 <
82 < /* 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);
81 > double  calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
82 > void    coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]);
83   double  radians(double degres);
84   double  degres(double radians);
85   void    theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);
86   double  integ_lv(float *lv,float *theta);
73 – float   *theta_ordered(char *filename);
74 – float   *phi_ordered(char *filename);
75 – void    skip_comments(FILE *fp);
87  
88 + void printdefaults();
89 + void check_sun_position();
90 + void computesky();
91 + void printhead(int ac, char** av);
92 + void usage_error(char* msg);
93 + void printsky();
94  
95 + FILE * frlibopen(char* fname);
96  
97   /* astronomy and geometry*/
98   double  get_eccentricity();
99   double  air_mass();
82 – double  get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
100  
101 <
102 < /* date*/
101 > double  solar_sunset(int month, int day);
102 > double  solar_sunrise(int month, int day);
103 > double  stadj();
104   int     jdate(int month, int day);
105  
88 –
89 –
90 –
91 –
92 – /* sun calculation constants */
93 – extern double  s_latitude;
94 – extern double  s_longitude;
95 – extern double  s_meridian;
96 –
106   const double    AU = 149597890E3;
107   const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
108 < const double    solar_constant_l = 127.5;   /* solar constant klux */
108 > const double    solar_constant_l = 127500;   /* solar constant lux */
109  
110   const double    half_sun_angle = 0.2665;
111   const double    half_direct_angle = 2.85;
112  
113 < const double    skyclearinf = 1.000;    /* limitations for the variation of the Perez parameters */
114 < const double    skyclearsup = 12.1;
113 > const double    skyclearinf = 1.0;          /* limitations for the variation of the Perez parameters */
114 > const double    skyclearsup = 12.01;
115   const double    skybriginf = 0.01;
116   const double    skybrigsup = 0.6;
117  
# Line 117 | Line 126 | double  altitude, azimuth;                     /* or solar angles */
126  
127  
128   /* definition of the sky conditions through the Perez parametrization */
129 < double  skyclearness, skybrightness;
130 < double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
131 < double  diffusilluminance, directilluminance, diffusirradiance, directirradiance;
132 < double  sunzenith, daynumber=150, atm_preci_water=2;
129 > double  skyclearness = 0;
130 > double  skybrightness = 0;
131 > double  solarradiance;
132 > double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
133 > double  sunzenith, daynumber, atm_preci_water=2;
134  
135 < double  diffnormalization, dirnormalization;
135 > /*double  sunaltitude_border = 0;*/
136 > double  diffnormalization = 0;
137 > double  dirnormalization = 0;
138   double  *c_perez;
139  
140 < int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
141 < int     input=0;        /*define the input for the calulation*/
140 > int     output=0;       /* define the unit of the output (sky luminance or radiance): */
141 >                        /* visible watt=0, solar watt=1, lumen=2 */
142 > int     input=0;        /* define the input for the calulation */
143  
144 + int     suppress_warnings=0;
145 +
146          /* default values */
147 < int  cloudy = 0;                                /* 1=standard, 2=uniform */
148 < int  dosun = 1;
147 > int     cloudy = 0;                             /* 1=standard, 2=uniform */
148 > int     dosun = 1;
149   double  zenithbr = -1.0;
150   double  betaturbidity = 0.1;
151   double  gprefl = 0.2;
# Line 138 | Line 153 | int    S_INTER=0;
153  
154          /* computed values */
155   double  sundir[3];
156 < double  groundbr;
156 > double  groundbr = 0;
157   double  F2;
158   double  solarbr = 0.0;
159   int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
160 + float   timeinterval = 0;
161  
162 < char  *progname;
163 < char  errmsg[128];
162 > char    *progname;
163 > char    errmsg[128];
164  
165 + double  st;
166  
167 < main(argc, argv)
168 < int  argc;
152 < char  *argv[];
167 >
168 > int main(int argc, char** argv)
169   {
170          int  i;
171  
172          progname = argv[0];
173          if (argc == 2 && !strcmp(argv[1], "-defaults")) {
174                  printdefaults();
175 <                exit(0);
175 >                return 0;
176          }
177          if (argc < 4)
178 <                userror("arg count");
178 >                usage_error("arg count");
179          if (!strcmp(argv[1], "-ang")) {
180                  altitude = atof(argv[2]) * (M_PI/180);
181                  azimuth = atof(argv[3]) * (M_PI/180);
# Line 167 | Line 183 | char  *argv[];
183          } else {
184                  month = atoi(argv[1]);
185                  if (month < 1 || month > 12)
186 <                        userror("bad month");
186 >                        usage_error("bad month");
187                  day = atoi(argv[2]);
188                  if (day < 1 || day > 31)
189 <                        userror("bad day");
189 >                        usage_error("bad day");
190                  hour = atof(argv[3]);
191                  if (hour < 0 || hour >= 24)
192 <                        userror("bad hour");
192 >                        usage_error("bad hour");
193                  tsolar = argv[3][0] == '+';
194          }
195          for (i = 4; i < argc; i++)
# Line 183 | Line 199 | char  *argv[];
199                                  cloudy = 0;
200                                  dosun = argv[i][0] == '+';
201                                  break;
186 –                        case 'r':
202                          case 'R':
203                                  u_solar = argv[i][1] == 'R' ? -1 : 1;
204                                  solarbr = atof(argv[++i]);
# Line 195 | Line 210 | char  *argv[];
210                          case 't':
211                                  betaturbidity = atof(argv[++i]);
212                                  break;
213 +                        case 'w':
214 +                                suppress_warnings = 1;
215 +                                break;                  
216                          case 'b':
217                                  zenithbr = atof(argv[++i]);
218                                  break;
# Line 210 | Line 228 | char  *argv[];
228                          case 'm':
229                                  s_meridian = atof(argv[++i]) * (M_PI/180);
230                                  break;
213 –
231                          
232                          case 'O':
233 <                                output = atof(argv[++i]);       /*define the unit of the output of the program :
234 <                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm)
235 <                                                                default is set to 0*/
233 >                                output = atof(argv[++i]);       /*define the unit of the output of the program:
234 >                                                                sky and sun luminance/radiance
235 >                                                                (0==W visible, 1==W solar radiation, 2==lm) */
236                                  break;
237                                  
238                          case 'P':
# Line 227 | Line 244 | char  *argv[];
244                          case 'W':                                       /* direct normal Irradiance [W/m^2] */
245                                  input = 1;                              /* diffuse horizontal Irrad. [W/m^2] */
246                                  directirradiance = atof(argv[++i]);
247 <                                diffusirradiance = atof(argv[++i]);
247 >                                diffuseirradiance = atof(argv[++i]);
248                                  break;
249                                  
250                          case 'L':                                       /* direct normal Illuminance [Lux] */
251                                  input = 2;                              /* diffuse horizontal Ill. [Lux] */
252                                  directilluminance = atof(argv[++i]);
253 <                                diffusilluminance = atof(argv[++i]);
253 >                                diffuseilluminance = atof(argv[++i]);
254                                  break;
255                          
256                          case 'G':                                       /* direct horizontal Irradiance [W/m^2] */
257                                  input = 3;                              /* diffuse horizontal Irrad. [W/m^2] */
258                                  directirradiance = atof(argv[++i]);
259 <                                diffusirradiance = atof(argv[++i]);
259 >                                diffuseirradiance = atof(argv[++i]);
260                                  break;
244 –                                
261                          
262 +                        case 'E':                                       /* Erbs model based on the */
263 +                                input = 4;                              /* global-horizontal irradiance [W/m^2] */
264 +                                globalirradiance = atof(argv[++i]);
265 +                                break;
266 +                        
267 +                        case 'i':
268 +                                timeinterval = atof(argv[++i]);
269 +                                break;
270 +                        
271 +                        
272                          default:
273                                  sprintf(errmsg, "unknown option: %s", argv[i]);
274 <                                userror(errmsg);
274 >                                usage_error(errmsg);
275                          }
276                  else
277 <                        userror("bad option");
277 >                        usage_error("bad option");
278  
279 <        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
280 <                fprintf(stderr,
255 <                    "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
279 >        if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180)
280 >                fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
281                      progname, (s_longitude-s_meridian)*12/M_PI);
282  
283  
284 <        /* allocation dynamique de memoire pour les pointeurs */
285 <        if ( (c_perez = malloc(5*sizeof(double))) == NULL )
286 <        {
262 <                fprintf(stderr,"Out of memory error in function main !");
263 <                exit(1);
264 <        }
284 >        /* dynamic memory allocation for the pointers */
285 >        if ( (c_perez = calloc(5, sizeof(double))) == NULL )
286 >        { fprintf(stderr,"Out of memory error in function main"); return 1; }
287  
288 <
288 >        
289          printhead(argc, argv);
268 –
290          computesky();
291          printsky();
292 +        return 0;
293  
272 –        exit(0);
294   }
295  
296  
297 < computesky()                    /* compute sky parameters */
297 >
298 >
299 >
300 > void computesky()
301   {
302  
303 <        /* new variables */
304 <        int     j, i;
305 <        float   *lv_mod;  /* 145 luminance values*/
306 <          /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
283 <        float   *theta_o, *phi_o, *coeff_perez;
303 >        int     j;
304 >        
305 >        float   *lv_mod;  /* 145 luminance values */
306 >        float   *theta_o, *phi_o;
307          double  dzeta, gamma;
285 –        double  diffusion;
308          double  normfactor;
309 +        double  erbs_s0, erbs_kt;
310  
311  
289 –
312          /* compute solar direction */
313 <
313 >                
314          if (month) {                    /* from date and time */
315                  int  jd;
316 <                double  sd, st;
316 >                double  sd;
317  
318                  jd = jdate(month, day);         /* Julian date */
319                  sd = sdec(jd);                  /* solar declination */
# Line 299 | Line 321 | computesky()                   /* compute sky parameters */
321                          st = hour;
322                  else
323                          st = hour + stadj(jd);
324 +                
325 +                                        
326 +                if(timeinterval) {
327 +                        
328 +                        if(timeinterval<0) {
329 +                        fprintf(stderr, "time interval negative\n");
330 +                        exit(1);
331 +                        }
332 +                                                                        
333 +                        if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) {                      
334 +                         st= (st+timeinterval/120+solar_sunrise(month,day))/2;
335 +                         if(suppress_warnings==0)
336 +                         { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
337 +                        }
338 +                
339 +                        if(fabs(solar_sunset(month,day)-st)<timeinterval/120) {
340 +                         st= (st-timeinterval/120+solar_sunset(month,day))/2;
341 +                         if(suppress_warnings==0)
342 +                         { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
343 +                        }
344 +                        
345 +                        if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) {
346 +                          if(suppress_warnings==0)
347 +                          { fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); }
348 +                         altitude = salt(sd, st);
349 +                         azimuth = sazi(sd, st);
350 +                         print_error_sky();
351 +                         exit(0);
352 +                        }
353 +                }
354 +                else
355 +                
356 +                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
357 +                        if(suppress_warnings==0)
358 +                        { fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); }
359 +                        altitude = salt(sd, st);
360 +                        azimuth = sazi(sd, st);
361 +                        print_error_sky();
362 +                        exit(0);
363 +                }
364 +                
365                  altitude = salt(sd, st);
366                  azimuth = sazi(sd, st);
367                  
368                  daynumber = (double)jdate(month, day);
369 <
369 >                
370          }
371 +        
372 +        
373 +
374 +
375 +                        
376          if (!cloudy && altitude > 87.*M_PI/180.) {
377 <                fprintf(stderr,
377 >                
378 >                if (suppress_warnings==0) {
379 >                    fprintf(stderr,
380                      "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
381                      progname);
382 <                printf(
313 <                    "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
382 >                }
383                  altitude = 87.*M_PI/180.;
384          }
385 +        
386 +        
387 +        
388          sundir[0] = -sin(azimuth)*cos(altitude);
389          sundir[1] = -cos(azimuth)*cos(altitude);
390          sundir[2] = sin(altitude);
# Line 320 | Line 392 | computesky()                   /* compute sky parameters */
392                  
393          /* calculation for the new functions */
394          sunzenith = 90 - altitude*180/M_PI;
395 <        
324 <        
395 >                        
396  
397 < /* compute the inputs for the calculation of the light distribution over the sky*/
398 <        if (input==0)
397 >        /* compute the inputs for the calculation of the light distribution over the sky*/
398 >        if (input==0)           /* P */
399                  {
400                  check_parametrization();
401 <                diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
401 >                diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
402                  directirradiance = direct_irradiance_from_sky_clearness();
403                  check_irradiances();
404                  
405                  if (output==0 || output==2)
406                          {
407 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
407 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
408                          directilluminance = directirradiance*direct_n_effi_PEREZ();
409                          check_illuminances();
410                          }
411                  }
412          
413  
414 <        else if (input==1)
414 >        else if (input==1)      /* W */
415                  {
416                  check_irradiances();
417                  skybrightness = sky_brightness();
418                  skyclearness =  sky_clearness();
419 +                
420                  check_parametrization();
421 <
421 >                                                        
422                  if (output==0 || output==2)
423                          {
424 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
424 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
425                          directilluminance = directirradiance*direct_n_effi_PEREZ();
426                          check_illuminances();
427                          }
# Line 357 | Line 429 | computesky()                   /* compute sky parameters */
429                  }
430                          
431          
432 <        else if (input==2)
432 >        else if (input==2)      /* L */
433                  {              
434                  check_illuminances();
435                  illu_to_irra_index();
# Line 365 | Line 437 | computesky()                   /* compute sky parameters */
437                  }
438                  
439  
440 <        else if (input==3)
440 >        else if (input==3)      /* G */
441                  {
442                          if (altitude<=0)
443                          {
444 <                        fprintf(stderr, "solar zenith angle larger than 90� \n the models used are not more valid\n");
445 <                        exit(1);
444 >                                if (suppress_warnings==0)
445 >                                     fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
446 >                                directirradiance = 0;
447 >                                diffuseirradiance = 0;
448 >                        } else {
449 >                        
450 >                                directirradiance=directirradiance/sin(altitude);
451                          }
452 <
376 <                directirradiance=directirradiance/sin(altitude);
452 >                                
453                  check_irradiances();
454                  skybrightness = sky_brightness();
455                  skyclearness =  sky_clearness();
# Line 381 | Line 457 | computesky()                   /* compute sky parameters */
457  
458                  if (output==0 || output==2)
459                          {
460 <                        diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
460 >                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
461                          directilluminance = directirradiance*direct_n_effi_PEREZ();
462                          check_illuminances();
463                          }
464  
465                  }
466  
391 –        
392 –        else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
467  
468 +        else if (input==4)      /* E */         /* Implementation of the Erbs model. W.Sprenger (04/13) */
469 +                {
470 +                        
471 +                        if (altitude<=0)
472 +                        {
473 +                                if (suppress_warnings==0 && globalirradiance > 50)
474 +                                        fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n");
475 +                                globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0;
476 +                        
477 +                        } else {
478 +                        
479 +                        erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude);
480 +                        
481 +                        if (globalirradiance>erbs_s0)
482 +                        {
483 +                                if (suppress_warnings==0)
484 +                                        fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
485 +                                globalirradiance=erbs_s0*0.999;                
486 +                        }
487 +                        
488 +                        erbs_kt=globalirradiance/erbs_s0;
489 +                        
490 +                        if (erbs_kt<=0.22)      diffuseirradiance=globalirradiance*(1-0.09*erbs_kt);
491 +                        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));
492 +                        else if (erbs_kt<1)     diffuseirradiance=globalirradiance*(0.165);
493 +                        
494 +                        directirradiance=globalirradiance-diffuseirradiance;
495 +                        
496 +                        printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
497 +                        printf("# WARNING: the -E option is only recommended for a rough estimation!\n");
498 +                        
499 +                        directirradiance=directirradiance/sin(altitude);
500 +                                                                                                                        
501 +                        }
502 +                        
503 +                check_irradiances();
504 +                skybrightness = sky_brightness();
505 +                skyclearness =  sky_clearness();
506 +                check_parametrization();
507  
508 +                if (output==0 || output==2)
509 +                        {
510 +                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
511 +                        directilluminance = directirradiance*direct_n_effi_PEREZ();
512 +                        check_illuminances();
513 +                        }
514 +
515 +                }
516 +                
517 +                
518 +                
519          
520 < /* normalization factor for the relative sky luminance distribution, diffuse part*/
520 >        else    { fprintf(stderr,"error at the input arguments"); exit(1); }
521  
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 –        }
522  
523 < /* read the coefficients for the Perez sky luminance model */
524 <        if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
525 <        {
408 <                fprintf(stderr,"lect_coeff_perez does not work\n");
409 <                exit(2);
410 <        }
411 <
523 >        
524 >        /* normalization factor for the relative sky luminance distribution, diffuse part*/
525 >        
526          if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
527          {
528                  fprintf(stderr,"Out of memory in function main");
# Line 416 | Line 530 | computesky()                   /* compute sky parameters */
530          }
531  
532          /* read the angles */
533 <        theta_o = theta_ordered("defangle.dat");
534 <        phi_o = phi_ordered("defangle.dat");
533 >        theta_o = defangle_theta;
534 >        phi_o = defangle_phi;
535 >        
536  
537 < /* parameters for the perez model */
537 >        /* parameters for the perez model */
538          coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
539  
540 < /*calculation of the modelled luminance */
540 >        
541 >        
542 >        /*calculation of the modelled luminance */
543          for (j=0;j<145;j++)
544          {
545                  theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
546 +                                
547                  *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
548 <                /*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
548 >                
549 >                /* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
550          }
551 <
551 >        
552          /* integration of luminance for the normalization factor, diffuse part of the sky*/
553 +        
554          diffnormalization = integ_lv(lv_mod, theta_o);
435 –        /*printf("perez integration %lf\n", diffnormalization);*/
555          
437 –        
556  
557  
558 < /*normalization coefficient in lumen or in watt*/
558 >        /*normalization coefficient in lumen or in watt*/
559          if (output==0)
560                  {
561 <                diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
561 >                diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
562                  }
563          else if (output==1)
564                  {
565 <                diffnormalization = diffusirradiance/diffnormalization;
565 >                diffnormalization = diffuseirradiance/diffnormalization;
566                  }
567          else if (output==2)
568                  {
569 <                diffnormalization = diffusilluminance/diffnormalization;
569 >                diffnormalization = diffuseilluminance/diffnormalization;
570                  }
571  
572 <        else    {fprintf(stderr,"output argument : wrong number"); exit(1);}
572 >        else    {fprintf(stderr,"Wrong output specification.\n"); exit(1);}
573  
574  
575  
576  
577 < /* calculation for the solar source */  
577 >        /* calculation for the solar source */  
578          if (output==0)
579                  solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
580                  
# Line 466 | Line 584 | computesky()                   /* compute sky parameters */
584          else
585                  solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
586          
469 –                        
587  
588  
589 < /* Compute the ground radiance */
590 < zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
591 < zenithbr*=diffnormalization;
475 < fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
589 >        /* Compute the ground radiance */
590 >        zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
591 >        zenithbr*=diffnormalization;
592          
593 < if (skyclearness==1)
593 >        if (skyclearness==1)
594          normfactor = 0.777778;
595                  
596 < if (skyclearness>=6)
596 >        if (skyclearness>=6)
597          {              
598          F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
599          normfactor = normsc()/F2/M_PI;
600          }
601  
602 < if ( (skyclearness>1) && (skyclearness<6) )
602 >        if ( (skyclearness>1) && (skyclearness<6) )
603          {
604          S_INTER=1;
605          F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
606          normfactor = normsc()/F2/M_PI;
607          }
608  
609 < groundbr = zenithbr*normfactor;
494 < printf("# Ground ambient level: %.1f\n", groundbr);
609 >        groundbr = zenithbr*normfactor;
610  
611 < if (dosun&&(skyclearness>1))
612 < groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
611 >        if (dosun&&(skyclearness>1))
612 >        groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
613  
614 < groundbr *= gprefl;
614 >        groundbr *= gprefl;
615  
616  
617 +                
618 +        if(*(c_perez+1)>0)
619 +        {
620 +          if(suppress_warnings==0)
621 +                {  fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));}  
622 +          print_error_sky();
623 +          exit(0);
624 +        }
625  
626 +
627   return;
628   }
629  
# Line 507 | Line 631 | return;
631  
632  
633  
634 + double solar_sunset(int month,int day)
635 + {
636 +     float W;
637 +     extern double s_latitude;
638 +     W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
639 +     return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
640 + }
641  
642  
643 < printsky()                      /* print out sky */
643 >
644 >
645 > double solar_sunrise(int month,int day)
646   {
647 +     float W;
648 +     extern double s_latitude;
649 +     W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
650 +     return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
651 + }
652 +
653 +
654 +
655 +
656 + void printsky()
657 + {      
658 +        
659 +        printf("# Local solar time: %.2f\n", st);
660 +        printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
661 +
662 +
663          if (dosun&&(skyclearness>1))
664 <                {              
664 >        {              
665                  printf("\nvoid light solar\n");
666                  printf("0\n0\n");
667                  printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
668                  printf("\nsolar source sun\n");
669                  printf("0\n0\n");
670                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
671 <                }
523 <                
524 <        if (dosun&&(skyclearness==1))
525 <                {              
671 >        } else if (dosun) {
672                  printf("\nvoid light solar\n");
673                  printf("0\n0\n");
674                  printf("3 0.0 0.0 0.0\n");
675                  printf("\nsolar source sun\n");
676                  printf("0\n0\n");
677                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
678 <                }
678 >        }
679          
680  
681          printf("\nvoid brightfunc skyfunc\n");
682          printf("2 skybright perezlum.cal\n");
683          printf("0\n");
684          printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
685 <            *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
686 <            sundir[0], sundir[1], sundir[2]);
685 >                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
686 >                sundir[0], sundir[1], sundir[2]);
687 >        
688   }
689  
690  
691 < printdefaults()                 /* print default values */
691 >
692 > void print_error_sky()
693   {
694 +
695 +
696 +        sundir[0] = -sin(azimuth)*cos(altitude);
697 +        sundir[1] = -cos(azimuth)*cos(altitude);
698 +        sundir[2] = sin(altitude);
699 +
700 +        printf("# Local solar time: %.2f\n", st);
701 +        printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
702 +
703 +        printf("\nvoid brightfunc skyfunc\n");
704 +        printf("2 skybright perezlum.cal\n");
705 +        printf("0\n");
706 +        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]);
707 + }
708 +        
709 +
710 +
711 +
712 +
713 + void printdefaults()                    /* print default values */
714 + {
715          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
716          if (zenithbr > 0.0)
717                  printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
# Line 554 | Line 723 | printdefaults()                        /* print default values */
723   }
724  
725  
726 < userror(msg)                    /* print usage error and quit */
727 < char  *msg;
726 >
727 >
728 > void usage_error(char* msg)                     /* print usage error and quit */
729   {
730          if (msg != NULL)
731 <                fprintf(stderr, "%s: Use error - %s\n", progname, msg);
732 <        fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
733 <        fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
731 >                fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
732 >        fprintf(stderr, "Usage: %s      month day hour    [...]\n", progname);
733 >        fprintf(stderr, "   or: %s -ang altitude azimuth  [...]\n", progname);
734 >        fprintf(stderr, "               followed by:      -P          epsilon delta [options]\n");
735 >        fprintf(stderr, "                        or:      [-W|-L|-G]  direct_value diffuse_value [options]\n");
736 >        fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
737 >        fprintf(stderr, "       Description:\n");
738          fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
739          fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
740          fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
741          fprintf(stderr, "       -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
742 +        fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
743 +        fprintf(stderr, "       Output specification with option:\n");
744          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");
745 +        fprintf(stderr, "       gendaylit version 2.4 (2013/09/04)  \n\n");
746          exit(1);
747   }
748  
749  
750  
751 < double
752 < normsc()                        /* compute normalization factor (E0*F2/L0) */
751 >
752 > double normsc()           /* compute normalization factor (E0*F2/L0) */
753   {
754          static double  nfc[2][5] = {
755                                  /* clear sky approx. */
# Line 595 | Line 772 | normsc()                       /* compute normalization factor (E0*F2/L0)
772  
773  
774  
775 < printhead(ac, av)               /* print command header */
776 < register int  ac;
777 < register char  **av;
775 >
776 >
777 > void printhead(int ac, char** av)               /* print command header */
778   {
779          putchar('#');
780          while (ac--) {
# Line 610 | Line 787 | register char  **av;
787  
788  
789  
613 – 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 – }
790  
791  
630 –
792   /* Perez models */
793  
794   /* Perez global horizontal luminous efficacy model */
# Line 637 | Line 798 | double glob_h_effi_PEREZ()
798          double  value;
799          double  category_bounds[10], a[10], b[10], c[10], d[10];
800          int     category_total_number, category_number, i;
801 <
802 <
803 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
804 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
805 <
801 >        
802 >        check_parametrization();
803 >        
804 >        
805 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
806 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */
807 >    
808 >    
809          /* initialize category bounds (clearness index bounds) */
810  
811          category_total_number = 8;
# Line 696 | Line 860 | fprintf(stderr, "Warning : skyclearness or skybrightne
860  
861  
862  
699 –
863          for (i=1; i<=category_total_number; i++)
864          {
865                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
# Line 710 | Line 873 | fprintf(stderr, "Warning : skyclearness or skybrightne
873   }
874  
875  
876 +
877 +
878   /* global horizontal diffuse efficacy model, according to PEREZ */
879   double glob_h_diffuse_effi_PEREZ()
880   {
# Line 717 | Line 882 | double glob_h_diffuse_effi_PEREZ()
882          double  category_bounds[10], a[10], b[10], c[10], d[10];
883          int     category_total_number, category_number, i;
884  
885 +        check_parametrization();
886          
887 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
888 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
889 <
887 >        
888 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
889 > fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
890 >    
891   /* initialize category bounds (clearness index bounds) */
892  
893          category_total_number = 8;
894  
895 + //XXX:  category_bounds > 0.1
896          category_bounds[1] = 1;
897          category_bounds[2] = 1.065;
898          category_bounds[3] = 1.230;
# Line 775 | Line 943 | fprintf(stderr, "Warning : skyclearness or skybrightne
943  
944  
945  
946 <
946 >        category_number = -1;
947          for (i=1; i<=category_total_number; i++)
948          {
949                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
950                          category_number = i;
951          }
952  
953 +        if (category_number == -1) {
954 +                if (suppress_warnings==0)
955 +                    fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
956 +                print_error_sky();
957 +                exit(0);
958 +        }
959 +                
960 +
961          value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) +
962              d[category_number]*log(skybrightness);
963  
964          return(value);
965 +
966   }
967  
968  
969 +
970 +
971 +
972 +
973   /* direct normal efficacy model, according to PEREZ */
974  
975   double direct_n_effi_PEREZ()
# Line 799 | Line 980 | double         category_bounds[10], a[10], b[10], c[10], d[10
980   int     category_total_number, category_number, i;
981  
982  
983 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
984 < fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
983 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
984 >   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/
985  
986  
987   /* initialize category bounds (clearness index bounds) */
# Line 874 | Line 1055 | return(value);
1055   /*check the range of epsilon and delta indexes of the perez parametrization*/
1056   void check_parametrization()
1057   {
1058 < if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
1058 >
1059 > if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1060                  {
1061 <                fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
1062 <                exit(1);        
1061 >
1062 > /*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1063 >                
1064 >                if (skyclearness<skyclearinf){
1065 >                        /* if (suppress_warnings==0)
1066 >                            fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1067 >                        skyclearness=skyclearinf;
1068                  }
1069 +                if (skyclearness>skyclearsup){
1070 +                        /* if (suppress_warnings==0)
1071 +                            fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1072 +                        skyclearness=skyclearsup-0.001;
1073 +                }
1074 +                if (skybrightness<skybriginf){
1075 +                        /* if (suppress_warnings==0)
1076 +                            fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1077 +                        skybrightness=skybriginf;
1078 +                }
1079 +                if (skybrightness>skybrigsup){
1080 +                        /* if (suppress_warnings==0)
1081 +                            fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1082 +                        skybrightness=skybrigsup;
1083 +                }
1084 +
1085 +        return; }
1086          else return;
1087   }
1088  
1089  
1090 < /* likelihood of the direct and diffuse components */
1090 >
1091 >
1092 >
1093 > /* validity of the direct and diffuse components */
1094   void    check_illuminances()
1095   {
1096 <        if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
1097 <        {
1098 <        fprintf(stderr,"direct or diffuse illuminances out of range\n");
1099 <        exit(1);
1096 >        if (directilluminance < 0) {
1097 >                if(suppress_warnings==0)
1098 >                { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1099 >                directilluminance = 0.0;
1100          }
1101 < return;
1101 >        if (diffuseilluminance < 0) {
1102 >                if(suppress_warnings==0)
1103 >                { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1104 >                diffuseilluminance = 0.0;
1105 >        }
1106 >        
1107 >        if (directilluminance+diffuseilluminance==0 && altitude > 0) {
1108 >                if(suppress_warnings==0)
1109 >                { fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); }
1110 >                print_error_sky();
1111 >                exit(0);
1112 >        }
1113 >        
1114 >        if (directilluminance > solar_constant_l) {
1115 >                if(suppress_warnings==0)
1116 >                { fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); }
1117 >                print_error_sky();
1118 >                exit(0);
1119 >        }
1120   }
1121  
1122  
1123   void    check_irradiances()
1124   {
1125 <        if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
1126 <        {
1127 <        fprintf(stderr,"direct or diffuse irradiances out of range\n");
1128 <        exit(1);
1129 <        }      
1130 < return;
1125 >        if (directirradiance < 0) {
1126 >                if(suppress_warnings==0)
1127 >                { fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); }
1128 >                directirradiance = 0.0;
1129 >        }
1130 >        if (diffuseirradiance < 0) {
1131 >                if(suppress_warnings==0)
1132 >                { fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); }
1133 >                diffuseirradiance = 0.0;
1134 >        }
1135 >        
1136 >        if (directirradiance+diffuseirradiance==0 && altitude > 0) {
1137 >                if(suppress_warnings==0)
1138 >                { fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); }
1139 >                print_error_sky();
1140 >                exit(0);
1141 >        }
1142 >        
1143 >        if (directirradiance > solar_constant_e) {
1144 >                if(suppress_warnings==0)
1145 >                { fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); }
1146 >                print_error_sky();
1147 >                exit(0);
1148 >        }
1149   }
1150          
1151  
# Line 912 | Line 1155 | double sky_brightness()
1155   {
1156   double value;
1157  
1158 < value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1158 > value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
1159  
1160   return(value);
1161   }
# Line 921 | Line 1164 | return(value);
1164   /* Perez sky's clearness */
1165   double sky_clearness()
1166   {
1167 < double value;
1167 >        double value;
1168  
1169 < 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) ;
1169 >        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) ;
1170  
1171 < return(value);
1171 >        return(value);
1172   }
1173  
1174  
1175  
1176   /* diffus horizontal irradiance from Perez sky's brightness */
1177 < double diffus_irradiance_from_sky_brightness()
1177 > double diffuse_irradiance_from_sky_brightness()
1178   {
1179          double value;
1180  
# Line 946 | Line 1189 | double direct_irradiance_from_sky_clearness()
1189   {
1190          double value;
1191  
1192 <        value = diffus_irradiance_from_sky_brightness();
1192 >        value = diffuse_irradiance_from_sky_brightness();
1193          value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
1194  
1195          return(value);
1196   }
1197  
1198  
1199 < void illu_to_irra_index(void)
1199 >
1200 >
1201 > void illu_to_irra_index()
1202   {
1203 < double  test1=0.1, test2=0.1;
1203 > double  test1=0.1, test2=0.1, d_eff;
1204   int     counter=0;      
1205  
1206 < diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
1207 < directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1206 > diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l);
1207 > directirradiance = directilluminance*solar_constant_e/(solar_constant_l);
1208   skyclearness =  sky_clearness();
1209   skybrightness = sky_brightness();
1210 < 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);*/
1210 > check_parametrization();
1211  
1212 <        test1=diffusirradiance;
1212 >
1213 > while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1214 >                || (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1215 >                || (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1216 >                 && !(counter==9) )
1217 >        {
1218 >        
1219 >        test1=diffuseirradiance;
1220          test2=directirradiance;
1221          counter++;
1222          
1223 <        diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
1224 <        directirradiance = directilluminance/direct_n_effi_PEREZ();
982 <        /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
1223 >        diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1224 >        d_eff = direct_n_effi_PEREZ();
1225          
1226 +        
1227 +        if (d_eff < 0.1)
1228 +                directirradiance = 0;
1229 +        else    
1230 +                directirradiance = directilluminance/d_eff;
1231 +        
1232          skybrightness = sky_brightness();
1233          skyclearness =  sky_clearness();
1234 <        if (skyclearness>12) skyclearness=12;
1235 <        if (skybrightness<0.05) skybrightness=0.01;
988 <
989 <                /*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
990 <
1234 >        check_parametrization();
1235 >                
1236          }
1237  
1238  
1239   return;
1240   }              
1241  
1242 <
998 < int lect_coeff_perez(char *filename,float **coeff_perez)
1242 > static int get_numlin(float epsilon)
1243   {
1244 <        FILE *fcoeff_perez;
1245 <        float temp;
1246 <        int i,j;
1247 <
1248 <        if ((fcoeff_perez = frlibopen(filename)) == NULL)
1249 <        {
1250 <                fprintf(stderr,"file %s cannot be opened\n", filename);
1251 <                return 1; /* il y a un probleme de fichier */
1252 <        }
1253 <        else
1254 <        {
1255 <                /*printf("file %s  open\n", filename);*/
1256 <        }
1257 <        
1258 <        skip_comments(fcoeff_perez);
1015 <
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 */
1244 >        if (epsilon < 1.065)
1245 >                return 0;
1246 >        else if (epsilon < 1.230)
1247 >                return 1;
1248 >        else if (epsilon < 1.500)
1249 >                return 2;
1250 >        else if (epsilon < 1.950)
1251 >                return 3;
1252 >        else if (epsilon < 2.800)
1253 >                return 4;
1254 >        else if (epsilon < 4.500)
1255 >                return 5;
1256 >        else if (epsilon < 6.200)
1257 >                return 6;
1258 >        return 7;
1259   }
1260  
1027 –
1028 –
1261   /* sky luminance perez model */
1262 < double calc_rel_lum_perez(double dzeta,double gamma,double Z,
1031 < double epsilon,double Delta,float *coeff_perez)
1262 > double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1263   {
1264 +                        
1265          float x[5][4];
1266          int i,j,num_lin;
1267          double c_perez[5];
1268  
1269          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1270          {
1271 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1271 >                fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
1272                  exit(1);
1273          }
1274  
# Line 1045 | Line 1277 | double epsilon,double Delta,float *coeff_perez)
1277          {
1278                  if ( Delta < 0.2 ) Delta = 0.2;
1279          }
1280 <
1281 <        if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
1282 <        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
1283 <        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 <
1280 >        
1281 >        
1282 >        num_lin = get_numlin(epsilon);
1283 >        
1284          for (i=0;i<5;i++)
1285                  for (j=0;j<4;j++)
1286                  {
1287                          x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1288 <                        /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1288 >                        /* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1289                  }
1290  
1291  
# Line 1086 | Line 1312 | double epsilon,double Delta,float *coeff_perez)
1312  
1313  
1314   /* coefficients for the sky luminance perez model */
1315 < void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
1315 > void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[])
1316   {
1317          float x[5][4];
1318          int i,j,num_lin;
1319  
1320          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1321          {
1322 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1322 >                fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n");
1323                  exit(1);
1324          }
1325  
# Line 1102 | Line 1328 | void coeff_lum_perez(double Z, double epsilon, double
1328          {
1329                  if ( Delta < 0.2 ) Delta = 0.2;
1330          }
1331 +        
1332 +        
1333 +        num_lin = get_numlin(epsilon);
1334  
1335 <        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;
1335 >        /*fprintf(stderr,"numlin %d\n", num_lin);*/
1336  
1337          for (i=0;i<5;i++)
1338                  for (j=0;j<4;j++)
# Line 1142 | Line 1364 | void coeff_lum_perez(double Z, double epsilon, double
1364   }
1365  
1366  
1367 +
1368   /* degrees into radians */
1369   double radians(double degres)
1370   {
1371          return degres*M_PI/180.0;
1372   }
1373  
1374 +
1375   /* radian into degrees */
1376   double degres(double radians)
1377   {
1378          return radians/M_PI*180.0;
1379   }
1380  
1381 +
1382   /* calculation of the angles dzeta and gamma */
1383   void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1384   {
# Line 1163 | Line 1388 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1388          else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
1389          {
1390                  printf("error in calculation of gamma (angle between point and sun");
1391 <                exit(3);
1391 >                exit(1);
1392          }
1393          else
1394                  *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
1395   }
1396  
1397  
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;
1398  
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 –        
1199 –        skip_comments(file_in);
1200 –
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 –        
1245 –        skip_comments(file_in);
1246 –
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 – /********************************************************************************/
1399   double integ_lv(float *lv,float *theta)
1400   {
1401          int i;
1402          double buffer=0.0;
1403 <
1403 >        
1404          for (i=0;i<145;i++)
1405 +        {
1406                  buffer += (*(lv+i))*cos(radians(*(theta+i)));
1407 <
1407 >        }
1408 >                        
1409          return buffer*2*M_PI/144;
1289 –
1410   }
1411  
1412  
1413  
1294 –
1295 –
1296 –
1414   /* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1415  
1416   double get_eccentricity()
# Line 1306 | Line 1423 | double get_eccentricity()
1423              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1424  
1425          return (E0);
1309 –
1426   }
1427  
1428  
# Line 1314 | Line 1430 | double get_eccentricity()
1430   double  air_mass()
1431   {
1432   double  m;
1317 –
1433   if (sunzenith>90)
1434          {
1435 <        fprintf(stderr, "solar zenith angle larger than 90� in fuction air_mass():\n the models used are not more valid\n");
1436 <        exit(1);
1435 >        if(suppress_warnings==0)
1436 >        { fprintf(stderr, "Warning: air mass has reached the maximal value\n"); }
1437 >        sunzenith=90;
1438          }
1323 –        
1439   m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1440   return(m);
1441   }
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 –
1442  
1443  
1444  

Diff Legend

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