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.12 by greg, Fri Aug 9 16:44:19 2013 UTC vs.
Revision 2.21 by greg, Thu Jan 28 19:03:15 2021 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 + *  version 2.6 (2021/01/29): dew point dependency added according to Perez publication 1990 (W -> atm_preci_water=exp(0.07*Td-0.075) ). by J. Wienold, EPFL
9   */
10  
11 + #define  _USE_MATH_DEFINES
12   #include  <stdio.h>
13   #include  <string.h>
14   #include  <math.h>
# Line 52 | Line 55 | float defangle_phi[] = {
55          90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
56          80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
57          240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};
58 + /* default values for Berlin */
59 + float   locus[] = {
60 + -4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786};
61  
62  
63  
64 < /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
64 > /* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
65   double sky_brightness();
66   double sky_clearness();
67  
# Line 88 | Line 94 | void printdefaults();
94   void check_sun_position();
95   void computesky();
96   void printhead(int ac, char** av);
97 < void userror(char* msg);
97 > void usage_error(char* msg);
98   void printsky();
99  
100   FILE * frlibopen(char* fname);
# Line 97 | Line 103 | FILE * frlibopen(char* fname);
103   double  get_eccentricity();
104   double  air_mass();
105  
106 < double  solar_sunset();
107 < double  solar_sunrise();
102 < double  stadj();
103 < int     jdate(int month, int day);
106 > double  solar_sunset(int month, int day);
107 > double  solar_sunrise(int month, int day);
108  
105 –
106 –
107 – /* sun calculation constants */
108 – extern double   s_latitude;
109 – extern double   s_longitude;
110 – extern double   s_meridian;
111 –
109   const double    AU = 149597890E3;
110   const double    solar_constant_e = 1367;    /* solar constant W/m^2 */
111 < const double    solar_constant_l = 127.5;   /* solar constant klux */
111 > const double    solar_constant_l = 127500;   /* solar constant lux */
112  
113   const double    half_sun_angle = 0.2665;
114   const double    half_direct_angle = 2.85;
115  
116 < const double    skyclearinf = 1.0;      /* limitations for the variation of the Perez parameters */
117 < const double    skyclearsup = 12.1;
116 > const double    skyclearinf = 1.0;          /* limitations for the variation of the Perez parameters */
117 > const double    skyclearsup = 12.01;
118   const double    skybriginf = 0.01;
119   const double    skybrigsup = 0.6;
120  
121  
122  
123   /* required values */
124 + int  year = 0;                                  /* year (optional) */
125   int     month, day;                             /* date */
126   double  hour;                                   /* time */
127   int     tsolar;                                 /* 0=standard, 1=solar */
# Line 136 | Line 134 | double skyclearness = 0;
134   double  skybrightness = 0;
135   double  solarradiance;
136   double  diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;
137 < double  sunzenith, daynumber, atm_preci_water=2;
137 > double  sunzenith, daynumber, atm_preci_water, Td=10.97353115;
138  
139   /*double  sunaltitude_border = 0;*/
140   double  diffnormalization = 0;
# Line 146 | Line 144 | double         *c_perez;
144   int     output=0;       /* define the unit of the output (sky luminance or radiance): */
145                          /* visible watt=0, solar watt=1, lumen=2 */
146   int     input=0;        /* define the input for the calulation */
147 <
147 > int     color_output=0;
148   int     suppress_warnings=0;
149  
150          /* default values */
# Line 157 | Line 155 | double  betaturbidity = 0.1;
155   double  gprefl = 0.2;
156   int     S_INTER=0;
157  
158 +
159          /* computed values */
160   double  sundir[3];
161   double  groundbr = 0;
# Line 168 | Line 167 | float  timeinterval = 0;
167   char    *progname;
168   char    errmsg[128];
169  
170 + double  st;
171  
172  
173 –
173   int main(int argc, char** argv)
174   {
175          int  i;
# Line 181 | Line 180 | int main(int argc, char** argv)
180                  return 0;
181          }
182          if (argc < 4)
183 <                userror("arg count");
183 >                usage_error("arg count");
184          if (!strcmp(argv[1], "-ang")) {
185                  altitude = atof(argv[2]) * (M_PI/180);
186                  azimuth = atof(argv[3]) * (M_PI/180);
# Line 189 | Line 188 | int main(int argc, char** argv)
188          } else {
189                  month = atoi(argv[1]);
190                  if (month < 1 || month > 12)
191 <                        userror("bad month");
191 >                        usage_error("bad month");
192                  day = atoi(argv[2]);
193                  if (day < 1 || day > 31)
194 <                        userror("bad day");
194 >                        usage_error("bad day");
195                  hour = atof(argv[3]);
196                  if (hour < 0 || hour >= 24)
197 <                        userror("bad hour");
197 >                        usage_error("bad hour");
198                  tsolar = argv[3][0] == '+';
199          }
200          for (i = 4; i < argc; i++)
201                  if (argv[i][0] == '-' || argv[i][0] == '+')
202                          switch (argv[i][1]) {
203 +                        case 'd':
204 +                                Td = atof(argv[++i]);
205 +                                if (Td < -40 || Td > 40) {
206 +                                        Td=10.97353115; }
207 +                                break;
208                          case 's':
209                                  cloudy = 0;
210                                  dosun = argv[i][0] == '+';
211                                  break;
212 +                        case 'y':
213 +                                year = atoi(argv[++i]);
214 +                                break;
215                          case 'R':
216                                  u_solar = argv[i][1] == 'R' ? -1 : 1;
217                                  solarbr = atof(argv[++i]);
# Line 213 | Line 220 | int main(int argc, char** argv)
220                                  cloudy = argv[i][0] == '+' ? 2 : 1;
221                                  dosun = 0;
222                                  break;
223 +                        case 'C':
224 +                                if (argv[i][2] == 'I' && argv[i][3] == 'E' ) {
225 +                                locus[0] = -4.607e9;
226 +                                locus[1] = 2.9678e6;
227 +                                locus[2] = 0.09911e3;
228 +                                locus[3] = 0.244063;
229 +                                locus[4] = -2.0064e9;
230 +                                locus[5] = 1.9018e6;
231 +                                locus[6] = 0.24748e3;
232 +                                locus[7] = 0.23704;
233 +                                locus[8] = -3.0;
234 +                                locus[9] = 2.87;
235 +                                locus[10] = -0.275;
236 +                                 }else{ color_output = 1;
237 +                                 }
238 +                                break;
239 +                        case 'l':
240 +                                locus[0] = atof(argv[++i]);
241 +                                locus[1] = atof(argv[++i]);
242 +                                locus[2] = atof(argv[++i]);
243 +                                locus[3] = atof(argv[++i]);
244 +                                locus[4] = locus[0];
245 +                                locus[5] = locus[1];
246 +                                locus[6] = locus[2];
247 +                                locus[7] = locus[3];
248 +                                locus[8] = atof(argv[++i]);
249 +                                locus[9] = atof(argv[++i]);
250 +                                locus[10] = atof(argv[++i]);
251 +                                break;
252 +
253                          case 't':
254                                  betaturbidity = atof(argv[++i]);
255                                  break;
# Line 236 | Line 273 | int main(int argc, char** argv)
273                                  break;
274                          
275                          case 'O':
276 <                                output = atof(argv[++i]);       /*define the unit of the output of the program :
277 <                                                                sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) */
276 >                                output = atof(argv[++i]);       /*define the unit of the output of the program:
277 >                                                                sky and sun luminance/radiance
278 >                                                                (0==W visible, 1==W solar radiation, 2==lm) */
279                                  break;
280                                  
281                          case 'P':
# Line 269 | Line 307 | int main(int argc, char** argv)
307                                  globalirradiance = atof(argv[++i]);
308                                  break;
309                          
272 –                        /*
273 –                        case 'l':
274 –                                sunaltitude_border = atof(argv[++i]);
275 –                                break;
276 –                        */
277 –                                
310                          case 'i':
311                                  timeinterval = atof(argv[++i]);
312                                  break;
# Line 282 | Line 314 | int main(int argc, char** argv)
314                          
315                          default:
316                                  sprintf(errmsg, "unknown option: %s", argv[i]);
317 <                                userror(errmsg);
317 >                                usage_error(errmsg);
318                          }
319                  else
320 <                        userror("bad option");
320 >                        usage_error("bad option");
321  
322 <        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
323 <                fprintf(stderr,
292 <                    "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
322 >        if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180)
323 >                fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n",
324                      progname, (s_longitude-s_meridian)*12/M_PI);
325  
326  
327          /* dynamic memory allocation for the pointers */
297 –        
328          if ( (c_perez = calloc(5, sizeof(double))) == NULL )
329 <        {
300 <                fprintf(stderr,"Out of memory error in function main");
301 <                return 1;
302 <        }
329 >        { fprintf(stderr,"Out of memory error in function main"); return 1; }
330  
331 +        
332 +        atm_preci_water=exp(0.07*Td-0.075);
333          printhead(argc, argv);
334          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 –        
335          printsky();
336          return 0;
337 +
338   }
339  
340  
341  
342  
343  
322 –
323 –
324 –
344   void computesky()
345   {
346  
# Line 337 | Line 356 | void computesky()
356          /* compute solar direction */
357                  
358          if (month) {                    /* from date and time */
359 <                int  jd;
341 <                double  sd, st;
359 >                double  sd;
360  
361 <                jd = jdate(month, day);         /* Julian date */
362 <                sd = sdec(jd);                  /* solar declination */
363 <                if (tsolar)                     /* solar time */
364 <                        st = hour;
365 <                else
366 <                        st = hour + stadj(jd);
367 <                                
368 <                
369 <                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
370 <                print_error_sky();
371 <                exit(1);
361 >                st = hour;
362 >                if (year) {                     /* Michalsky algorithm? */
363 >                        double  mjd = mjdate(year, month, day, hour);
364 >                        if (tsolar)
365 >                                sd = msdec(mjd, NULL);
366 >                        else
367 >                                sd = msdec(mjd, &st);
368 >                } else {
369 >                        int  jd = jdate(month, day);    /* Julian date */
370 >                        sd = sdec(jd);                  /* solar declination */
371 >                        if (!tsolar)                    /* get solar time? */
372 >                                st = hour + stadj(jd);
373                  }
374 <                
356 <                                        
374 >                                                        
375                  if(timeinterval) {
376                          
377                          if(timeinterval<0) {
378                          fprintf(stderr, "time interval negative\n");
379                          exit(1);
380                          }
381 <                                        
382 <                        if(fabs(solar_sunrise(month,day)-st)<timeinterval/60) {                
383 <                        
384 <                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
385 <                        st= (st+timeinterval/120+solar_sunrise(month,day))/2;
381 >                                                                        
382 >                        if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) {                      
383 >                         st= (st+timeinterval/120+solar_sunrise(month,day))/2;
384 >                         if(suppress_warnings==0)
385 >                         { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
386                          }
387                  
388 <                        if(fabs(solar_sunset(month,day)-st)<timeinterval/60) {
389 <                        fprintf(stderr, "Solar position corrected at %d %d %.3f\n",month,day,hour);
390 <                        st= (st-timeinterval/120+solar_sunset(month,day))/2;
388 >                        if(fabs(solar_sunset(month,day)-st)<timeinterval/120) {
389 >                         st= (st-timeinterval/120+solar_sunset(month,day))/2;
390 >                         if(suppress_warnings==0)
391 >                         { fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); }
392                          }
393 +                        
394 +                        if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) {
395 +                          if(suppress_warnings==0)
396 +                          { fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); }
397 +                         altitude = salt(sd, st);
398 +                         azimuth = sazi(sd, st);
399 +                         print_error_sky();
400 +                         exit(0);
401 +                        }
402                  }
403 +                else
404                  
405 +                if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) {
406 +                        if(suppress_warnings==0)
407 +                        { fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); }
408 +                        altitude = salt(sd, st);
409 +                        azimuth = sazi(sd, st);
410 +                        print_error_sky();
411 +                        exit(0);
412 +                }
413                  
414                  altitude = salt(sd, st);
415                  azimuth = sazi(sd, st);
# Line 384 | Line 421 | void computesky()
421          
422  
423  
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 –
424                          
425          if (!cloudy && altitude > 87.*M_PI/180.) {
426                  
# Line 407 | Line 432 | void computesky()
432                  altitude = 87.*M_PI/180.;
433          }
434          
435 +        
436 +        
437          sundir[0] = -sin(azimuth)*cos(altitude);
438          sundir[1] = -cos(azimuth)*cos(altitude);
439          sundir[2] = sin(altitude);
# Line 414 | Line 441 | void computesky()
441                  
442          /* calculation for the new functions */
443          sunzenith = 90 - altitude*180/M_PI;
444 <        
418 <        
444 >                        
445  
446          /* compute the inputs for the calculation of the light distribution over the sky*/
447          if (input==0)           /* P */
# Line 439 | Line 465 | void computesky()
465                  check_irradiances();
466                  skybrightness = sky_brightness();
467                  skyclearness =  sky_clearness();
468 +                
469                  check_parametrization();
470 <
470 >                                                        
471                  if (output==0 || output==2)
472                          {
473                          diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
# Line 464 | Line 491 | void computesky()
491                          if (altitude<=0)
492                          {
493                                  if (suppress_warnings==0)
494 <                                     fprintf(stderr, "Warning: solar zenith angle larger than 90 degrees; using zero irradiance to proceed\n");
494 >                                     fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n");
495                                  directirradiance = 0;
496                                  diffuseirradiance = 0;
497                          } else {
# Line 504 | Line 531 | void computesky()
531                          {
532                                  if (suppress_warnings==0)
533                                          fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n");
507 –                                
534                                  globalirradiance=erbs_s0*0.999;                
535                          }
536                          
# Line 517 | Line 543 | void computesky()
543                          directirradiance=globalirradiance-diffuseirradiance;
544                          
545                          printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance);
546 <                        printf("# WARNING: the -E option is only recommended for a rough estimation!");
546 >                        printf("# WARNING: the -E option is only recommended for a rough estimation!\n");
547                          
548                          directirradiance=directirradiance/sin(altitude);
549                                                                                                                          
# Line 540 | Line 566 | void computesky()
566                  
567                  
568          
569 <        else    {fprintf(stderr,"error at the input arguments"); exit(1);}
569 >        else    { fprintf(stderr,"error at the input arguments"); exit(1); }
570  
571  
572          
# Line 562 | Line 588 | void computesky()
588  
589          
590          
565 –        
566 –        
591          /*calculation of the modelled luminance */
592          for (j=0;j<145;j++)
593          {
# Line 609 | Line 633 | void computesky()
633          else
634                  solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
635          
612 –                        
636  
637  
638 < /* Compute the ground radiance */
639 < zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
640 < zenithbr*=diffnormalization;
638 >        /* Compute the ground radiance */
639 >        zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
640 >        zenithbr*=diffnormalization;
641          
642 < if (skyclearness==1)
642 >        if (skyclearness==1)
643          normfactor = 0.777778;
644                  
645 < if (skyclearness>=6)
645 >        if (skyclearness>=6)
646          {              
647          F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
648          normfactor = normsc()/F2/M_PI;
649          }
650  
651 < if ( (skyclearness>1) && (skyclearness<6) )
651 >        if ( (skyclearness>1) && (skyclearness<6) )
652          {
653          S_INTER=1;
654          F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
655          normfactor = normsc()/F2/M_PI;
656          }
657  
658 < groundbr = zenithbr*normfactor;
658 >        groundbr = zenithbr*normfactor;
659  
660 < if (dosun&&(skyclearness>1))
660 >        if (dosun&&(skyclearness>1))
661          groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];                
662  
663 < groundbr *= gprefl;
663 >        groundbr *= gprefl;
664  
665  
666 +                
667 +        if(*(c_perez+1)>0)
668 +        {
669 +          if(suppress_warnings==0)
670 +                {  fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));}  
671 +          print_error_sky();
672 +          exit(0);
673 +        }
674  
675 +
676   return;
677   }
678  
679  
680  
681  
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 –        
682  
663 –
664 –
665 –
666 –
667 –
683   double solar_sunset(int month,int day)
684   {
685       float W;
# Line 674 | Line 689 | double solar_sunset(int month,int day)
689   }
690  
691  
692 +
693 +
694   double solar_sunrise(int month,int day)
695   {
696       float W;
# Line 685 | Line 702 | double solar_sunrise(int month,int day)
702  
703  
704  
705 + void printsky()
706 + {      
707 +        
708 +        printf("# Local solar time: %.2f\n", st);
709 +        printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
710 +        printf("# epsilon, delta, atmospheric precipitable water content : %.4f %.4f %.4f \n", skyclearness, skybrightness,atm_preci_water );
711  
712  
690 –
691 –
692 –
693 –
694 –
695 – void printsky()                 /* print out sky */
696 – {
713          if (dosun&&(skyclearness>1))
714          {              
715                  printf("\nvoid light solar\n");
# Line 710 | Line 726 | void printsky()                        /* print out sky */
726                  printf("0\n0\n");
727                  printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
728          }
729 + /* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */      
730 +        if  (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 )  
731 +        {
732 +        fprintf(stderr, "       warning: sky clearness(epsilon)= %f \n",skyclearness);
733 +        fprintf(stderr, "       warning: intermediate sky!! \n");
734 +        fprintf(stderr, "       warning: color model for intermediate sky pending  \n");
735 +        fprintf(stderr, "       warning: no color output ! \n");
736 +        color_output=0;
737 +        }
738 +        if (color_output==1)
739 +        {
740 +        printf("\nvoid colorfunc skyfunc\n");
741 +        printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n");
742 +        printf("0\n");
743 +        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,
744 +                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
745 +                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]);
746 +        }else{
747 +        printf("\nvoid brightfunc skyfunc\n");
748 +        printf("2 skybright perezlum.cal\n");
749 +        printf("0\n");
750 +        printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
751 +                *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4),
752 +                sundir[0], sundir[1], sundir[2]);
753 +         }
754          
755 + }
756  
757 +
758 +
759 + void print_error_sky()
760 + {
761 +
762 +
763 +        sundir[0] = -sin(azimuth)*cos(altitude);
764 +        sundir[1] = -cos(azimuth)*cos(altitude);
765 +        sundir[2] = sin(altitude);
766 +
767 +        printf("# Local solar time: %.2f\n", st);
768 +        printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI);
769 +
770          printf("\nvoid brightfunc skyfunc\n");
771          printf("2 skybright perezlum.cal\n");
772          printf("0\n");
773 <        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 <        
773 >        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]);
774   }
775 +        
776  
777  
778 +
779 +
780   void printdefaults()                    /* print default values */
781   {
782          printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
# Line 735 | Line 790 | void printdefaults()                   /* print default values */
790   }
791  
792  
793 < void userror(char* msg)                 /* print usage error and quit */
793 >
794 >
795 > void usage_error(char* msg)                     /* print usage error and quit */
796   {
797          if (msg != NULL)
798                  fprintf(stderr, "%s: Use error - %s\n\n", progname, msg);
799 <        fprintf(stderr, "Usage: %s      month day hour    [...]\n", progname);
800 <        fprintf(stderr, "   or: %s -ang altitude azimuth  [...]\n", progname);
799 >        fprintf(stderr, "Usage: %s      month day hour [-y year]        [...]\n", progname);
800 >        fprintf(stderr, "   or: %s -ang altitude azimuth                [...]\n", progname);
801          fprintf(stderr, "               followed by:      -P          epsilon delta [options]\n");
802          fprintf(stderr, "                        or:      [-W|-L|-G]  direct_value diffuse_value [options]\n");
803 <        fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
804 <        fprintf(stderr, "       Description:\n");
803 >        fprintf(stderr, "                        or:      -E          global_irradiance [options]\n\n");
804 >        fprintf(stderr, "    Description:\n");
805          fprintf(stderr, "       -P epsilon delta  (these are the Perez parameters) \n");
806          fprintf(stderr, "       -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
807          fprintf(stderr, "       -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
# Line 752 | Line 809 | void userror(char* msg)                        /* print usage error and qui
809          fprintf(stderr, "       -E global-horizontal-irradiance (W/m^2)\n\n");
810          fprintf(stderr, "       Output specification with option:\n");
811          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");
812 <        fprintf(stderr, "       gendaylit version 2.3 (2013/08/08)  \n\n");
812 >        fprintf(stderr, "       gendaylit version 2.5 (2018/04/18)  \n\n");
813          exit(1);
814   }
815  
816  
817  
818 +
819   double normsc()           /* compute normalization factor (E0*F2/L0) */
820   {
821          static double  nfc[2][5] = {
# Line 781 | Line 839 | double normsc()                  /* compute normalization fac
839  
840  
841  
842 +
843 +
844   void printhead(int ac, char** av)               /* print command header */
845   {
846          putchar('#');
# Line 794 | Line 854 | void printhead(int ac, char** av)              /* print command he
854  
855  
856  
857 +
858 +
859   /* Perez models */
860  
861   /* Perez global horizontal luminous efficacy model */
# Line 865 | Line 927 | double glob_h_effi_PEREZ()
927  
928  
929  
868 –
930          for (i=1; i<=category_total_number; i++)
931          {
932                  if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
# Line 879 | Line 940 | double glob_h_effi_PEREZ()
940   }
941  
942  
943 +
944 +
945   /* global horizontal diffuse efficacy model, according to PEREZ */
946   double glob_h_diffuse_effi_PEREZ()
947   {
# Line 886 | Line 949 | double glob_h_diffuse_effi_PEREZ()
949          double  category_bounds[10], a[10], b[10], c[10], d[10];
950          int     category_total_number, category_number, i;
951  
889 –
890 –
891 –        
952          check_parametrization();
953          
954          
955   /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
956 <     fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
956 > fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
957      
898 –
899 –
958   /* initialize category bounds (clearness index bounds) */
959  
960          category_total_number = 8;
# Line 952 | Line 1010 | double glob_h_diffuse_effi_PEREZ()
1010  
1011  
1012  
955 –
1013          category_number = -1;
1014          for (i=1; i<=category_total_number; i++)
1015          {
# Line 962 | Line 1019 | double glob_h_diffuse_effi_PEREZ()
1019  
1020          if (category_number == -1) {
1021                  if (suppress_warnings==0)
1022 <                    fprintf(stderr, "ERROR: Model parameters out of range, skyclearness = %lf \n", skyclearness);
1022 >                    fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
1023                  print_error_sky();
1024 <                exit(1);
1024 >                exit(0);
1025          }
1026                  
1027  
# Line 977 | Line 1034 | double glob_h_diffuse_effi_PEREZ()
1034  
1035  
1036  
1037 +
1038 +
1039 +
1040   /* direct normal efficacy model, according to PEREZ */
1041  
1042   double direct_n_effi_PEREZ()
# Line 987 | Line 1047 | double         category_bounds[10], a[10], b[10], c[10], d[10
1047   int     category_total_number, category_number, i;
1048  
1049  
1050 < if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
1051 <   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");
1050 > /*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
1051 >   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/
1052  
1053  
1054   /* initialize category bounds (clearness index bounds) */
# Line 1062 | Line 1122 | return(value);
1122   /*check the range of epsilon and delta indexes of the perez parametrization*/
1123   void check_parametrization()
1124   {
1125 +
1126   if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
1127                  {
1128  
1129   /*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
1130 +                
1131                  if (skyclearness<skyclearinf){
1132 <                        if (suppress_warnings==0)
1133 <                            /* fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1132 >                        /* if (suppress_warnings==0)
1133 >                            fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
1134                          skyclearness=skyclearinf;
1135                  }
1136                  if (skyclearness>skyclearsup){
1137 <                        if (suppress_warnings==0)
1138 <                            /* fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1139 <                        skyclearness=skyclearsup-0.1;
1137 >                        /* if (suppress_warnings==0)
1138 >                            fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
1139 >                        skyclearness=skyclearsup-0.001;
1140                  }
1141                  if (skybrightness<skybriginf){
1142 <                        if (suppress_warnings==0)
1143 <                            /* fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1142 >                        /* if (suppress_warnings==0)
1143 >                            fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
1144                          skybrightness=skybriginf;
1145                  }
1146                  if (skybrightness>skybrigsup){
1147 <                        if (suppress_warnings==0)
1148 <                            /* fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1147 >                        /* if (suppress_warnings==0)
1148 >                            fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
1149                          skybrightness=skybrigsup;
1150                  }
1151  
# Line 1092 | Line 1154 | if (skyclearness<skyclearinf || skyclearness>skyclears
1154   }
1155  
1156  
1157 +
1158 +
1159 +
1160   /* validity of the direct and diffuse components */
1161   void    check_illuminances()
1162   {
1163          if (directilluminance < 0) {
1164 <                fprintf(stderr,"WARNING: direct illuminance < 0. Using 0.0\n");
1164 >                if(suppress_warnings==0)
1165 >                { fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); }
1166                  directilluminance = 0.0;
1167          }
1168          if (diffuseilluminance < 0) {
1169 <                fprintf(stderr,"WARNING: diffuse illuminance < 0. Using 0.0\n");
1169 >                if(suppress_warnings==0)
1170 >                { fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); }
1171                  diffuseilluminance = 0.0;
1172          }
1173 <        if (directilluminance > solar_constant_l*1000.0) {
1174 <                fprintf(stderr,"ERROR: direct illuminance exceeds solar constant\n");
1175 <                exit(1);
1173 >        
1174 >        if (directilluminance+diffuseilluminance==0 && altitude > 0) {
1175 >                if(suppress_warnings==0)
1176 >                { fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); }
1177 >                print_error_sky();
1178 >                exit(0);
1179          }
1180 +        
1181 +        if (directilluminance > solar_constant_l) {
1182 +                if(suppress_warnings==0)
1183 +                { fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); }
1184 +                print_error_sky();
1185 +                exit(0);
1186 +        }
1187   }
1188  
1189  
1190   void    check_irradiances()
1191   {
1192          if (directirradiance < 0) {
1193 <                fprintf(stderr,"WARNING: direct irradiance < 0. Using 0.0\n");
1193 >                if(suppress_warnings==0)
1194 >                { fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); }
1195                  directirradiance = 0.0;
1196          }
1197          if (diffuseirradiance < 0) {
1198 <                fprintf(stderr,"WARNING: diffuse irradiance < 0. Using 0.0\n");
1198 >                if(suppress_warnings==0)
1199 >                { fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); }
1200                  diffuseirradiance = 0.0;
1201          }
1202 +        
1203 +        if (directirradiance+diffuseirradiance==0 && altitude > 0) {
1204 +                if(suppress_warnings==0)
1205 +                { fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); }
1206 +                print_error_sky();
1207 +                exit(0);
1208 +        }
1209 +        
1210          if (directirradiance > solar_constant_e) {
1211 <                fprintf(stderr,"ERROR: direct irradiance exceeds solar constant\n");
1212 <                exit(1);
1211 >                if(suppress_warnings==0)
1212 >                { fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); }
1213 >                print_error_sky();
1214 >                exit(0);
1215          }
1216   }
1217          
# Line 1181 | Line 1270 | void illu_to_irra_index()
1270   double  test1=0.1, test2=0.1, d_eff;
1271   int     counter=0;      
1272  
1273 < diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l*1000);
1274 < directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
1273 > diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l);
1274 > directirradiance = directilluminance*solar_constant_e/(solar_constant_l);
1275   skyclearness =  sky_clearness();
1276   skybrightness = sky_brightness();
1277   check_parametrization();
# Line 1210 | Line 1299 | while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(d
1299          skybrightness = sky_brightness();
1300          skyclearness =  sky_clearness();
1301          check_parametrization();
1213 –        
1214 –        /*fprintf(stderr,"skyclearness = %lf, skybrightness = %lf, directirradiance = %lf, diffuseirradiance = %lf\n",skyclearness, skybrightness, directirradiance, diffuseirradiance);*/
1302                  
1303          }
1304  
# Line 1248 | Line 1335 | double calc_rel_lum_perez(double dzeta,double gamma,do
1335  
1336          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1337          {
1338 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez!\n");
1338 >                fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
1339                  exit(1);
1340          }
1341  
# Line 1299 | Line 1386 | void coeff_lum_perez(double Z, double epsilon, double
1386  
1387          if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
1388          {
1389 <                fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
1389 >                fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n");
1390                  exit(1);
1391          }
1392  
# Line 1348 | Line 1435 | void coeff_lum_perez(double Z, double epsilon, double
1435   /* degrees into radians */
1436   double radians(double degres)
1437   {
1438 <        return degres*M_PI/180.0;
1438 >        return degres*(M_PI/180.);
1439   }
1440  
1441  
1442   /* radian into degrees */
1443   double degres(double radians)
1444   {
1445 <        return radians/M_PI*180.0;
1445 >        return radians*(180./M_PI);
1446   }
1447  
1448  
# Line 1368 | Line 1455 | void theta_phi_to_dzeta_gamma(double theta,double phi,
1455          else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
1456          {
1457                  printf("error in calculation of gamma (angle between point and sun");
1458 <                exit(3);
1458 >                exit(1);
1459          }
1460          else
1461                  *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
# Line 1386 | Line 1473 | double integ_lv(float *lv,float *theta)
1473                  buffer += (*(lv+i))*cos(radians(*(theta+i)));
1474          }
1475                          
1476 <        return buffer*2*M_PI/144;
1476 >        return buffer*(2.*M_PI/145.);
1477   }
1478  
1479  
# Line 1412 | Line 1499 | double         air_mass()
1499   double  m;
1500   if (sunzenith>90)
1501          {
1502 <        fprintf(stderr, "Solar zenith angle larger than 90 degrees in function air_mass()\n");
1503 <        exit(1);
1502 >        if(suppress_warnings==0)
1503 >        { fprintf(stderr, "Warning: air mass has reached the maximal value\n"); }
1504 >        sunzenith=90;
1505          }
1506   m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
1507   return(m);

Diff Legend

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