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

Comparing src/gen/gendaylit.c (file contents):
Revision 2.12 by greg, Fri Aug 9 16:44:19 2013 UTC vs.
Revision 2.18 by greg, Thu Nov 7 23:15:06 2019 UTC

# Line 4 | Line 4
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 + *  print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018
8   */
9  
10 + #define  _USE_MATH_DEFINES
11   #include  <stdio.h>
12   #include  <string.h>
13   #include  <math.h>
# Line 52 | Line 54 | float defangle_phi[] = {
54          90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
55          80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
56          240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
57 + /* default values for Berlin */
58 + float   locus[] = {
59 + -4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786};
60  
61  
62  
63 < /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
63 > /* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
64   double sky_brightness();
65   double sky_clearness();
66  
# Line 88 | Line 93 | void printdefaults();
93   void check_sun_position();
94   void computesky();
95   void printhead(int ac, char** av);
96 < void userror(char* msg);
96 > void usage_error(char* msg);
97   void printsky();
98  
99   FILE * frlibopen(char* fname);
# Line 97 | Line 102 | FILE * frlibopen(char* fname);
102   double  get_eccentricity();
103   double  air_mass();
104  
105 < double  solar_sunset();
106 < double  solar_sunrise();
102 < double  stadj();
103 < int     jdate(int month, int day);
105 > double  solar_sunset(int month, int day);
106 > double  solar_sunrise(int month, int day);
107  
105 –
106 –
107 – /* sun calculation constants */
108 – extern double   s_latitude;
109 – extern double   s_longitude;
110 – extern double   s_meridian;
111 –
108   const double    AU = 149597890E3;
109   const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
110 < const double    solar_constant_l = 127.5;   /* solar constant klux */
110 > const double    solar_constant_l = 127500;   /* solar constant lux */
111  
112   const double    half_sun_angle = 0.2665;
113   const double    half_direct_angle = 2.85;
114  
115 < const double    skyclearinf = 1.0;      /* limitations for the variation of the Perez parameters */
116 < const double    skyclearsup = 12.1;
115 > const double    skyclearinf = 1.0;          /* limitations for the variation of the Perez parameters */
116 > const double    skyclearsup = 12.01;
117   const double    skybriginf = 0.01;
118   const double    skybrigsup = 0.6;
119  
120  
121  
122   /* required values */
123 + int  year = 0;                                  /* year (optional) */
124   int     month, day;                             /* date */
125   double  hour;                                   /* time */
126   int     tsolar;                                 /* 0=standard, 1=solar */
# Line 146 | Line 143 | double         *c_perez;
143   int     output=0;       /* define the unit of the output (sky luminance or radiance): */
144                          /* visible watt=0, solar watt=1, lumen=2 */
145   int     input=0;        /* define the input for the calulation */
146 <
146 > int     color_output=0;
147   int     suppress_warnings=0;
148  
149          /* default values */
# Line 157 | Line 154 | double  betaturbidity = 0.1;
154   double  gprefl = 0.2;
155   int     S_INTER=0;
156  
157 +
158          /* computed values */
159   double  sundir[3];
160   double  groundbr = 0;
# Line 168 | Line 166 | float  timeinterval = 0;
166   char    *progname;
167   char    errmsg[128];
168  
169 + double  st;
170  
171  
173 –
172   int main(int argc, char** argv)
173   {
174          int  i;
# Line 181 | Line 179 | int main(int argc, char** argv)
179                  return 0;
180          }
181          if (argc < 4)
182 <                userror("arg count");
182 >                usage_error("arg count");
183          if (!strcmp(argv[1], "-ang")) {
184                  altitude = atof(argv[2]) * (M_PI/180);
185                  azimuth = atof(argv[3]) * (M_PI/180);
# Line 189 | Line 187 | int main(int argc, char** argv)
187          } else {
188                  month = atoi(argv[1]);
189                  if (month < 1 || month > 12)
190 <                        userror("bad month");
190 >                        usage_error("bad month");
191                  day = atoi(argv[2]);
192                  if (day < 1 || day > 31)
193 <                        userror("bad day");
193 >                        usage_error("bad day");
194                  hour = atof(argv[3]);
195                  if (hour < 0 || hour >= 24)
196 <                        userror("bad hour");
196 >                        usage_error("bad hour");
197                  tsolar = argv[3][0] == '+';
198          }
199          for (i = 4; i < argc; i++)
# Line 205 | Line 203 | int main(int argc, char** argv)
203                                  cloudy = 0;
204                                  dosun = argv[i][0] == '+';
205                                  break;
206 +                        case 'y':
207 +                                year = atoi(argv[++i]);
208 +                                break;
209                          case 'R':
210                                  u_solar = argv[i][1] == 'R' ? -1 : 1;
211                                  solarbr = atof(argv[++i]);
# Line 213 | Line 214 | int main(int argc, char** argv)
214                                  cloudy = argv[i][0] == '+' ? 2 : 1;
215                                  dosun = 0;
216                                  break;
217 +                        case 'C':
218 +                                if (argv[i][2] == 'I' && argv[i][3] == 'E' ) {
219 +                                locus[0] = -4.607e9;
220 +                                locus[1] = 2.9678e6;
221 +                                locus[2] = 0.09911e3;
222 +                                locus[3] = 0.244063;
223 +                                locus[4] = -2.0064e9;
224 +                                locus[5] = 1.9018e6;
225 +                                locus[6] = 0.24748e3;
226 +                                locus[7] = 0.23704;
227 +                                locus[8] = -3.0;
228 +                                locus[9] = 2.87;
229 +                                locus[10] = -0.275;
230 +                                 }else{ color_output = 1;
231 +                                 }
232 +                                break;
233 +                        case 'l':
234 +                                locus[0] = atof(argv[++i]);
235 +                                locus[1] = atof(argv[++i]);
236 +                                locus[2] = atof(argv[++i]);
237 +                                locus[3] = atof(argv[++i]);
238 +                                locus[4] = locus[0];
239 +                                locus[5] = locus[1];
240 +                                locus[6] = locus[2];
241 +                                locus[7] = locus[3];
242 +                                locus[8] = atof(argv[++i]);
243 +                                locus[9] = atof(argv[++i]);
244 +                                locus[10] = atof(argv[++i]);
245 +                                break;
246 +
247                          case 't':
248                                  betaturbidity = atof(argv[++i]);
249                                  break;
# Line 236 | Line 267 | int main(int argc, char** argv)
267                                  break;
268                          
269                          case 'O':
270 <                                output = atof(argv[++i]);       /*define the unit of the output of the program :
271 <                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
270 >                                output = atof(argv[++i]);       /*define the unit of the output of the program:
271 >                                                                sky and sun luminance/radiance
272 >                                                                (0==W visible, 1==W solar radiation, 2==lm) */
273                                  break;
274                                  
275                          case 'P':
# Line 269 | Line 301 | int main(int argc, char** argv)
301                                  globalirradiance = atof(argv[++i]);
302                                  break;
303                          
272 –                        /*
273 –                        case 'l':
274 –                                sunaltitude_border = atof(argv[++i]);
275 –                                break;
276 –                        */
277 –                                
304                          case 'i':
305                                  timeinterval = atof(argv[++i]);
306                                  break;
# Line 282 | Line 308 | int main(int argc, char** argv)
308                          
309                          default:
310                                  sprintf(errmsg, "unknown option: %s", argv[i]);
311 <                                userror(errmsg);
311 >                                usage_error(errmsg);
312                          }
313                  else
314 <                        userror("bad option");
314 >                        usage_error("bad option");
315  
316 <        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
317 <                fprintf(stderr,
292 <                    "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
316 >        if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180)
317 >                fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
318                      progname, (s_longitude-s_meridian)*12/M_PI);
319  
320  
321          /* dynamic memory allocation for the pointers */
297 –        
322          if ( (c_perez = calloc(5, sizeof(double))) == NULL )
323 <        {
300 <                fprintf(stderr,"Out of memory error in function main");
301 <                return 1;
302 <        }
323 >        { fprintf(stderr,"Out of memory error in function main"); return 1; }
324  
325 +        
326          printhead(argc, argv);
327          computesky();
306 –        
307 –        if(*(c_perez+1)>0)
308 –        {
309 –          fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));      
310 –          print_error_sky();
311 –          exit(1);
312 –        }
313 –        
328          printsky();
329          return 0;
330 +
331   }
332  
333  
334  
335  
336  
322 –
323 –
324 –
337   void computesky()
338   {
339  
# Line 337 | Line 349 | void computesky()
349          /* compute solar direction */
350                  
351          if (month) {                    /* from date and time */
352 <                int  jd;
341 <                double  sd, st;
352 >                double  sd;
353  
354 <                jd = jdate(month, day);         /* Julian date */
355 <                sd = sdec(jd);                  /* solar declination */
356 <                if (tsolar)                     /* solar time */
357 <                        st = hour;
358 <                else
359 <                        st = hour + stadj(jd);
360 <                                
361 <                
362 <                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
363 <                print_error_sky();
364 <                exit(1);
354 >                st = hour;
355 >                if (year) {                     /* Michalsky algorithm? */
356 >                        double  mjd = mjdate(year, month, day, hour);
357 >                        if (tsolar)
358 >                                sd = msdec(mjd, NULL);
359 >                        else
360 >                                sd = msdec(mjd, &st);
361 >                } else {
362 >                        int  jd = jdate(month, day);    /* Julian date */
363 >                        sd = sdec(jd);                  /* solar declination */
364 >                        if (!tsolar)                    /* get solar time? */
365 >                                st = hour + stadj(jd);
366                  }
367 <                
356 <                                        
367 >                                                        
368                  if(timeinterval) {
369                          
370                          if(timeinterval<0) {
371                          fprintf(stderr, "time interval negative\n");
372                          exit(1);
373                          }
374 <                                        
375 <                        if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) {                
376 <                        
377 <                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
378 <                        st= (st+timeinterval/120+solar_sunrise(month,day))/2;
374 >                                                                        
375 >                        if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) {                      
376 >                         st= (st+timeinterval/120+solar_sunrise(month,day))/2;
377 >                         if(suppress_warnings==0)
378 >                         { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
379                          }
380                  
381 <                        if(fabs(solar_sunset(month,day)-st)<timeinterval/60) {
382 <                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
383 <                        st= (st-timeinterval/120+solar_sunset(month,day))/2;
381 >                        if(fabs(solar_sunset(month,day)-st)<timeinterval/120) {
382 >                         st= (st-timeinterval/120+solar_sunset(month,day))/2;
383 >                         if(suppress_warnings==0)
384 >                         { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
385                          }
386 +                        
387 +                        if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) {
388 +                          if(suppress_warnings==0)
389 +                          { fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); }
390 +                         altitude = salt(sd, st);
391 +                         azimuth = sazi(sd, st);
392 +                         print_error_sky();
393 +                         exit(0);
394 +                        }
395                  }
396 +                else
397                  
398 +                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
399 +                        if(suppress_warnings==0)
400 +                        { fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); }
401 +                        altitude = salt(sd, st);
402 +                        azimuth = sazi(sd, st);
403 +                        print_error_sky();
404 +                        exit(0);
405 +                }
406                  
407                  altitude = salt(sd, st);
408                  azimuth = sazi(sd, st);
# Line 384 | Line 414 | void computesky()
414          
415  
416  
387 –        /* if loop for the -l option. W.Sprenger (01/2013) */
388 –        /*
389 –        if (altitude*180/M_PI < sunaltitude_border) {
390 –        
391 –        if (suppress_warnings==0) {
392 –            fprintf(stderr, "Warning: sun altitude (%.3f degrees) below the border (%.3f degrees)\n",altitude*180/M_PI,sunaltitude_border);
393 –                                   }
394 –        print_error_sky();
395 –        exit(1);
396 –        }
397 –        */
398 –
417                          
418          if (!cloudy && altitude > 87.*M_PI/180.) {
419                  
# Line 407 | Line 425 | void computesky()
425                  altitude = 87.*M_PI/180.;
426          }
427          
428 +        
429 +        
430          sundir[0] = -sin(azimuth)*cos(altitude);
431          sundir[1] = -cos(azimuth)*cos(altitude);
432          sundir[2] = sin(altitude);
# Line 414 | Line 434 | void computesky()
434                  
435          /* calculation for the new functions */
436          sunzenith = 90 - altitude*180/M_PI;
437 <        
418 <        
437 >                        
438  
439          /* compute the inputs for the calculation of the light distribution over the sky*/
440          if (input==0)           /* P */
# Line 439 | Line 458 | void computesky()
458                  check_irradiances();
459                  skybrightness = sky_brightness();
460                  skyclearness =  sky_clearness();
461 +                
462                  check_parametrization();
463 <
463 >                                                        
464                  if (output==0 || output==2)
465                          {
466                          diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
# Line 464 | Line 484 | void computesky()
484                          if (altitude<=0)
485                          {
486                                  if (suppress_warnings==0)
487 <                                     fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n");
487 >                                     fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
488                                  directirradiance = 0;
489                                  diffuseirradiance = 0;
490                          } else {
# Line 504 | Line 524 | void computesky()
524                          {
525                                  if (suppress_warnings==0)
526                                          fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
507 –                                
527                                  globalirradiance=erbs_s0*0.999;                
528                          }
529                          
# Line 517 | Line 536 | void computesky()
536                          directirradiance=globalirradiance-diffuseirradiance;
537                          
538                          printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
539 <                        printf("# WARNING: the -E option is only recommended for a rough estimation!");
539 >                        printf("# WARNING: the -E option is only recommended for a rough estimation!\n");
540                          
541                          directirradiance=directirradiance/sin(altitude);
542                                                                                                                          
# Line 540 | Line 559 | void computesky()
559                  
560                  
561          
562 <        else    {fprintf(stderr,"error at the input arguments"); exit(1);}
562 >        else    { fprintf(stderr,"error at the input arguments"); exit(1); }
563  
564  
565          
# Line 562 | Line 581 | void computesky()
581  
582          
583          
565 –        
566 –        
584          /*calculation of the modelled luminance */
585          for (j=0;j<145;j++)
586          {
# Line 609 | Line 626 | void computesky()
626          else
627                  solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
628          
612 –                        
629  
630  
631 < /* Compute the ground radiance */
632 < zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
633 < zenithbr*=diffnormalization;
631 >        /* Compute the ground radiance */
632 >        zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
633 >        zenithbr*=diffnormalization;
634          
635 < if (skyclearness==1)
635 >        if (skyclearness==1)
636          normfactor = 0.777778;
637                  
638 < if (skyclearness>=6)
638 >        if (skyclearness>=6)
639          {              
640          F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
641          normfactor = normsc()/F2/M_PI;
642          }
643  
644 < if ( (skyclearness>1) && (skyclearness<6) )
644 >        if ( (skyclearness>1) && (skyclearness<6) )
645          {
646          S_INTER=1;
647          F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
648          normfactor = normsc()/F2/M_PI;
649          }
650  
651 < groundbr = zenithbr*normfactor;
651 >        groundbr = zenithbr*normfactor;
652  
653 < if (dosun&&(skyclearness>1))
653 >        if (dosun&&(skyclearness>1))
654          groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
655  
656 < groundbr *= gprefl;
656 >        groundbr *= gprefl;
657  
658  
659 +                
660 +        if(*(c_perez+1)>0)
661 +        {
662 +          if(suppress_warnings==0)
663 +                {  fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));}  
664 +          print_error_sky();
665 +          exit(0);
666 +        }
667  
668 +
669   return;
670   }
671  
672  
673  
674  
650 – void print_error_sky()
651 – {
652 –        sundir[0] = -sin(azimuth)*cos(altitude);
653 –        sundir[1] = -cos(azimuth)*cos(altitude);
654 –        sundir[2] = sin(altitude);
655 –        
656 –        printf("\nvoid brightfunc skyfunc\n");
657 –        printf("2 skybright perezlum.cal\n");
658 –        printf("0\n");
659 –        printf("10 0.00 0.00  0.000 0.000 0.000 0.000 0.000  %f %f %f \n", sundir[0], sundir[1], sundir[2]);
660 – }
661 –        
675  
663 –
664 –
665 –
666 –
667 –
676   double solar_sunset(int month,int day)
677   {
678       float W;
# Line 674 | Line 682 | double solar_sunset(int month,int day)
682   }
683  
684  
685 +
686 +
687   double solar_sunrise(int month,int day)
688   {
689       float W;
# Line 685 | Line 695 | double solar_sunrise(int month,int day)
695  
696  
697  
698 + void printsky()
699 + {      
700 +        
701 +        printf("# Local solar time: %.2f\n", st);
702 +        printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
703  
704  
690 –
691 –
692 –
693 –
694 –
695 – void printsky()                 /* print out sky */
696 – {
705          if (dosun&&(skyclearness>1))
706          {              
707                  printf("\nvoid light solar\n");
# Line 710 | Line 718 | void printsky()                        /* print out sky */
718                  printf("0\n0\n");
719                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
720          }
721 + /* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */      
722 +        if  (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 )  
723 +        {
724 +        fprintf(stderr, "       warning: sky clearness(epsilon)= %f \n",skyclearness);
725 +        fprintf(stderr, "       warning: intermediate sky!! \n");
726 +        fprintf(stderr, "       warning: color model for intermediate sky pending  \n");
727 +        fprintf(stderr, "       warning: no color output ! \n");
728 +        color_output=0;
729 +        }
730 +        if (color_output==1)
731 +        {
732 +        printf("\nvoid colorfunc skyfunc\n");
733 +        printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n");
734 +        printf("0\n");
735 +        printf("22 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", diffnormalization, groundbr,
736 +                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
737 +                sundir[0], sundir[1], sundir[2],skyclearness,locus[0],locus[1],locus[2],locus[3],locus[4],locus[5],locus[6],locus[7],locus[8],locus[9],locus[10]);
738 +        }else{
739 +        printf("\nvoid brightfunc skyfunc\n");
740 +        printf("2 skybright perezlum.cal\n");
741 +        printf("0\n");
742 +        printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
743 +                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
744 +                sundir[0], sundir[1], sundir[2]);
745 +         }
746          
747 + }
748  
749 +
750 +
751 + void print_error_sky()
752 + {
753 +
754 +
755 +        sundir[0] = -sin(azimuth)*cos(altitude);
756 +        sundir[1] = -cos(azimuth)*cos(altitude);
757 +        sundir[2] = sin(altitude);
758 +
759 +        printf("# Local solar time: %.2f\n", st);
760 +        printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
761 +
762          printf("\nvoid brightfunc skyfunc\n");
763          printf("2 skybright perezlum.cal\n");
764          printf("0\n");
765 <        printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
719 <                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
720 <                sundir[0], sundir[1], sundir[2]);
721 <        
765 >        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]);
766   }
767 +        
768  
769  
770 +
771 +
772   void printdefaults()                    /* print default values */
773   {
774          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
# Line 735 | Line 782 | void printdefaults()                   /* print default values */
782   }
783  
784  
785 < void userror(char* msg)                 /* print usage error and quit */
785 >
786 >
787 > void usage_error(char* msg)                     /* print usage error and quit */
788   {
789          if (msg != NULL)
790                  fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
# Line 752 | Line 801 | void userror(char* msg)                        /* print usage error and qui
801          fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
802          fprintf(stderr, "       Output specification with option:\n");
803          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");
804 <        fprintf(stderr, "       gendaylit version 2.3 (2013/08/08)  \n\n");
804 >        fprintf(stderr, "       gendaylit version 2.5 (2018/04/18)  \n\n");
805          exit(1);
806   }
807  
808  
809  
810 +
811   double normsc()           /* compute normalization factor (E0*F2/L0) */
812   {
813          static double  nfc[2][5] = {
# Line 781 | Line 831 | double normsc()                  /* compute normalization fac
831  
832  
833  
834 +
835 +
836   void printhead(int ac, char** av)               /* print command header */
837   {
838          putchar('#');
# Line 794 | Line 846 | void printhead(int ac, char** av)              /* print command he
846  
847  
848  
849 +
850 +
851   /* Perez models */
852  
853   /* Perez global horizontal luminous efficacy model */
# Line 865 | Line 919 | double glob_h_effi_PEREZ()
919  
920  
921  
868 –
922          for (i=1; i<=category_total_number; i++)
923          {
924                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
# Line 879 | Line 932 | double glob_h_effi_PEREZ()
932   }
933  
934  
935 +
936 +
937   /* global horizontal diffuse efficacy model, according to PEREZ */
938   double glob_h_diffuse_effi_PEREZ()
939   {
# Line 886 | Line 941 | double glob_h_diffuse_effi_PEREZ()
941          double  category_bounds[10], a[10], b[10], c[10], d[10];
942          int     category_total_number, category_number, i;
943  
889 –
890 –
891 –        
944          check_parametrization();
945          
946          
947   /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
948 <     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
948 > fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
949      
898 –
899 –
950   /* initialize category bounds (clearness index bounds) */
951  
952          category_total_number = 8;
# Line 952 | Line 1002 | double glob_h_diffuse_effi_PEREZ()
1002  
1003  
1004  
955 –
1005          category_number = -1;
1006          for (i=1; i<=category_total_number; i++)
1007          {
# Line 962 | Line 1011 | double glob_h_diffuse_effi_PEREZ()
1011  
1012          if (category_number == -1) {
1013                  if (suppress_warnings==0)
1014 <                    fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness);
1014 >                    fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
1015                  print_error_sky();
1016 <                exit(1);
1016 >                exit(0);
1017          }
1018                  
1019  
# Line 977 | Line 1026 | double glob_h_diffuse_effi_PEREZ()
1026  
1027  
1028  
1029 +
1030 +
1031 +
1032   /* direct normal efficacy model, according to PEREZ */
1033  
1034   double direct_n_effi_PEREZ()
# Line 987 | Line 1039 | double         category_bounds[10], a[10], b[10], c[10], d[10
1039   int     category_total_number, category_number, i;
1040  
1041  
1042 < if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
1043 <   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");
1042 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
1043 >   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/
1044  
1045  
1046   /* initialize category bounds (clearness index bounds) */
# Line 1062 | Line 1114 | return(value);
1114   /*check the range of epsilon and delta indexes of the perez parametrization*/
1115   void check_parametrization()
1116   {
1117 +
1118   if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1119                  {
1120  
1121   /*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1122 +                
1123                  if (skyclearness<skyclearinf){
1124 <                        if (suppress_warnings==0)
1125 <                            /* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1124 >                        /* if (suppress_warnings==0)
1125 >                            fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1126                          skyclearness=skyclearinf;
1127                  }
1128                  if (skyclearness>skyclearsup){
1129 <                        if (suppress_warnings==0)
1130 <                            /* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1131 <                        skyclearness=skyclearsup-0.1;
1129 >                        /* if (suppress_warnings==0)
1130 >                            fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1131 >                        skyclearness=skyclearsup-0.001;
1132                  }
1133                  if (skybrightness<skybriginf){
1134 <                        if (suppress_warnings==0)
1135 <                            /* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1134 >                        /* if (suppress_warnings==0)
1135 >                            fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1136                          skybrightness=skybriginf;
1137                  }
1138                  if (skybrightness>skybrigsup){
1139 <                        if (suppress_warnings==0)
1140 <                            /* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1139 >                        /* if (suppress_warnings==0)
1140 >                            fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1141                          skybrightness=skybrigsup;
1142                  }
1143  
# Line 1092 | Line 1146 | if (skyclearness<skyclearinf || skyclearness>skyclears
1146   }
1147  
1148  
1149 +
1150 +
1151 +
1152   /* validity of the direct and diffuse components */
1153   void    check_illuminances()
1154   {
1155          if (directilluminance < 0) {
1156 <                fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
1156 >                if(suppress_warnings==0)
1157 >                { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1158                  directilluminance = 0.0;
1159          }
1160          if (diffuseilluminance < 0) {
1161 <                fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n");
1161 >                if(suppress_warnings==0)
1162 >                { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1163                  diffuseilluminance = 0.0;
1164          }
1165 <        if (directilluminance > solar_constant_l*1000.0) {
1166 <                fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n");
1167 <                exit(1);
1165 >        
1166 >        if (directilluminance+diffuseilluminance==0 && altitude > 0) {
1167 >                if(suppress_warnings==0)
1168 >                { fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); }
1169 >                print_error_sky();
1170 >                exit(0);
1171          }
1172 +        
1173 +        if (directilluminance > solar_constant_l) {
1174 +                if(suppress_warnings==0)
1175 +                { fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); }
1176 +                print_error_sky();
1177 +                exit(0);
1178 +        }
1179   }
1180  
1181  
1182   void    check_irradiances()
1183   {
1184          if (directirradiance < 0) {
1185 <                fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n");
1185 >                if(suppress_warnings==0)
1186 >                { fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); }
1187                  directirradiance = 0.0;
1188          }
1189          if (diffuseirradiance < 0) {
1190 <                fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n");
1190 >                if(suppress_warnings==0)
1191 >                { fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); }
1192                  diffuseirradiance = 0.0;
1193          }
1194 +        
1195 +        if (directirradiance+diffuseirradiance==0 && altitude > 0) {
1196 +                if(suppress_warnings==0)
1197 +                { fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); }
1198 +                print_error_sky();
1199 +                exit(0);
1200 +        }
1201 +        
1202          if (directirradiance > solar_constant_e) {
1203 <                fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n");
1204 <                exit(1);
1203 >                if(suppress_warnings==0)
1204 >                { fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); }
1205 >                print_error_sky();
1206 >                exit(0);
1207          }
1208   }
1209          
# Line 1181 | Line 1262 | void illu_to_irra_index()
1262   double  test1=0.1, test2=0.1, d_eff;
1263   int     counter=0;      
1264  
1265 < diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000);
1266 < directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1265 > diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l);
1266 > directirradiance = directilluminance*solar_constant_e/(solar_constant_l);
1267   skyclearness =  sky_clearness();
1268   skybrightness = sky_brightness();
1269   check_parametrization();
# Line 1210 | Line 1291 | while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(d
1291          skybrightness = sky_brightness();
1292          skyclearness =  sky_clearness();
1293          check_parametrization();
1213 –        
1214 –        /*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/
1294                  
1295          }
1296  
# Line 1248 | Line 1327 | double calc_rel_lum_perez(double dzeta,double gamma,do
1327  
1328          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1329          {
1330 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n");
1330 >                fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
1331                  exit(1);
1332          }
1333  
# Line 1299 | Line 1378 | void coeff_lum_perez(double Z, double epsilon, double
1378  
1379          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1380          {
1381 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1381 >                fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n");
1382                  exit(1);
1383          }
1384  
# Line 1368 | Line 1447 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1447          else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
1448          {
1449                  printf("error in calculation of gamma (angle between point and sun");
1450 <                exit(3);
1450 >                exit(1);
1451          }
1452          else
1453                  *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
# Line 1412 | Line 1491 | double         air_mass()
1491   double  m;
1492   if (sunzenith>90)
1493          {
1494 <        fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1495 <        exit(1);
1494 >        if(suppress_warnings==0)
1495 >        { fprintf(stderr, "Warning: air mass has reached the maximal value\n"); }
1496 >        sunzenith=90;
1497          }
1498   m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1499   return(m);

Diff Legend

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