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.11 by greg, Tue Apr 30 17:05:27 2013 UTC vs.
Revision 2.13 by greg, Wed Aug 14 17:11:43 2013 UTC

# Line 1 | Line 1
1 – #ifndef lint
2 – static const char RCSid[] = "$Id$";
3 – #endif
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
10 – *
11 – *      24.1.2006                some adjustments for cygwin compilation, inclusion of RADIANCE3.7 libraries, by J. Wienold
12 – *               2011/10/08      [email protected]:
13 – *                                               - integrated coeff_perez.dat and defangles.dat
14 – *                                               - avoid some segfaults caused by out of range parameters and
15 – *                                               - numerically dangerous range checks
7   */
8  
18 – /*
19 – *  gendaylit.c         program to generate the angular distribution of the daylight.
20 – *                      Our zenith is along the Z-axis, the X-axis
21 – *                      points east, and the Y-axis points north.
22 – */
23 –
9   #define _USE_MATH_DEFINES
25 –
10   #include  <stdio.h>
11   #include  <string.h>
12   #include  <math.h>
13   #include  <stdlib.h>
14  
15   #include  "color.h"
16 + #include  "sun.h"
17   #include  "paths.h"
18  
19   #define  DOT(v1,v2)     (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
# Line 38 | Line 23 | double  normsc();
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,0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,
27 <        -1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156,6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285,
28 <        -1.1000,-0.2515,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,-2.6204,-0.0156,0.1597,0.4199,-0.5562,
29 <        -0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297,33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329,
30 <        -0.6000,-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,1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294,
31 <        -1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750,14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844,
32 <        -1.0000,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,0.2656,1.0468,-0.3788,-2.4517,1.4656,
33 <        -1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625,-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636};
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 < float defangle_theta[] = {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, 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, 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, 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, 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, 24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};
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[] = {0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
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  
# Line 62 | Line 64 | double sky_clearness();
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  
66 – /* 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();
# Line 74 | Line 78 | void   check_illuminances();
78   void    illu_to_irra_index();
79   void    print_error_sky();
80  
81 <
78 < /* Perez sky luminance model */
79 < double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
80 <                double epsilon,double Delta,float coeff_perez[]);
81 < /* coefficients for the sky luminance perez model */
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);
# Line 86 | Line 86 | void   theta_phi_to_dzeta_gamma(double theta,double phi,
86   double  integ_lv(float *lv,float *theta);
87  
88   void printdefaults();
89 + void check_sun_position();
90   void computesky();
91   void printhead(int ac, char** av);
92   void userror(char* msg);
# Line 97 | Line 98 | FILE * frlibopen(char* fname);
98   double  get_eccentricity();
99   double  air_mass();
100  
101 < extern int jdate(int month, int day);
102 < extern double stadj(int  jd);
103 < extern double sdec(int  jd);
104 < extern double salt(double sd, double st);
104 < extern double sazi(double sd, double st);
101 > double  solar_sunset();
102 > double  solar_sunrise();
103 > double  stadj();
104 > int     jdate(int month, int day);
105  
106  
107 +
108   /* sun calculation constants */
109 < extern double  s_latitude;
110 < extern double  s_longitude;
111 < extern double  s_meridian;
109 > extern double   s_latitude;
110 > extern double   s_longitude;
111 > extern double   s_meridian;
112  
113   const double    AU = 149597890E3;
114   const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
# Line 116 | Line 117 | const double   solar_constant_l = 127.5;   /* solar co
117   const double    half_sun_angle = 0.2665;
118   const double    half_direct_angle = 2.85;
119  
120 < const double    skyclearinf = 1.000;    /* limitations for the variation of the Perez parameters */
120 > const double    skyclearinf = 1.0;      /* limitations for the variation of the Perez parameters */
121   const double    skyclearsup = 12.1;
122   const double    skybriginf = 0.01;
123   const double    skybrigsup = 0.6;
# Line 134 | Line 135 | double  altitude, azimuth;                     /* or solar angles */
135   /* definition of the sky conditions through the Perez parametrization */
136   double  skyclearness = 0;
137   double  skybrightness = 0;
138 < double  solarradiance;  /*radiance of the sun disk and of the circumsolar area*/
139 < double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance;
140 < double  sunzenith, daynumber=150, atm_preci_water=2;
138 > double  solarradiance;
139 > double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
140 > double  sunzenith, daynumber, atm_preci_water=2;
141  
142 < double sunaltitude_border = 0;
142 > /*double  sunaltitude_border = 0;*/
143   double  diffnormalization = 0;
144 < double dirnormalization = 0;
144 > double  dirnormalization = 0;
145   double  *c_perez;
146  
147 < int     output=0;       /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
148 < int     input=0;        /*define the input for the calulation*/
147 > int     output=0;       /* define the unit of the output (sky luminance or radiance): */
148 >                        /* visible watt=0, solar watt=1, lumen=2 */
149 > int     input=0;        /* define the input for the calulation */
150  
151   int     suppress_warnings=0;
152  
153          /* default values */
154 < int  cloudy = 0;                                /* 1=standard, 2=uniform */
155 < int  dosun = 1;
154 > int     cloudy = 0;                             /* 1=standard, 2=uniform */
155 > int     dosun = 1;
156   double  zenithbr = -1.0;
157   double  betaturbidity = 0.1;
158   double  gprefl = 0.2;
# Line 162 | Line 164 | double  groundbr = 0;
164   double  F2;
165   double  solarbr = 0.0;
166   int     u_solar = 0;                            /* -1=irradiance, 1=radiance */
167 + float   timeinterval = 0;
168  
169 < char  *progname;
170 < char  errmsg[128];
169 > char    *progname;
170 > char    errmsg[128];
171  
172  
173 +
174 +
175   int main(int argc, char** argv)
176   {
177          int  i;
# Line 201 | Line 206 | int main(int argc, char** argv)
206                                  cloudy = 0;
207                                  dosun = argv[i][0] == '+';
208                                  break;
204 –                        case 'r':
209                          case 'R':
210                                  u_solar = argv[i][1] == 'R' ? -1 : 1;
211                                  solarbr = atof(argv[++i]);
# Line 233 | Line 237 | int main(int argc, char** argv)
237                                  break;
238                          
239                          case 'O':
240 <                                output = atoi(argv[++i]);       /*define the unit of the output of the program :
240 >                                output = atof(argv[++i]);       /*define the unit of the output of the program :
241                                                                  sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
242                                  break;
243                                  
# Line 261 | Line 265 | int main(int argc, char** argv)
265                                  diffuseirradiance = atof(argv[++i]);
266                                  break;
267                          
268 <                        case 'l':
268 >                        case 'E':                                       /* Erbs model based on the */
269 >                                input = 4;                              /* global-horizontal irradiance [W/m^2] */
270 >                                globalirradiance = atof(argv[++i]);
271 >                                break;
272 >                        
273 >                        /*
274 >                        case 'l':
275                                  sunaltitude_border = atof(argv[++i]);
276                                  break;
277 +                        */
278 +                                
279 +                        case 'i':
280 +                                timeinterval = atof(argv[++i]);
281 +                                break;
282                          
283                          
284                          default:
# Line 279 | Line 294 | int main(int argc, char** argv)
294                      progname, (s_longitude-s_meridian)*12/M_PI);
295  
296  
297 <        /* allocation dynamique de memoire pour les pointeurs */
297 >        /* dynamic memory allocation for the pointers */
298 >        
299          if ( (c_perez = calloc(5, sizeof(double))) == NULL )
300          {
301 <                fprintf(stderr,"Out of memory error in function main !");
301 >                fprintf(stderr,"Out of memory error in function main");
302                  return 1;
303          }
304  
305          printhead(argc, argv);
290 –
306          computesky();
307 +        
308 +        if(*(c_perez+1)>0)
309 +        {
310 +          fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));      
311 +          print_error_sky();
312 +          exit(1);
313 +        }
314 +        
315          printsky();
293 –
316          return 0;
317   }
318  
319  
320 < void computesky()                       /* compute sky parameters */
320 >
321 >
322 >
323 >
324 >
325 >
326 > void computesky()
327   {
328  
301 –        /* new variables */
329          int     j;
330 <        float   *lv_mod;  /* 145 luminance values*/
331 <          /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
330 >        
331 >        float   *lv_mod;  /* 145 luminance values */
332          float   *theta_o, *phi_o;
333          double  dzeta, gamma;
334          double  normfactor;
335 +        double  erbs_s0, erbs_kt;
336  
337  
310 –
338          /* compute solar direction */
339 <
339 >                
340          if (month) {                    /* from date and time */
341                  int  jd;
342                  double  sd, st;
# Line 320 | Line 347 | void computesky()                      /* compute sky parameters */
347                          st = hour;
348                  else
349                          st = hour + stadj(jd);
350 +                                
351 +                
352 +                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
353 +                print_error_sky();
354 +                exit(1);
355 +                }
356 +                
357 +                                        
358 +                if(timeinterval) {
359 +                        
360 +                        if(timeinterval<0) {
361 +                        fprintf(stderr, "time interval negative\n");
362 +                        exit(1);
363 +                        }
364 +                                        
365 +                        if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) {                
366 +                        
367 +                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
368 +                        st= (st+timeinterval/120+solar_sunrise(month,day))/2;
369 +                        }
370 +                
371 +                        if(fabs(solar_sunset(month,day)-st)<timeinterval/60) {
372 +                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
373 +                        st= (st-timeinterval/120+solar_sunset(month,day))/2;
374 +                        }
375 +                }
376 +                
377 +                
378                  altitude = salt(sd, st);
379                  azimuth = sazi(sd, st);
380                  
381                  daynumber = (double)jdate(month, day);
382 <
382 >                
383          }
384          
385          
386 <        
387 <        
388 <        
389 <        /* if loop for the -l option. 01/2013 Sprenger  */
335 <        
386 >
387 >
388 >        /* if loop for the -l option. W.Sprenger (01/2013) */
389 >        /*
390          if (altitude*180/M_PI < sunaltitude_border) {
391          
392 <        if (suppress_warnings==0)
393 <            fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
392 >        if (suppress_warnings==0) {
393 >            fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
394 >                                   }
395          print_error_sky();
396 <        exit(0);
396 >        exit(1);
397          }
398 <        
399 <        
398 >        */
399 >
400                          
346 –        
347 –                        
401          if (!cloudy && altitude > 87.*M_PI/180.) {
402                  
403                  if (suppress_warnings==0) {
# Line 366 | Line 419 | void computesky()                      /* compute sky parameters */
419          
420  
421          /* compute the inputs for the calculation of the light distribution over the sky*/
422 <        if (input==0)
422 >        if (input==0)           /* P */
423                  {
424                  check_parametrization();
425                  diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
# Line 382 | Line 435 | void computesky()                      /* compute sky parameters */
435                  }
436          
437  
438 <        else if (input==1)
438 >        else if (input==1)      /* W */
439                  {
440                  check_irradiances();
441                  skybrightness = sky_brightness();
# Line 399 | Line 452 | void computesky()                      /* compute sky parameters */
452                  }
453                          
454          
455 <        else if (input==2)
455 >        else if (input==2)      /* L */
456                  {              
457                  check_illuminances();
458                  illu_to_irra_index();
# Line 407 | Line 460 | void computesky()                      /* compute sky parameters */
460                  }
461                  
462  
463 <        else if (input==3)
463 >        else if (input==3)      /* G */
464                  {
465                          if (altitude<=0)
466                          {
# Line 416 | Line 469 | void computesky()                      /* compute sky parameters */
469                                  directirradiance = 0;
470                                  diffuseirradiance = 0;
471                          } else {
472 <                                directirradiance=directirradiance/sin(altitude);
472 >                        
473 >                                directirradiance=directirradiance/sin(altitude);
474                          }
475 +                                
476                  check_irradiances();
477                  skybrightness = sky_brightness();
478                  skyclearness =  sky_clearness();
# Line 432 | Line 487 | void computesky()                      /* compute sky parameters */
487  
488                  }
489  
490 +
491 +        else if (input==4)      /* E */         /* Implementation of the Erbs model. W.Sprenger (04/13) */
492 +                {
493 +                        
494 +                        if (altitude<=0)
495 +                        {
496 +                                if (suppress_warnings==0 && globalirradiance > 50)
497 +                                        fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n");
498 +                                globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0;
499 +                        
500 +                        } else {
501 +                        
502 +                        erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude);
503 +                        
504 +                        if (globalirradiance>erbs_s0)
505 +                        {
506 +                                if (suppress_warnings==0)
507 +                                        fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
508 +                                
509 +                                globalirradiance=erbs_s0*0.999;                
510 +                        }
511 +                        
512 +                        erbs_kt=globalirradiance/erbs_s0;
513 +                        
514 +                        if (erbs_kt<=0.22)      diffuseirradiance=globalirradiance*(1-0.09*erbs_kt);
515 +                        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));
516 +                        else if (erbs_kt<1)     diffuseirradiance=globalirradiance*(0.165);
517 +                        
518 +                        directirradiance=globalirradiance-diffuseirradiance;
519 +                        
520 +                        printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
521 +                        printf("# WARNING: the -E option is only recommended for a rough estimation!");
522 +                        
523 +                        directirradiance=directirradiance/sin(altitude);
524 +                                                                                                                        
525 +                        }
526 +                        
527 +                check_irradiances();
528 +                skybrightness = sky_brightness();
529 +                skyclearness =  sky_clearness();
530 +                check_parametrization();
531 +
532 +                if (output==0 || output==2)
533 +                        {
534 +                        diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
535 +                        directilluminance = directirradiance*direct_n_effi_PEREZ();
536 +                        check_illuminances();
537 +                        }
538 +
539 +                }
540 +                
541 +                
542 +                
543          
544 <        else    {fprintf(stderr,"error in giving the input arguments"); exit(1);}
544 >        else    {fprintf(stderr,"error at the input arguments"); exit(1);}
545  
546  
547          
548          /* normalization factor for the relative sky luminance distribution, diffuse part*/
549 <
549 >        
550          if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
551          {
552                  fprintf(stderr,"Out of memory in function main");
# Line 448 | Line 556 | void computesky()                      /* compute sky parameters */
556          /* read the angles */
557          theta_o = defangle_theta;
558          phi_o = defangle_phi;
559 +        
560  
561          /* parameters for the perez model */
562          coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
563  
564 +        
565 +        
566 +        
567 +        
568          /*calculation of the modelled luminance */
569          for (j=0;j<145;j++)
570          {
571                  theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
572 +                                
573                  *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
574 <                // printf("theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));
574 >                
575 >                /* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */
576          }
577 <
577 >        
578          /* integration of luminance for the normalization factor, diffuse part of the sky*/
579 +        
580          diffnormalization = integ_lv(lv_mod, theta_o);
465 –        /*printf("perez integration %lf\n", diffnormalization);*/
581          
467 –        
582  
583  
584          /*normalization coefficient in lumen or in watt*/
# Line 548 | Line 662 | void print_error_sky()
662          
663  
664  
665 +
666 +
667 +
668 +
669 + double solar_sunset(int month,int day)
670 + {
671 +     float W;
672 +     extern double s_latitude;
673 +     W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
674 +     return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
675 + }
676 +
677 +
678 + double solar_sunrise(int month,int day)
679 + {
680 +     float W;
681 +     extern double s_latitude;
682 +     W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day))));
683 +     return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15));
684 + }
685 +
686 +
687 +
688 +
689 +
690 +
691 +
692 +
693 +
694 +
695 +
696   void printsky()                 /* print out sky */
697   {
698          if (dosun&&(skyclearness>1))
# Line 567 | Line 712 | void printsky()                        /* print out sky */
712                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
713          }
714          
715 +
716          printf("\nvoid brightfunc skyfunc\n");
717          printf("2 skybright perezlum.cal\n");
718          printf("0\n");
719          printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
720                  *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
721                  sundir[0], sundir[1], sundir[2]);
722 +        
723   }
724  
725  
# Line 592 | Line 739 | void printdefaults()                   /* print default values */
739   void userror(char* msg)                 /* print usage error and quit */
740   {
741          if (msg != NULL)
742 <                fprintf(stderr, "%s: Use error - %s\n", progname, msg);
743 <        fprintf(stderr, "Usage: %s month day hour [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname);
744 <        fprintf(stderr, "or:    %s -ang altitude azimuth [-P|-W|-L|-G] direct_value diffuse_value [options]\n", progname);
742 >                fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
743 >        fprintf(stderr, "Usage: %s      month day hour    [...]\n", progname);
744 >        fprintf(stderr, "   or: %s -ang altitude azimuth  [...]\n", progname);
745 >        fprintf(stderr, "               followed by:      -P          epsilon delta [options]\n");
746 >        fprintf(stderr, "                        or:      [-W|-L|-G]  direct_value diffuse_value [options]\n");
747 >        fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
748 >        fprintf(stderr, "       Description:\n");
749          fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
750          fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
751          fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
752          fprintf(stderr, "       -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
753 +        fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
754 +        fprintf(stderr, "       Output specification with option:\n");
755          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");
756 <        fprintf(stderr, "       gendaylit version 2.00 (2013/01/28)  \n");
756 >        fprintf(stderr, "       gendaylit version 2.3 (2013/08/08)  \n\n");
757          exit(1);
758   }
759  
# Line 642 | Line 795 | void printhead(int ac, char** av)              /* print command he
795  
796  
797  
645 –
646 –
798   /* Perez models */
799  
800   /* Perez global horizontal luminous efficacy model */
# Line 653 | Line 804 | double glob_h_effi_PEREZ()
804          double  value;
805          double  category_bounds[10], a[10], b[10], c[10], d[10];
806          int     category_total_number, category_number, i;
807 <
808 <
809 < if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
810 <     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n");
811 <
807 >        
808 >        check_parametrization();
809 >        
810 >        
811 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
812 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */
813 >    
814 >    
815          /* initialize category bounds (clearness index bounds) */
816  
817          category_total_number = 8;
# Line 733 | Line 887 | double glob_h_diffuse_effi_PEREZ()
887          double  category_bounds[10], a[10], b[10], c[10], d[10];
888          int     category_total_number, category_number, i;
889  
890 +
891 +
892          
893 < if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
894 <     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_effi_PEREZ \n");
893 >        check_parametrization();
894 >        
895 >        
896 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
897 >     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
898 >    
899 >
900  
901   /* initialize category bounds (clearness index bounds) */
902  
903          category_total_number = 8;
904  
905 < //XXX: category_bounds > 0.1
905 > //XXX:  category_bounds > 0.1
906          category_bounds[1] = 1;
907          category_bounds[2] = 1.065;
908          category_bounds[3] = 1.230;
# Line 802 | Line 963 | if ((skyclearness<skyclearinf || skyclearness>skyclear
963  
964          if (category_number == -1) {
965                  if (suppress_warnings==0)
966 <                    fprintf(stderr, "ERROR: Model parameters out of range\n");
966 >                    fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness);
967                  print_error_sky();
968                  exit(1);
969          }
# Line 812 | Line 973 | if ((skyclearness<skyclearinf || skyclearness>skyclear
973              d[category_number]*log(skybrightness);
974  
975          return(value);
976 +
977   }
978  
979  
980 +
981   /* direct normal efficacy model, according to PEREZ */
982  
983   double direct_n_effi_PEREZ()
# Line 905 | Line 1068 | if (skyclearness<skyclearinf || skyclearness>skyclears
1068  
1069   /*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1070                  if (skyclearness<skyclearinf){
908 –                        skyclearness=skyclearinf;
1071                          if (suppress_warnings==0)
1072 <                            fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness);
1072 >                            /* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1073 >                        skyclearness=skyclearinf;
1074                  }
1075                  if (skyclearness>skyclearsup){
913 –                        skyclearness=skyclearsup-0.1;
1076                          if (suppress_warnings==0)
1077 <                            fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness);
1077 >                            /* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1078 >                        skyclearness=skyclearsup-0.1;
1079                  }
1080                  if (skybrightness<skybriginf){
918 –                        skybrightness=skybriginf;
1081                          if (suppress_warnings==0)
1082 <                            fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness);
1082 >                            /* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1083 >                        skybrightness=skybriginf;
1084                  }
1085                  if (skybrightness>skybrigsup){
923 –                        skybrightness=skybrigsup;
1086                          if (suppress_warnings==0)
1087 <                            fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness);
1087 >                            /* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1088 >                        skybrightness=skybrigsup;
1089                  }
1090  
1091          return; }
# Line 1012 | Line 1175 | double direct_irradiance_from_sky_clearness()
1175   }
1176  
1177  
1178 +
1179 +
1180   void illu_to_irra_index()
1181   {
1182   double  test1=0.1, test2=0.1, d_eff;
# Line 1022 | Line 1187 | directirradiance = directilluminance*solar_constant_e/
1187   skyclearness =  sky_clearness();
1188   skybrightness = sky_brightness();
1189   check_parametrization();
1190 <        
1190 >
1191 >
1192   while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
1193 <                || skyclearness>skyclearinf || skyclearness<skyclearsup
1194 <                || skybrightness>skybriginf || skybrightness<skybrigsup )
1195 <                 && !(counter==5) )
1193 >                || (!(skyclearness<skyclearinf || skyclearness>skyclearsup))
1194 >                || (!(skybrightness<skybriginf || skybrightness>skybrigsup)) )
1195 >                 && !(counter==9) )
1196          {
1197 <
1197 >        
1198          test1=diffuseirradiance;
1199          test2=directirradiance;
1200          counter++;
1201          
1202          diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ();
1203          d_eff = direct_n_effi_PEREZ();
1204 +        
1205 +        
1206          if (d_eff < 0.1)
1207                  directirradiance = 0;
1208 <        else
1208 >        else    
1209                  directirradiance = directilluminance/d_eff;
1210          
1211          skybrightness = sky_brightness();
1212          skyclearness =  sky_clearness();
1213          check_parametrization();
1214          
1215 +        /*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/
1216 +                
1217          }
1218  
1219  
# Line 1072 | Line 1242 | static int get_numlin(float epsilon)
1242   /* sky luminance perez model */
1243   double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
1244   {
1245 +                        
1246          float x[5][4];
1247          int i,j,num_lin;
1248          double c_perez[5];
# Line 1087 | Line 1258 | double calc_rel_lum_perez(double dzeta,double gamma,do
1258          {
1259                  if ( Delta < 0.2 ) Delta = 0.2;
1260          }
1261 <
1261 >        
1262 >        
1263          num_lin = get_numlin(epsilon);
1264 <
1264 >        
1265          for (i=0;i<5;i++)
1266                  for (j=0;j<4;j++)
1267                  {
1268                          x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
1269 <                        /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
1269 >                        /* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
1270                  }
1271  
1272  
# Line 1137 | Line 1309 | void coeff_lum_perez(double Z, double epsilon, double
1309          {
1310                  if ( Delta < 0.2 ) Delta = 0.2;
1311          }
1312 <
1312 >        
1313 >        
1314          num_lin = get_numlin(epsilon);
1315  
1316 <        //fprintf(stderr,"numlin %d\n", num_lin);
1316 >        /*fprintf(stderr,"numlin %d\n", num_lin);*/
1317  
1318          for (i=0;i<5;i++)
1319                  for (j=0;j<4;j++)
# Line 1172 | Line 1345 | void coeff_lum_perez(double Z, double epsilon, double
1345   }
1346  
1347  
1348 +
1349   /* degrees into radians */
1350   double radians(double degres)
1351   {
1352          return degres*M_PI/180.0;
1353   }
1354  
1355 +
1356   /* radian into degrees */
1357   double degres(double radians)
1358   {
1359          return radians/M_PI*180.0;
1360   }
1361  
1362 +
1363   /* calculation of the angles dzeta and gamma */
1364   void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
1365   {
# Line 1201 | Line 1377 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1377  
1378  
1379  
1204 – /********************************************************************************/
1205 – /*      Fonction: integ_lv                                                      */
1206 – /*                                                                              */
1207 – /*      In: float *lv,*theta                                                    */
1208 – /*          int sun_pos                                                         */
1209 – /*                                                                              */
1210 – /*      Out: double                                                             */
1211 – /*                                                                              */
1212 – /*      Update: 29/08/93                                                        */
1213 – /*                                                                              */
1214 – /*      Rem:                                                                    */
1215 – /*                                                                              */
1216 – /*      But: calcul l'integrale de luminance relative sans la dir. du soleil    */
1217 – /*                                                                              */
1218 – /********************************************************************************/
1380   double integ_lv(float *lv,float *theta)
1381   {
1382          int i;
1383          double buffer=0.0;
1384 <
1384 >        
1385          for (i=0;i<145;i++)
1386 +        {
1387                  buffer += (*(lv+i))*cos(radians(*(theta+i)));
1388 <
1388 >        }
1389 >                        
1390          return buffer*2*M_PI/144;
1228 –
1391   }
1392  
1393  
1394  
1233 –
1234 –
1235 –
1395   /* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
1396  
1397   double get_eccentricity()
# Line 1245 | Line 1404 | double get_eccentricity()
1404              0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
1405  
1406          return (E0);
1248 –
1407   }
1408  
1409  
# Line 1253 | Line 1411 | double get_eccentricity()
1411   double  air_mass()
1412   {
1413   double  m;
1256 –
1414   if (sunzenith>90)
1415          {
1416          fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1417          exit(1);
1418          }
1262 –        
1419   m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1420   return(m);
1421   }

Diff Legend

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